# Vibration Signal Denoising Using Deep Learning

Yuyan Wu
Department of Civil and Environmental Engineering
Stanford University
<EMAIL_ADDRESS>
Senyang Jiang
Department of Computer Science
Stanford University
<EMAIL_ADDRESS>
Yousef Liang
Department of Computer Science
Stanford University
<EMAIL_ADDRESS>
Yuyan Wu SUNet ID: wuyuyan, Department of Civil and Environmental Engineering,
Stanford University<EMAIL_ADDRESS>Senyang Jiang SUNet ID: senyangj,
Department of Computer Science, Stanford University<EMAIL_ADDRESS>Youzhi Liang SUNet ID: youzhil, Department of Computer Science, Stanford
University<EMAIL_ADDRESS>
##  1 Introduction

Structure vibration signals induced by footsteps are widely used for tasks
like occupant identification [1], localization [2], human activity inference
[3], structure health monitoring [4] and so on. The vibration signals are
collected as time series with amplitude values. However, the collected signals
are always noisy in practice due to the influence of environmental noise,
electromagnetic interference and other factors. The presence of noise affects
the process of signal analysis, thus affecting the accuracy and error of the
final tasks.

###  1.1 Related Works

Currently, there are a lot of previous works related to signal denoising.
These denoising methods are widely used for the denoising of speech signals,
bioacoustic signals and electrocardiogram (ECG) signals. The methods can be
divided into conventional methods and deep learning methods. The conventional
methods include linear filtering, spectrum subtraction and statistical-based
approaches for noise removal such as wiener filtering. Deep learning-based
signal denoising are popular in recent years for its ability to model non-
stationary noise environments. There are many neural network structures
proposed in previous works like RNNoise [5], Demucs [6], DCCRN[7]. Many deep
learning methods are tried for structure vibration signals denoising in
previous research like [8] and [9]. These works demonstrate the superiority of
deep learning for signal-denoising tasks.

Our task is to reduce the noise from the footstep-induced vibration signals.
In this project, we are considering two kinds of noises. One is the stationary
noise due to circuit interference, which includes Gaussian-distribution noise,
white noise and so on. Another is the non-stationary noises such as the
vibrations caused by items dropping, door closing, music noise, etc.

##  2 Dataset

We cannot directly get the vibration signal without the noise because there
always exists electrical noise when collecting data. So we consider using the
simulation method to get a clean signal and add the common distribution of
noises to it for training and testing. To get the noise distribution from item
dropping, we collected the floor vibration signals caused by item dropping to
the ground.

###  2.1 Abaqus Synthetic Dataset

We used Abaqus to build a model to simulate the vibration signal generated by
the floor after applying an impact force. We change the sensing location, the
thickness of the plate, and the force and generated a series of signals. The
model built and the signal simulated is shown in Fig. 1 (a).

Figure 1: (a) The Abaqus model and the generated vibration signal (b) Clean
signal from Abacus vs. noisy signal after applying Gaussian noise (c) An
example of the signals generated based on Eqn. 1, displacement overlaid with a
supplemental Gaussian noise, $w(t)$ as a function of time, $t$ [s].

Using the synthetic clean signal from Abaqus, we can apply noise on top of it
to simulate the noisy signals we would get in real scenarios. For this paper,
we start with Gaussian noises. For each set of the signal, we generate a
Gaussian noise signal of the same length, with a scale selected uniformly
between 0.001 and 0.002.

Each set of signals has nearly 1,000,000 data points, which is too long to be
used as input to the RNN model. As a result, we split each set of signals into
multiple segments of length 500, and use these segments as inputs to the RNN
model. The size of data we have at the end is 112224 segments of length 500.
Fig. 1 (b) is an example of a segment of clean signal and its corresponding
noisy signal after applying the generated Gaussian noise.

###  2.2 PDE/ODE Synthetic Dataset

In light of the high expense of experiments and numerical simulation using
Abacus, we utilized PDE/ODE solver (scipy.integrate) to generate 100k
synthetic time series ($\mathcal{R}^{100,000\times 500})$ with supplementary
noise.

The dynamics of structural vibrations induced by footsteps can be simplified
as the impulse response of a Kirchhoff-Love plate subject to Dirichlet
boundary conditions. Under the assumption of linearity for isotropic plates,
we model the induced vibrations as the solutions of a system of PDEs [10]:

$D_{i}\nabla^{2}\nabla^{2}w_{i}(\mathbf{x},t)-T_{i}\nabla^{2}w_{i}(\mathbf{x},t)=\delta(\mathbf{x},t)-\rho_{i}h_{i}\ddot{w}_{i}(\mathbf{x},t)-K_{i}\dot{w}_{i}(\mathbf{x},t),$
(1)

where $w_{i}(\mathbf{x},t)$ is the transverse deflection,
$D_{i}\sim\mathcal{N}(\mu_{D},\sigma_{D})$,
$T_{i}\sim\mathcal{N}(\mu_{T},\sigma_{T})$,
$\rho_{i}h_{i}\sim\mathcal{N}(\mu_{\rho h},\sigma_{\rho h})$ and
$K_{i}\sim\mathcal{U}(a_{U},b_{U})$ are all parameters of the structure,
governing the dynamical response subject to impulse $\delta(\mathbf{x},t)$.
Each synthetic time series is obtained by $\sum
w_{j}(\mathbf{x},t)+\epsilon_{j}(t)$, where $\epsilon_{j}S(t)$ represents the
supplementary noise, simulating the sensor noise aforementioned, with an
example shown in Fig. 1 (c).

###  2.3 Experimental Dataset

We collected 900 floor vibration signals caused by item dropping to the ground
and 900 floor vibration signals induced by the footstep. The sample of the
signal is shown in the input of Fig 4. This dataset is used for the task of
reducing the noise caused by item dropping.

In addition, we collected 300s’ floor vibration signal of music noise as a
non-stationary noise dataset. This dataset is achieved by placing a speaker on
the floor near the sensor that plays music continuously. We constituted a
dataset for signal separation. The source signal is generated by superimposing
the music noise of a random section with the floor vibration signal induced by
footsteps. The footstep-induced floor vibration signal before adding the music
noise is set as the target signal. This dataset totally contains 3430 samples.

##  3 Tasks, Methods and Results

###  3.1 Task 1: Noise Reduction for Stationary Noises

a. PDE/ODE Synthetic Dataset

We built a ’many-to-many’ LSTM model with a 0.2 drop-out to filter the signal
with Gaussian noise, using a mean square error as the cost function. As shown
in Fig. 2 (a), the loss reduces from 0.37 to 0.16 at the epoch of 100. The
inset shows an example of a filtered and actual signal, suggesting a well-
filtered signal using RNN for the regime of large amplitudes. To further
examine the error, Fig. 2 (b) shows a discrete prediction for larger positive
amplitudes (around 5 to 7.5) and mild heteroskedasticity due to a larger
observed error for small amplitudes (around -2.5 to 2.5).

Figure 2: (a) Training loss versus epoch. (Inset) At epoch = 100, an example
in test set illustrating the unfiltered, filtered, and actual signal (b)
Signal with noise versus actual signal. (Inset) Histogram of the error
distribution. Bi-LSTM Fourier denotes the bi-directional LSTM with additional
frequency-domain input using FFT.

Aiming for higher prediction accuracy of filtered signals, we further explore
feature engineering techniques and alternative RNN architectures, serving as
an alternative to 1D-CNN for biometric signal processing [11]. In light of a
time series representing vibrational signals, features presented in the
frequency domain captures the underlying physics governed by the fundamental
structures. We extract the magnitude and phase using a Fast Fourier Transform
to serve as additional inputs into a bi-directional LSTM neural network (see
Fig. 3). As shown in Fig. 2, the loss vs epoch and error distributions show an
over 50% error reduction using the new architecture, which agrees with the
quantitative metrics, including RMSE, MAE and $R^{2}$-score listed in Table.1.

Figure 3: (a) Architecture for bi-directional LSTM with additional frequency-
domain input. (b) Table for metrics used to evaluate the unfiltered signal and
filtered signal using two LSTM architectures. Bi-LSTM Fourier denotes the bi-
directional LSTM with additional frequency-domain input using FFT.

###  3.2 Task 2: Noise Reduction for Non-Stationary Noises

####  3.2.1 Situation 1: Interference vibration noises do not overlap with the
effective signal (Take Item Dropping as an Example)

To simplify the question, we first assume that floor vibration signal pulses
caused by footstep and item dropping do not overlap with each other. So our
task becomes recognizing the item-dropping-induced vibration signal pulses.
Firstly, we extract features using wavelet transform. This method was chosen
because the wavelet transform has a different time-frequency resolution for
the low-frequency region and the high-frequency region. We build a simple
convolutional neural network as shown in Fig. 4, using Categorical
Crossentropy as the loss function to recognize these two kinds of signals. We
shuffled the dataset and divided it into a training set and a test set in a
9:1 ratio. The recognition result reaches 99%. This shows that our model
performs well in the task of distinguishing between the footstep-induced and
item dropping-induced floor vibration signals.

Figure 4: The CNN structure used for classification of footstep-induced and
item dropping-induced floor vibrations

####  3.2.2 Situation 2: Interference vibration noises overlap with the
effective signal (Take Music Noise as an Example)

Now, we are considering the situation when the non-stationary noises are
overlapped with the effective signals. In this task, we use the dataset with
music noise described in 2.3. The network used here is the U-Net. We use the
wandb package for hyperparameter selection. We select the hyperparameter with
the best result. (num filters = 16, kernel size = 5, batch size = 16, epoch
=100) The structure of the network is described in 5 (a). Finally, the mean
absolute error (MAE) after denoising by U-Net is 0.028. (Before Denoising:
0.06). The denoising result is shown in 5 (c). We can see that the network is
capable of separating the vibration signal from the music noise

Figure 5: (a) The U-Net Structure used for signal separation task. (b)
Validation Set MAE value for different hyperparameter settings. Our final
structure is set according to the green line’s setting (The line with the
minimum validation loss) (c) Signal denoising result

###  3.3 Task 3: General Noise Reduction using RNNoise

As described in related works 1.1, RNNoise [5] is a neural network
architecture used for enhancing human voice from background noises. It
combines classical signal processing with modern deep learning methods. In
this paper, we explore the effectiveness of RNNoise in reducing both
stationary and non-stationary noises, possibly using it as a baseline to
compare with other end-to-end deep learning architectures presented in this
paper.
The source code repo 111Github link: https://github.com/xiph/rnnoise for the
RNNoise paper provides us with a tool for generating the training data. Using
the dataset we created previously, We prepared 27 minutes of clean structural
vibration signals and 30 minutes of Gaussian noises and pure music, and the
tool from RNNoise repo will mix the clean signal and noises to create the
training data.
We test the trained RNNoise model with signals with either Gaussian noises or
Music noises. Fig. 6 shows a segment of the noisy, filtered, and clean signal
on the same diagram.

Figure 6: unfiltered, filtered, actual signal with (a) Gaussian Noise (b)
Music Noise

We also calculated the RMSE (root mean squared error) between the model output
and the ground truth clean signal on Gaussian and Music noises. The RMSE for
Gaussian noise removal is $0.0228$, and the RMSE for music noise removal is
$0.0310$. In comparison, when classical low-pass filters are used for removing
the Gaussian noises and music noises, it can only reach an RMSE value of 0.014
and 0.16 for the Gaussian noises and music noises correspondingly. For
stationary noises like Gaussian noise, the classical low-pass filter has a
better performance than RNNoise. However, for music noise which is non-
stationary, the classical filter performs a lot worse than RNNoise. This shows
that the classical filters are not suitable for filtering non-smooth noise.
RNNoise is more effective in reducing Gaussian noise than music noise. This is
expected since the pattern of music is less predictable than Gaussian noises.
Lastly, we convert the output signal to playable .wav format. By listening to
the output signal, we found Gaussian noises are imperceptible anymore, but we
can still hear the background music, even though its volume has been reduced.

##  4 Summary and Conclusion

In this project, we mainly explore the denoising methods for footstep-induced
vibration signals. We have considered different kinds of noise including
stationary noises such as gaussian noises and non-stationary noises such as
item-dropping vibration noise and music noises.

For the non-stationary noises, we are considering 2 situations. In the first
situation, we are considering some instant noises, taking item-dropping
induced vibration signals as the example. A simple CNN is built to classify
the vibration signals induced by item dropping and footsteps. In this task,
the accuracy we got is 99$\%$. Going further, we are considering the non-
stationary noises that overlaps with the footstep-induced vibration signals,
taking music noise as an example. In this task, we built a U-Net structure and
successfully reduced the music noise from the footstep-induced vibration
signals. The MAE loss drops from 0.06 to 0.028 after denoising using U-Net.

For stationary noises, we further investigate the LSTM architecture using
ODE/PDE synthetic data with Gaussian noise. A LSTM model with dropout reduces
the noise by approximately 70% based on RMSE/MAE metrics. Bi-Directional LSTM
with additional frequency-domain features demonstrates much higher denoising
performance with an $R^{2}$-score over 0.996 and further reduced the
heteroskedasticity over the magnitude of signal.

Finally, we explored using popular denoising architecture (RNNoise) to perform
filtering for both stationary and non-stationary noises. Even though this
model is mainly used for enhancing human voices, through retraining the model
from scratch using our data, we are able to obtain pretty low errors on
removing Gaussian and music noises. By converting the output signal to audio
file, we are able to qualitatively verify that the model is able to filter out
Gaussian noises, and increase signal-to-noise ratio on music noises.

## References

  * [1] Shijia Pan, Tong Yu, Mostafa Mirshekari, Jonathon Fagert, Amelie Bonde, Ole J Mengshoel, Hae Young Noh, and Pei Zhang.  Footprintid: Indoor pedestrian identification through ambient structural vibration sensing.  Proceedings of the ACM on Interactive, Mobile, Wearable and Ubiquitous Technologies, 1(3):1–31, 2017.
  * [2] Mostafa Mirshekari, Shijia Pan, Jonathon Fagert, Eve M Schooler, Pei Zhang, and Hae Young Noh.  Occupant localization using footstep-induced structural vibration.  Mechanical Systems and Signal Processing, 112:77–97, 2018.
  * [3] Shijia Pan, Mostafa Mirshekari, Jonathon Fagert, Ceferino Gabriel Ramirez, Albert Jin Chung, Chih Chi Hu, John Paul Shen, Pei Zhang, and Hae Young Noh.  Characterizing human activity induced impulse and slip-pulse excitations through structural vibration.  Journal of Sound and Vibration, 414:61–80, 2018.
  * [4] Diogo Montalvao, Nuno Manuel Mendes Maia, and António Manuel Relógio Ribeiro.  A review of vibration-based structural health monitoring with special emphasis on composite materials.  Shock and vibration digest, 38(4):295–324, 2006.
  * [5] Jean-Marc Valin.  A hybrid dsp/deep learning approach to real-time full-band speech enhancement.  In 2018 IEEE 20th international workshop on multimedia signal processing (MMSP), pages 1–5. IEEE, 2018.
  * [6] Alexandre Défossez, Nicolas Usunier, Léon Bottou, and Francis Bach.  Music source separation in the waveform domain.  arXiv preprint arXiv:1911.13254, 2019.
  * [7] Yanxin Hu, Yun Liu, Shubo Lv, Mengtao Xing, Shimin Zhang, Yihui Fu, Jian Wu, Bihong Zhang, and Lei Xie.  Dccrn: Deep complex convolution recurrent network for phase-aware speech enhancement.  arXiv preprint arXiv:2008.00264, 2020.
  * [8] Qingsong Xiong, Haibei Xiong, Cheng Yuan, and Qingzhao Kong.  A novel deep convolutional image-denoiser network for structural vibration signal denoising.  Engineering Applications of Artificial Intelligence, 117:105507, 2023.
  * [9] Gao Fan, Jun Li, and Hong Hao.  Vibration signal denoising for structural health monitoring by residual convolutional neural networks.  Measurement, 157:107651, 2020.
  * [10] Duong TA Nguyen, Longfei Li, and Hangjie Ji.  Stable and accurate numerical methods for generalized kirchhoff–love plates.  Journal of Engineering Mathematics, 130(1):1–26, 2021.
  * [11] Zeyu Bai, Ruizhi Yang, and Youzhi Liang.  Mental task classification using electroencephalogram signal.  arXiv preprint arXiv:1910.03023, 2019.


# A Survey of Computational Tools in Solar Physics

Monica G. Bobra  Stuart Mumford  Russell J. Hewett  Steven D. Christe  Kevin
Reardon  Sabrina Savage  Jack Ireland  Tiago M. D. Pereira  Bin Chen  David
Pérez-Suárez  W.W. Hansen Experimental Physics Laboratory, Stanford
University, Stanford, CA 94305, USA  SP2RC, School of Mathematics and
Statistics, The University of Sheffield, Sheffield, UK  Aperio Software Ltd.,
Headingley Enterprise and Arts Centre, Bennett Road, Leeds, UK  Department of
Mathematics, Virginia Polytechnic Institute and State University, Blacksburg,
VA 24061-0123, USA  NASA Goddard Space Flight Center, Greenbelt, MD 20771, USA
National Solar Observatory, Boulder, CO 80303, USA  NASA Marshall Space Flight
Center, Huntsville, AL 35812, USA  Institute of Theoretical Astrophysics,
University of Oslo, Oslo, Norway  Rosseland Centre for Solar Physics,
University of Oslo, Oslo, Norway  Center for Solar-Terrestrial Research, New
Jersey Institute of Technology, Newark, NJ 07102-1982, USA  University College
London, London, UK

###### keywords:

Instrumentation and Data Management

##  1 Introduction

S-Introduction The SunPy Project (The SunPy Community et al., 2020)
facilitates and promotes the use and development of community-led, free, and
open source111According to the Open Source Initiative, stewards of the Open
Source Definition (available at opensource.org/osd), open source software
consists of source code under an open source license. Open source licenses
“allow modifications and derived works, and must allow them to be distributed
under the same terms as the license of the original software.” In addition,
open source software must not discriminate against persons, groups, or fields
and the associated licenses must be non-specific, non-restrictive, and
technology-neutral. data-analysis software for solar physics based on the
scientific Python environment. To better understand the software and hardware
preferences of the solar-physics community, the Project developed a
13-question survey (reproduced in Appendix A-SurveyQs) and disseminated it
internationally222The UK Solar Physics, European Physical Society’s Solar
Physics Division, American Astronomical Society’s Solar Physics Division, the
solar-physics subdivision of the Astronomical Society of Japan, the
Astronomical Society of India, and the Brazilian Astronomical Society
organizations advertised the survey to their members. The SunPy Project also
advertised the survey on the @SunPyProject Twitter account and sent it to the
sunpy and sunpy-dev e-mail lists, both of which are public. over a six-month
period between 7 February 2019 and 28 July 2019.

Many of the survey questions were similar (and in some cases, identical) to
those posed by Momcheva and Tollerud (2015) in an informal survey of 1142
members of the astrophysics community. The SunPy Project did this deliberately
to compare software preferences between the solar and astrophysics
communities.

This article presents the survey results, derived from analyzing 364 responses
from community members across 35 countries. All of the survey responses, along
with the code (Reback et al., 2020; Caswell et al., 2020; Waskom et al., 2020;
van der Walt, Colbert, and Varoquaux, 2011; Bobra, Mumford, and Pereira, 2020)
to analyze these data and produce the figures in this article, are publicly
available at github.com/sunpy/survey.

##  2 Demographics

S-Demographics

Since the SunPy Project relies largely on volunteer efforts, we chose to
construct and disseminate this survey ourselves (instead of going through a
formal channel such as the Statistical Research Center at the American
Institute of Physics). As a result, we recognize that this survey may suffer
from coverage error.

Our survey garnered 368 responses. Most of the survey respondents fit into one
of four career stages: 56% ($n$=205) described themselves as a faculty member,
staff scientist, or researcher, 15% ($n$=53) as a postdoc, 23% ($n$=84) as an
undergraduate or graduate student, and 6% ($n$=22) as a software or instrument
developer. This adds up to $n$=364. Four respondents did not fit into any
career stage, and we dropped their responses from our analysis.

Community members across 35 countries333See the analysis code, available at
github.com/sunpy/survey, for a full list. responded to our survey. About
three-quarters of the respondents came from the US, UK, Germany, India, and
Japan. Together, these five countries include about 1150 solar
physicists444The Solar Physics Division of the American Astronomical Society
includes 521 members (private communication, S. Savage, 28 January 2020). The
UK Solar Physics community estimates “over 150 scientists” on its website,
uksolphys.org, as of 27 January 2020. The European Solar Physics Division
counts 222 members (private communication, T.M.D. Pereira, 28 January 2020).
The Astronomical Society of India includes approximately 100 solar physicists
(private communication, D. Banerjee, 30 January 2020). The communication
newsletter for the solar-physics subdivision of the Astronomical Society of
Japan, called Renraku-kai, counts about 150 subscribers (private
communication, K. Hayashi, 27 January 2020).; therefore, our survey sampled
roughly a quarter of the solar-physics community. Our results are based on the
assumption that our sample is representative of the solar-physics community
overall.

We asked respondents to identify all of the areas of research relevant to
their career. Most respondents identified multiple sub-disciplines of
expertise. We found that 76% ($n$=275) work with space-based observational
data, 46% ($n$=169) work with ground-based observational data, and 26%
($n$=93) work on building instruments. A vast majority of respondents, 82%,
work with ground-based or space-based data. 29% ($n$=105) identified theory as
a relevant sub-discipline, and 47% ($n$=171) identified numerical simulations.

Most of the survey respondents (82%) chose to answer an optional question
about whether they self-identified as an underrepresented minority; 16% of
this subset (13% of the total sample) said yes. 79% of respondents chose to
answer another optional question about whether they self-identified as a
underrepresented gender identity; 11% of this subset (9% of the total sample)
said yes.

##  3 Software Tools

S-Software

In our survey of the solar-physics community, we found that 99$\pm$0.5% of
respondents use software in their research555While three respondents did
apparently indicate that they do not use software in their research, their
further answers on the survey about software package usage suggest that those
might have been erroneous responses.. In a survey of the astrophysics
community, Momcheva and Tollerud (2015) found that 100% of respondents use
software in their research.

We asked users to list all of the scientific software tools, including
programming languages, software development tools, and data-analysis
frameworks, that they utilized within the last year. We summarized their
responses in Figure F-4a. We found that 66% of respondents use the Python
scientific software stack and 73% use IDL666Where relevant, we supplied our
counting error for non-demographic software and hardware related questions
(Questions 6 – 12). For Question 6, we report $\sqrt{3}/364$, or 0.5%, as the
percentage error in the number of no responses. Since this question required
respondents to pick one response from a binary choice, we apply that same
uncertainty to the yes responses. For Questions 7, 8, 10, and 11, which
required respondents to pick only one response from a list of options, we
quantified the percent error in each response simply by applying the square-
root rule for counting experiments (Taylor, 1997). For Questions 9 and 12,
which allowed respondents to select as many options as they liked, we do not
calculate a percent error.. Overall, respondents listed 42 different software
tools and the average respondent used five tools in the past year.

We observe a stark contrast in usage between the two primary data-analysis
languages in solar-physics research, Python and IDL, when viewed by respondent
career stage. The earlier the career stage, the greater the percentage of
Python users: 59% of faculty, staff scientists, and researchers, 75% of
postdocs, and 79% of students use Python. The earlier the career stage, the
fewer IDL users: 78% of faculty, staff scientists, and researchers, 75% of
postdocs, and 60% of students use IDL.

Of course, these tools are not necessarily used in isolation – about half
(45%) of respondents use both Python and IDL. Figure F-4c shows that 28% of
respondents use IDL exclusively (in other words, they use IDL and do not use
Python), while 21% use Python exclusively. The ratio of exclusive IDL users to
exclusive Python users is roughly 2:1 for faculty, staff, and research
scientists and the opposite, 1:2, for students.

Figure 10 of Momcheva and Tollerud (2015) shows that Python is not only the
most popular programming language within their sample of the astrophysics
community, but it is also the most popular within every individual career
category. Our survey results show that Python is the most popular programming
language only among students; IDL and Python are at parity for postdocs, and
IDL is more popular than Python for faculty, staff scientists, researchers,
software developers, and instrument developers. In this respect, the
astrophysics and solar-physics communities differ widely: 78% of solar-physics
faculty, staff scientists, and researchers in our sample use IDL777The use of
IDL by the solar-physics community may be explained partly by how instrument
teams provide their data. Many instrument teams provide data that have been
calibrated to a low level, plus software that allows the data to be further
calibrated for scientific use. The advantage of this model of scientific-data
provision is that as knowledge of the instrument improves over time, the
software can be updated to provide better high-level science-ready data
products. A side-effect of this model of scientific-data provision is that
scientific use of the data requires use of a particular package/language.
Since many instrument teams chose to take advantage of the significant
functionality provided by the SolarSoftWare (SSW: Freeland and Handy 1998)
package, much of the software required to create higher-level data products is
written in the primary language of SSW: IDL. Hence the model of scientific-
data provision may explain why IDL is used by a significant proportion of
respondents., compared with 44% of astrophysics faculty and scientists sampled
by Momcheva and Tollerud (2015).

The two groups of respondents share the same statistics, however, when it
comes to writing software. In both the astrophysics and solar-physics
communities, roughly a third of respondents write their own software most of
the time (see Figure F-3 of this article and Figure 3 of Momcheva and
Tollerud, 2015). Furthermore, about 90% of respondents in both communities
often or occasionally write their own software (see the same figures).

Figure 1.: Summary of results for survey Question 9 “Which of the following
[software tools] have you personally utilized in your work within the last
year?” Results are grouped by self-identified career stage (Question 2).
Respondents listed 42 different software tools; only tools used by 5% or more
of respondents are shown.F-4a

Figure 2.: Comparison of respondents that report using Python or IDL
exclusively by reported career role.F-4c

Figure 3.: Comparison of respondent’s software development and use activities
by reported career role, with uncertainty estimate.F-3

##  4 Education and Training

S-Education Although 99$\pm$0.5% of respondents use software in their research
and 91$\pm$5% often or occasionally write their own software, 63$\pm$4% of
respondents have not had any formal training (e.g. computer-science courses)
at an undergraduate or graduate level. We found that people who write mostly
their own software are no better trained than everyone else: 44$\pm$6% of
people who write their own software reported “a lot (e.g. computer science
courses)” of formal training, compared with 37$\pm$3% overall. We also found
that students today are twice as likely to have a lot of formal training in
programming compared with faculty, researchers, and staff scientists (see
Figure F-2). The amount of training does not vary with area of expertise; each
sub-discipline shows roughly the same amount of formal training as the general
population (37$\pm$3%).

Figure 4.: Comparison of respondent’s formal computer-science education
activities (at both undergraduate and graduate level) by reported career role,
with uncertainty estimate.F-2

##  5 Hardware Tools

S-Hardware We also found that most respondents utilize consumer hardware to
run software for solar-physics research. Although 82% of respondents work with
space-based or ground-based data, and some of these missions (e.g. the Solar
Dynamics Observatory and Daniel K. Inouye Solar Telescope) produce terabytes
of data per day, 14% use a regional or national cluster888We recognize that
some countries, such as the United States, require citizenship or permanent
residence status to use these clusters. and 5% use a commercial cloud provider
(see Figure F-7a). 29% use exclusively a laptop or desktop. The community puts
considerable effort into maintaining clusters and workstations, with 40% of
respondents using a shared workstation, 51% using a local cluster, and 96%
using a laptop or desktop.

These percentages vary significantly by sub-discipline. A larger percentage of
respondents in the numerical simulations and theory sub-disciplines use local
clusters (63% and 60%, respectively, compared with 51% overall) and regional
or national clusters (26% and 26%, respectively, compared with 14% overall).

Figure 5.: Responses to Question 12, related to computer resource and hardware
usage, broken down by career role (Question 2).F-7a

Figure 6.: Responses to Question 12, related to computer resource and hardware
usage, broken down by solar-physics research area (Question 1).F-7b

##  6 Citing Scientific Software

S-Citations Figure F-5 shows that 73$\pm$4% of respondents cite scientific
software in their research, although only 42$\pm$3% do so routinely. Roughly a
quarter (27$\pm$3%) never cite scientific software in their research. When
asked why, about half (53$\pm$8%) responded that they do not know how to
appropriately cite scientific software (see Figure F-6); we note that only
4$\pm$1% of respondents do not think software belongs in citations.

Figure 7.: Responses to Question 10, “Have you cited software papers in your
published research?”F-5

Figure 8.: Responses to Question 11, “Why haven’t you cited software in your
research?”, for those that responded “No” to Question 10.F-6

##  7 Discussion

S-Discussion

Scientific software is an indispensable component of the modern scientific
research workflow (Rüde et al., 2018). Virtually all of the solar-physics
community uses software in their research. Based on this fact, we find three
of the statistics presented in this article worrisome. First, similar to the
astrophysics community999Momcheva and Tollerud (2015) found that only 8$\pm$1%
of the astrophysics community received substantial training; however, their
question did not define “a lot” or “a little”., a significant fraction of the
solar-physics community (63$\pm$4% of respondents) have not taken any
computer-science courses at an undergraduate or graduate level. Second, most
of the solar-physics community (82% of respondents) works with space-based or
ground-based facilities, several of which produce terabyte- or petabyte-sized
data sets, and nearly a third of the community (29% of respondents) uses
exclusively a laptop or desktop to run software for solar-physics research. It
is unclear whether the computing power offered by laptops and desktops limit
the type of scientific endeavors in solar-physics. Finally, less than half of
the community (42$\pm$3% of respondents) routinely cites scientific software
in their research.

The United States National Academies of Sciences and Medicine (2018) report
entitled Software Policy Options for NASA Earth and Space Sciences recognizes
the lack of education in software development among scientists. The report
recommends initiating and sponsoring “programs to educate and train
researchers in open source best practices,” suggesting topics such as “export
controls, licensing and intellectual property, workflows, and software
development.” This includes sponsoring community members to attend conferences
about open source software development, such as Python in Astronomy
(openastronomy.org/pyastro) or Scientific Computing with Python
(conference.scipy.org), take online courses about software development,
available on learning platforms such as Coursera (coursera.org) and edX
(edx.org), join workshops like those led by The Carpentries (carpentries.org),
and develop training programs, such as the Large Synoptic Survey Telescope’s
Data Science Fellowship program (astrodatascience.org). Our findings in
Section S-Education show that the solar physics community could benefit
immensely from education and training in open source software.

The Ford Foundation’s report, entitled Roads and Bridges: The Unseen Labor
Behind Our Digital Infrastructure (Eghbal, 2016), also suggests “expanding the
pool of contributors so that more people, and more types of people, can build
and sustain public software together.” Increasing the diversity of the talent
pool, which is still lacking in the solar-physics community, will help sustain
a long-term future for open source software in solar-physics.

However, maximizing the scientific return of large data sets, such as those
produced by the Solar Dynamics Observatory and the Daniel K. Inouye Solar
Telescope, requires both skill in software development and computational
resources. The United States National Academies of Sciences and Medicine
(2020) report entitled Progress Toward Implementation of the 2013 Decadal
Survey for Solar and Space Physics: A Midterm Assessment and the Kavli
Foundation series of workshops called Petabytes To Science (Bauer et al.,
2019) recommend adopting science platforms, which co-locate both data and
computational resources required to analyze these data. In this paradigm,
users run software in an external computing environment where the data lives,
instead of moving the data to a desktop or laptop where the software lives.
The astrophysics community already developed several science platforms, such
as the ASTRO Data Lab (datalab.noao.edu), run by the NSF’s National Optical-
Infrared Astronomy Research Laboratory. We encourage the solar-physics
community to fund the development of science platforms so that scientists are
not restricted by the computational power of consumer hardware for analyses
involving terabytes of data.

Finally, we recognize that software development, and hardware development,
takes a vast amount of time. This time is rarely recognized by the academic
community, which largely rewards publications. Therefore, we encourage the
community to publish scientific software (by submitting articles that describe
research software to refereed journals and archiving this software in publicly
available digital repositories; see guides.github.com/activities/citable-
code), cite scientific software (see Appendix A-Cite about how to cite
scientific software), and count scientific software as a co-equal research
artifact when considering career evaluation. This has two benefits: it gives
academic credit and career recognition to those who write software and it
makes it easier to reproduce studies in solar-physics.

Some of the earliest advocates for scientific reproducibility, Claerbout and
Karrenbach (1992) and Buckheit and Donoho (1995), suggested that a journal
article “about computational science in a scientific publication is not the
scholarship itself, it is merely advertising of the scholarship.” The actual
scholarship, they argue, is the code and development environment used to
generate the results. Preserving these elements of scholarship require tools
like version control, which create snapshots of software or data as they
change over time. At the moment, less than half the community (44% of
respondents) uses version control101010We found that 44% of respondents
selected the option “Github (or similar)” in Question 9. However, we realize
this option is ambiguous. In retrospect, we should have provided “Git, Github,
or similar” instead of “Github (or similar)” as an option in Question 9.. The
United States National Academies of Sciences and Medicine (2019) report
entitled Reproducibility and Replicability in Science recommends that
“researchers should convey clear, specific, and complete information about any
computational methods and data products that support their published results
in order to enable other researchers to repeat the analysis,” including the
data, study methods, and computational environment.

Scientists make a critical choice when selecting a computational environment,
because the quality of our tools informs the quality of our research. A large
fraction of the community uses the Python scientific-software stack (66% of
respondents). This number will only grow over time, since Python is the most
popular programming language among students in the solar-physics community
(79% of students who took our survey use Python).

There are a number of reasons why the Python scientific-software stack is
growing in prominence both in the solar-physics community and many other
scientific disciplines111111The number of contributors to the SunPy codebase
grew by an average rate of one per month since 2011 (see The SunPy Community
et al., 2020, Figure 1). According to the 2019 Stack Overflow developer
survey, Python is the fastest-growing major programming language today (see
insights.stackoverflow.com/survey/2019); furthermore, most universities use
Python to teach computer science (Guo, 2014).. Interoperability between many
packages for numerical methods, plotting, astronomy, statistics, and computing
(e.g. Virtanen et al., 2020; van der Walt, Colbert, and Varoquaux, 2011;
McKinney, 2010; Hunter, 2007; The Astropy Collaboration et al., 2018;
VanderPlas et al., 2012; Pedregosa et al., 2011; Rocklin, 2015) allows
researchers to write code with relative speed and ease. The rise of more than
fifty packages in heliophysics alone (see heliopython.org) enables
interdisciplinary analysis across traditionally isolated fields. The open-
development model121212An open-development model goes beyond providing open
source software, it also includes making project-level decisions in publicly-
visible and accessible spaces, such as mailing lists, and inviting input from
the user and developer communities (Tollerud et al., 2019)., adopted by most
of the scientific Python ecosystem, improves the longevity of software since
anyone can contribute to the codebase and no single institution or person
controls the software.

For these reasons, the United States National Academies of Sciences and
Medicine (2018) report entitled Software Policy Options for NASA Earth and
Space Sciences recommends that the “NASA Science Mission Directorate should
explicitly recognize the scientific value of open source software and
incentivize its development and support, with the goal that open source
science software becomes routine scientific practice.” As the SunPy Advisory
Board, we endorse this recommendation not only for the NASA Science Mission
Directorate but for scientific funding agencies worldwide.

### Acknowledgments

We would like to thank Ivelina Momcheva, Erik Tollerud, and Nabil Freij for
their help and guidance on this project, the numerous professional astronomy
societies who helped publicize this survey, and everyone who took this survey.

## Disclosure of Potential Conflicts of Interest

The authors, all members of the SunPy Advisory Board serving in a volunteer,
unpaid, capacity, declare that they have no conflicts of interest.

##  Appendix A Survey Questions

A-SurveyQs The full contents of the survey, distributed as a Google Form,
appear below. All the responses to the Question 13, an optional question which
solicited general, free-form comments, are publicly available at
github.com/sunpy/survey.

  1. 1.

Which of these areas of solar physics do you work in? Check all that apply.

     * $\square$

Observational (Space-Based)

     * $\square$

Observational (Ground-Based)

     * $\square$

Numerical Simulations

     * $\square$

Theory

     * $\square$

Instrumentation

  2. 2.

How would you describe the stage of your career?

     * $\ocircle$

Undergraduate student

     * $\ocircle$

Graduate student

     * $\ocircle$

Postdoc

     * $\ocircle$

Faculty, Staff Scientist, Researcher

     * $\ocircle$

Software Developer

     * $\ocircle$

Instrument Developer

     * $\ocircle$

Retired

     * $\ocircle$

My role is something other than solar physics or software development

     * $\ocircle$

Other (Respondents can enter their own description.)

  3. 3.

What country is your institution in? (Respondents check appropriate country
from a list of options.)

  4. 4.

Do you self-identify as one or more underrepresented minorities in solar
physics? This question is optional.

     * $\ocircle$

Yes

     * $\ocircle$

No

  5. 5.

Do you self-identify as a unrepresented gender identity in Solar Physics? This
question is optional.

     * $\ocircle$

Yes

     * $\ocircle$

No

  6. 6.

Do you use software in your research?

     * $\ocircle$

Yes

     * $\ocircle$

No

  7. 7.

Have you had formal training in programming?

     * $\ocircle$

Yes, a lot (e.g. CS courses at an undergraduate or graduate level)

     * $\ocircle$

Yes, a little (e.g. online classes, books, workshops)

     * $\ocircle$

No

  8. 8.

Which of the following statements is most applicable to you?

     * $\ocircle$

I write mostly my own software.

     * $\ocircle$

I mostly use software written by others.

     * $\ocircle$

Somewhere in between.

  9. 9.

Which of the following have you personally utilized in your work within the
last year? Check all that apply.

     * $\square$

IDL

     * $\square$

SolarSoft

     * $\square$

Python

     * $\square$

SunPy

     * $\square$

Shell Scripting

     * $\square$

C

     * $\square$

C++

     * $\square$

Fortran

     * $\square$

IRAF

     * $\square$

Perl

     * $\square$

Javascript

     * $\square$

Julia

     * $\square$

MATLAB

     * $\square$

Java

     * $\square$

R

     * $\square$

SQL

     * $\square$

Ruby

     * $\square$

HTML/CSS

     * $\square$

Spreadsheets (e.g. Excel)

     * $\square$

Mathematica

     * $\square$

MPI

     * $\square$

Github (or similar)

     * $\square$

Other (Respondents can enter their own description.)

  10. 10.

Have you cited software papers in your published research?

     * $\ocircle$

Yes

     * $\ocircle$

Sometimes

     * $\ocircle$

No

  11. 11.

If ‘No’ for the previous question: Why haven’t you cited software in your
research?

     * $\ocircle$

I am not sure how to appropriately cite software

     * $\ocircle$

I do not think it is necessary

     * $\ocircle$

I do not think software belongs in citations

  12. 12.

On which of these have you run software for solar-physics research?

     * $\square$

Laptop / Desktop computer

     * $\square$

Shared workstation

     * $\square$

Local Cluster

     * $\square$

Regional or National Cluster

     * $\square$

GPU

     * $\square$

Commercial cloud

  13. 13.

Do you have any comments? (This is a free form response; comments are not
required. Please feel free to give us feedback about topics like: version
control, collaborative coding platforms such as Github, standard or best
practices in coding, operating systems, text editors, or your personal
experience with writing code and releasing software, or general thoughts about
SunPy).

##  Appendix B Citing Scientific Software

A-Cite To cite scientific software, please follow these two steps:

  1. 1.

Cite the refereed journal article describing the research software. To find
this article, visit the website for a software package and look for citation
instructions. For example, the SunPy website includes dedicated citation
instructions and an associated BibTex entry.

  2. 2.

Cite the software archive. Publicly available digital repositories, such as
Zenodo, issue a Digital Object Identifier (DOI) for archived software. (Some
institutions also provide digital repositories as part of their library
system.) Generally, open source software projects in the Python scientific
stack will archive their software every time that they release a new version.
For example, the SunPy Github page includes the Zenodo DOI for the most recent
release (as of this writing, v1.1.1); clicking on it leads to the Zenodo
deposit, which provides an associated BibTex entry.

Many projects release multiple versions of software per year, but they only
write a refereed journal article once in a while (for example, the SunPy
Project published an article about the v1.0 release, but they will not publish
an article about the v1.1.1 release). Therefore, creating reproducible results
requires citing both the journal article and the software archive. Here is an
example: This research used version 1.1.1 (Mumford et al., 2020) of the SunPy
open source software package (The SunPy Community et al., 2020).

## References

  * Bauer et al. (2019) Bauer, A.E., Bellm, E.C., Bolton, A.S., Chaudhuri, S., Connolly, A.J., Cruz, K.L., Desai, V., Drlica-Wagner, A., Economou, F., Gaffney, N., Kavelaars, J., Kinney, J., Li, T.S., Lundgren, B., Margutti, R., Narayan, G., Nord, B., Norman, D.J., O’Mullane, W., Padhi, S., Peek, J.E.G., Schafer, C., Schwamb, M.E., Smith, A.M., Tollerud, E.J., Weijmans, A.-M., Szalay, A.S.: 2019, Petabytes to Science. arXiv e-prints, arXiv:1905.05116. ADS.
  * Bobra, Mumford, and Pereira (2020) Bobra, M., Mumford, S., Pereira, T.M.D.: 2020, sunpy/survey: survey v0.2.0 (2020-03-09), Zenodo. DOI.
  * Buckheit and Donoho (1995) Buckheit, J., Donoho, D.L.: 1995, Wavelab and reproducible research, Springer-Verlag. https://statweb.stanford.edu/~wavelab.
  * Caswell et al. (2020) Caswell, T.A., Droettboom, M., Lee, A., Hunter, J., Firing, E., Stansby, D., Klymak, J., Hoffmann, T., de Andrade, E.S., Varoquaux, N., Nielsen, J.H., Root, B., Elson, P., May, R., Dale, D., Lee, J.-J., Seppänen, J.K., McDougall, D., Straw, A., Hobson, P., Gohlke, C., Yu, T.S., Ma, E., Vincent, A.F., Silvester, S., Moad, C., Kniazev, N., Ivanov, P., Ernest, E., Katins, J.: 2020, matplotlib/matplotlib v3.1.3, Zenodo. DOI.
  * Claerbout and Karrenbach (1992) Claerbout, J.F., Karrenbach, M.: 1992, Electronic documents give reproducible research a new meaning. In: SEG Technical Program Expanded Abstracts 1992, Society of Exploration Geophysicists, ???, 601\.
  * Eghbal (2016) Eghbal, N.: 2016, Roads and bridges: The unseen labor behind our digital infrastructure. https://www.fordfoundation.org/work/learning/research-reports/roads-and-bridges-the-unseen-labor-behind-our-digital-infrastructure.
  * Freeland and Handy (1998) Freeland, S., Handy, B.N.: 1998, Data analysis with the SolarSoft system. Solar Phys. 182, 497\. DOI.
  * Guo (2014) Guo, P.: 2014, Python is now the most popular introductory teaching language at top u.s. universities. https://cacm.acm.org/blogs/blog-cacm/176450.
  * Hunter (2007) Hunter, J.D.: 2007, Matplotlib: A 2d graphics environment. Computing in Science & Engineering 9, 90\. DOI.
  * McKinney (2010) McKinney, W.: 2010, Data structures for statistical computing in python. In: van der Walt, S., Millman, J. (eds.) Proc. 9th Python in Science Conf., 51\.
  * Momcheva and Tollerud (2015) Momcheva, I., Tollerud, E.: 2015, Software Use in Astronomy: An Informal Survey, arXiv:1507.03989. ADS.
  * Mumford et al. (2020) Mumford, S.J., Christe, S., Freij, N., Mayer, F., Hughitt, K., Ryan, D.F., Liedtke, S., Shih, A.Y., Pérez-Suárez, D., Chakraborty, P., I, V.K., Inglis, A.R., Pattnaik, P., Sipőcz, B.M., Sharma, R., Leonard, D., Hewett, R.J., Alex-Ian-Hamilton, Stansby, D., Panda, A., Earnshaw, M., Choudhary, N., Kumar, A., Hayes, L., Chanda, P., Haque, M.A., Konge, S., mdmueller, Kirk, M., haathi: 2020, Sunpy, Zenodo. DOI.
  * National Academies of Sciences and Medicine (2018) National Academies of Sciences, E., Medicine: 2018, Open source software policy options for nasa earth and space sciences, The National Academies Press, Washington, DC. ISBN 978-0-309-48271-4. DOI.
  * National Academies of Sciences and Medicine (2019) National Academies of Sciences, E., Medicine: 2019, Reproducibility and replicability in science, The National Academies Press, Washington, DC. ISBN 978-0-309-48616-3. DOI.
  * National Academies of Sciences and Medicine (2020) National Academies of Sciences, E., Medicine: 2020, Progress toward implementation of the 2013 decadal survey for solar and space physics: A midterm assessment, The National Academies Press, Washington, DC. ISBN 978-0-309-67127-9. DOI.
  * Pedregosa et al. (2011) Pedregosa, F., Varoquaux, G., Gramfort, A., Michel, V., Thirion, B., Grisel, O., Blondel, M., Prettenhofer, P., Weiss, R., Dubourg, V., Vanderplas, J., Passos, A., Cournapeau, D., Brucher, M., Perrot, M., Duchesnay, E.: 2011, Scikit-learn: Machine learning in Python. Journal of Machine Learning Research 12, 2825\.
  * Reback et al. (2020) Reback, J., McKinney, W., jbrockmendel, den Bossche, J.V., Augspurger, T., Cloud, P., gfyoung, Sinhrks, Klein, A., Roeschke, M., Tratner, J., She, C., Hawkins, S., Ayd, W., Petersen, T., Schendel, J., Hayden, A., Garcia, M., MomIsBestFriend, Jancauskas, V., Battiston, P., Seabold, S., chris-b1, h-vetinari, Hoyer, S., Overmeire, W., alimcmaster1, Mehyar, M., Dong, K., Whelan, C.: 2020, pandas-dev/pandas: Pandas 1.0.1, Zenodo. DOI.
  * Rocklin (2015) Rocklin, M.: 2015, Dask: Parallel Computation with Blocked algorithms and Task Scheduling. In: Huff, K., Bergstra, J. (eds.) Proc. of the 14th Python in Science Conf., 126\. DOI.
  * Rüde et al. (2018) Rüde, U., Willcox, K., McInnes, L.C., Sterck, H.D.: 2018, Research and education in computational science and engineering. SIAM Review 60, 707\. DOI. https://doi.org/10.1137/16M1096840.
  * Taylor (1997) Taylor, J.: 1997, An introduction to error analysis: the study of uncertainties in physical measurements, University Science Books, Sausalito, CA.
  * The Astropy Collaboration et al. (2018) The Astropy Collaboration, Price-Whelan, A.M., Sipőcz, B.M., Günther, H.M., Lim, P.L., Crawford, S.M., Conseil, S., Shupe, D.L., Craig, M.W., Dencheva, N., Ginsburg, A., VanderPlas, J.T., Bradley, L.D., Pérez-Suárez, D., de Val-Borro, M., Aldcroft, T.L., Cruz, K.L., Robitaille, T.P., Tollerud, E.J., Ardelean, C., Babej, T., Bach, Y.P., Bachetti, M., Bakanov, A.V., Bamford, S.P., Barentsen, G., Barmby, P., Baumbach, A., Berry, K.L., Biscani, F., Boquien, M., Bostroem, K.A., Bouma, L.G., Brammer, G.B., Bray, E.M., Breytenbach, H., Buddelmeijer, H., Burke, D.J., Calderone, G., Rodríguez, J.L.C., Cara, M., Cardoso, J.V.M., Cheedella, S., Copin, Y., Corrales, L., Crichton, D., D’Avella, D., Deil, C., Depagne, É., Dietrich, J.P., Donath, A., Droettboom, M., Earl, N., Erben, T., Fabbro, S., Ferreira, L.A., Finethy, T., Fox, R.T., Garrison, L.H., Gibbons, S.L.J., Goldstein, D.A., Gommers, R., Greco, J.P., Greenfield, P., Groener, A.M., Grollier, F., Hagen, A., Hirst, P., Homeier, D., Horton, A.J., Hosseinzadeh, G., Hu, L., Hunkeler, J.S., Ivezić, Ž., Jain, A., Jenness, T., Kanarek, G., Kendrew, S., Kern, N.S., Kerzendorf, W.E., Khvalko, A., King, J., Kirkby, D., Kulkarni, A.M., Kumar, A., Lee, A., Lenz, D., Littlefair, S.P., Ma, Z., Macleod, D.M., Mastropietro, M., McCully, C., Montagnac, S., Morris, B.M., Mueller, M., Mumford, S.J., Muna, D., Murphy, N.A., Nelson, S., Nguyen, G.H., Ninan, J.P., Nšthe, M., Ogaz, S., Oh, S., Parejko, J.K., Parley, N., Pascual, S., Patil, R., Patil, A.A., Plunkett, A.L., Prochaska, J.X., Rastogi, T., Janga, V.R., Sabater, J., Sakurikar, P., Seifert, M., Sherbert, L.E., Sherwood-Taylor, H., Shih, A.Y., Sick, J., Silbiger, M.T., Singanamalla, S., Singer, L.P., Sladen, P.H., Sooley, K.A., Sornarajah, S., Streicher, O., Teuben, P., Thomas, S.W., Tremblay, G.R., Turner, J.E.H., Terrón, V., van Kerkwijk, M.H., de la Vega, A., Watkins, L.L., Weaver, B.A., Whitmore, J.B., Woillez, J., Zabalza, V.: 2018, The astropy project: Building an open-science project and status of the v2.0 core package. The Astronomical Journal 156, 123\. DOI.
  * The SunPy Community et al. (2020) The SunPy Community, Barnes, W.T., Bobra, M.G., Christe, S.D., Freij, N., Hayes, L.A., Ireland, J., Mumford, S., Perez-Suarez, D., Ryan, D.F., Shih, A.Y., Chanda, P., Glogowski, K., Hewett, R., Hughitt, V.K., Hill, A., Hiware, K., Inglis, A., Kirk, M.S.F., Konge, S., Mason, J.P., Maloney, S.A., Murray, S.A., Panda, A., Park, J., Pereira, T.M.D., Reardon, K., Savage, S., Sipőcz, B.M., Stansby, D., Jain, Y., Taylor, G., Yadav, T., Rajul, Dang, T.K.: 2020, The SunPy Project: Open Source Development and Status of the Version 1.0 Core Package. Astrophys. J. 890, 68\. DOI. ADS.
  * Tollerud et al. (2019) Tollerud, E., Smith, A., Price-Whelan, A., Cruz, K., Norman, D., Narayan, G., Mumford, S., Allen, A., Chan, C.-k., Cherinka, B., Drlica-Wagner, A., Foreman-Mackey, D., Ginsburg, A., Gradvhol, A., Harrington, J., Hogg, D., Jartaltepe, J., Kinney, J., Merchant, N., Momcheva, I., Murphy, N., Peek, J., Peeples, M.S., Pickering, T., Rodriguez, D., Shamir, L., Sinha, M., Sipőcz, B., Sobeck, J., Sosey, M., Stevance, H., Teuben, P., Vohl, D., Weiner, B., Aldcroft, T., Allen, A., Alpaslan, M., Anderson, L., Barentsen, G., Bektesevic, D., Benavides, J., Berriman, B., Blanton, M., Bosch, J., Bouquin, D., Bradley, L., Bryan, G., Burke, D., Burns, K., Buzasi, D., Cabral, J.B., Cardoso, J.V.d.M., Chen, B., Clarkson, W., Collins, M., Corrales, L., Craig, M., Crawford, S., Domagal-Goldman, S., Dong, C., Durbin, M., Faherty, J.K., Farr, W., Forschini, L., Golkhou, V.Z., Günther, H.M., Hafok, H., Hahn, C., Hathi, N., Hedges, C., Huang, S., Hummels, C., Hunt, E., Huppenkothen, D., Juneau, S., van Kerkwijk, M., Kerzendorf, W., Laginja, I., Law, C., de Leon, J., Li, T., Lim, P.L., Malz, A.I., Mao, Y.-Y., Melchior, P., Merin, B., Miller, B., Modjaz, M., Morton, T., Mullally, S., Ogando, R., Parejko, J.K., Paz, D., Pearson, S., Pontoppidan, K., Pope, B., Rapetti, D., Rawls, M., Read, J., Robitaille, T., Rudnick, G., Sharma, S., Sharma, S., Shupe, D., Speagle, J., Starkenburg, T., Stasyszyn, F., Streicher, O., Tremblay, G., Villaescusa-Navarro, F., Vos, J.M., Weaver, B.A., Weltman, A., Wetzel, A., Williams, P.K.G., Winkel, B.: 2019, Sustaining Community-Driven Software for Astronomy in the 2020s. In: BAAS 51, 180\. ADS.
  * van der Walt, Colbert, and Varoquaux (2011) van der Walt, S., Colbert, S.C., Varoquaux, G.: 2011, The numpy array: A structure for efficient numerical computation. Computing in Science Engineering 13, 22\. DOI.
  * VanderPlas et al. (2012) VanderPlas, J., Connolly, A.J., Ivezic, Z., Gray, A.: 2012, Introduction to astroML: Machine learning for astrophysics. In: Proceedings of Conference on Intelligent Data Understanding (CIDU), Boulder, CO, 47\. IEEE. DOI. ADS.
  * Virtanen et al. (2020) Virtanen, P., Gommers, R., Oliphant, T.E., Haberland, M., Reddy, T., Cournapeau, D., Burovski, E., Peterson, P., Weckesser, W., Bright, J., van der Walt, S.J., Brett, M., Wilson, J., Jarrod Millman, K., Mayorov, N., Nelson, A.R.J., Jones, E., Kern, R., Larson, E., Carey, C., Polat, İ., Feng, Y., Moore, E.W., Vand erPlas, J., Laxalde, D., Perktold, J., Cimrman, R., Henriksen, I., Quintero, E.A., Harris, C.R., Archibald, A.M., Ribeiro, A.H., Pedregosa, F., van Mulbregt, P., SciPy 1.0 Contributors: 2020, SciPy 1.0: Fundamental Algorithms for Scientific Computing in Python. Nature Methods 17, 261\. DOI. ADS.
  * Waskom et al. (2020) Waskom, M., Botvinnik, O., Ostblom, J., Lukauskas, S., Hobson, P., MaozGelbart, Gemperline, D.C., Augspurger, T., Halchenko, Y., Cole, J.B., Warmenhoven, J., de Ruiter, J., Pye, C., Hoyer, S., Vanderplas, J., Villalba, S., Kunter, G., Quintero, E., Bachant, P., Martin, M., Meyer, K., Swain, C., Miles, A., Brunner, T., O’Kane, D., Yarkoni, T., Williams, M.L., Evans, C.: 2020, mwaskom/seaborn: v0.10.0 (january 2020), Zenodo. DOI. https://doi.org/10.5281/zenodo.3629446.


# Determining the parameters of the spiral arms of the Galaxy from kinematic
tracers based on Gaia DR3 data

S. I. Denyshchenko,1 P. N. Fedorov,1 V. S. Akhmetov,1,2 A. B. Velichko1,3 and
A. M. Dmytrenko1
1Institute of astronomy of V. N. Karazin Kharkiv national university, Svobody
sq. 4, 61022, Kharkiv, Ukraine
2INAF-Osservatorio Astrofisico di Torino, Via Osservatorio 20, Pino Torinese,
Turin, I-10025, Italy
3Department of Astronomy, University of Geneva, Chemin Pegasi 51, 1290
Versoix, Switzerland<EMAIL_ADDRESS><EMAIL_ADDRESS><EMAIL_ADDRESS>(VSA)

(Accepted XXX. Received YYY; in original form ZZZ)

###### Abstract

We present the results of determining the parameters of the spiral arms of the
Galaxy using the stars Gaia DR3, whose absolute magnitude is $M_{G}<$ 4, and
which allow tracing spiral arms at large distances from the Sun. As tracers of
spiral arms, we use the centroids of stellar spherical regions with a radius
of 0.5 kpc, in which the deformation velocities along the coordinate axis $R$
are insignificant. These kinematic tracers cover the Galactic plane within the
Galactocentric coordinate ranges $140^{\circ}<\theta<220^{\circ}$ and 4 kpc
$<R<$ 14 kpc. The numerical values of the pitch angles of the spirals and
their Galactocentric distances to the point of intersection of the spiral with
the direction of the Galactic center - the Sun are in good agreement with the
results of other authors. By extrapolating beyond the data we have, we present
a schematic four-arm global pattern, consisting of the Scutum–Centaurus,
Sagitarius-Carina, Perseus, Norma-Outer arms, as well as the local arm Orion.
The uncertainties of the determined spiral parameters confirm that the
structures identified are not false, but are reliable from the statistical
point of view.

###### keywords:

methods: data analysis – Galaxy: structure – Galaxy: kinematics and dynamics –
stars: kinematics and dynamics

††pubyear: 2023††pagerange: Determining the parameters of the spiral arms of
the Galaxy from kinematic tracers based on Gaia DR3 data–LABEL:lastpage

##  1 Introduction

The nature of the formation of the spiral structure of the galaxy is still not
determined. There are many assumptions about the causes of arms and spiral
patterns in the Galaxy (Toomre, 1981; Sellwood & Binney, 2002; Masset &
Tagger, 1997; Merrifield et al., 2006; Patsis, 2006; Quillen et al., 2011)

According to the classical wave theory, proposed by Lin & Shu (1964), spiral
arms are waves of increased density rotating around the center of the Galaxy
as a solid body at a constant angular velocity (i.e., pattern velocity),
despite the differential rotation of stars and interstellar medium, that
results in long lifetime of spiral arms (see also Shu, 2016). To date, many
works have been devoted to the study of the spiral structure, using various
methods and objects (Lin & Shu, 1964; Lin et al., 1969; Georgelin & Georgelin,
1976; Taylor & Cordes, 1993; Russeil, 2003; Paladini et al., 2004; Popova &
Loktin, 2005; Dias & Lépine, 2005; Levine et al., 2006; Moitinho et al., 2006;
Vázquez et al., 2008; Gerhard, 2011; Efremov, 2011; Bobylev & Bajkova, 2014;
Hou & Han, 2014, 2015; Dambis et al., 2015; Reid et al., 2019; Xu et al.,
2018; Xu et al., 2021; Poggio et al., 2021; Hao et al., 2021). Most of these
methods are based on the analysis of trigonometric parallaxes and radial
velocities of young objects such as OB stars, open clusters, hydrogen clouds,
HII regions, giant molecular clouds, methanol masers, etc. However, the spiral
structures of the Galaxy contain not only young objects, but also most other
types of objects, including red giants and subgiants. For example, Junqueira
et al. (2015) used a sample of giants of different ages for a new method for
determining the velocity of the Milky Way spiral pattern based on the
interaction between spiral arms and disk stars. The main problems in the
studies given above are the difficulty in determining trigonometric
parallaxes, their significant errors, as well as the limited number of samples
of the objects under study.

The third release of the Gaia space mission, which contains high-precision
proper motions, parallaxes, radial velocities and other astrometric and
astrophysical parameters for stars of various types, makes it possible to
analyze in great details the structural features of an extended region of our
Galaxy. A large sample of objects can be used not only to analyze the spatial
structure of the Galaxy from positions of the objects, but also to perform a
detailed kinematic analysis of Galactic various regions, that makes it
possible to obtain additional information about the structural features of
these regions. A frequently used physical model for studying stellar
kinematics is the Ogorodnikov-Milne (O–M) model (Ogorodnikov, 1965; Clube,
1972; Du Mont, 1977; Miyamoto & Soma, 1993; Miyamoto & Zhu, 1998). This model
makes it possible to study the velocity field in a deformable stellar system.
Fedorov et al. (2021, 2023) proposed a method for obtaining kinematic
parameters within the framework of the O–M model for stellar systems whose
centroids are located in the Galactic plane. In these works, the region of
determination of the kinematic parameters of stellar systems has been extended
to heliocentric distances up to 10 kpc.

One of the main goals of the spiral structure theory, which is based on
observational data, is to explain the stability of spiral patterns. Why do
these patterns exist throughout the many revolutions of the Galaxy, while the
differential galactic rotation tends to destroy these patterns? Obviously, if
the spiral pattern existed for a short period of time, we would not observe
them in a large number of disk galaxies, since differential rotation deforms
any structure in the disk, causing it to be destroyed during one or two
rotations of the disk. One of the possible explanations for the stability of
spiral patterns is the absence of contraction-expansion deformations that
destroy their shape in those regions of the Galaxy where spiral arms are
located.

From a kinematic point of view, this means that in the field of stellar
velocities described by the O–M model, the diagonal components of the
deformation velocity tensor should be close to zero. In this case, the
deformation velocities of the stellar system along the coordinate axes
(relative contraction-expansion velocities) will be small on average or will
be completely absent, which probably affects the stability of the stellar
structure configuration. The results of determining the parameters of spiral
arms given in this work are based on the assumptions mentioned above and are
verified by comparing them with the results obtained in various works by other
authors.

##  2 Kinematic parameters and data used

In our previous papers (Fedorov et al., 2021, 2023; Dmytrenko et al., 2023),
the kinematic parameters of the O–M model for stellar systems contained in
spherical regions of a 1 kpc radius, whose centroids are located in the
Galactic plane, have been determined in local coordinate systems. The
orientation of the local coordinate system with the origin at an arbitrary
point of the Galactic plane corresponds to the following conditions: the
$x^{\prime}$ axis is directed to the center of the Galaxy, the $y^{\prime}$
axis is in the direction of the Galaxy rotation, and the $z^{\prime}$ axis is
perpendicular to the Galactic plane. As shown earlier (Dmytrenko et al.,
2023), in the local coordinate system, the components $M^{+}_{13}$ and
$M^{+}_{23}$ of the velocity deformation tensors for these stellar systems are
close to zero over the entire range of changes in the Galactocentric
cylindrical coordinates $R$ and $\theta$, which indicates the motion of stars,
mainly in the Galactic plane. Therefore, in this work, we consider the tensor
$M^{+}$ as flat, i.e., having only four components
$M^{+}_{11},M^{+}_{12},M^{+}_{21},M^{+}_{22}$, of which only three are
independent ($M^{+}_{12}=M^{+}_{21}$).

Figure 1: The O–M model parameter $M^{+}_{11}$ taken from Fedorov et al.
(2023) as a function of the Galactic coordinates.

Figure 2: The O–M model parameter $M^{+}_{22}$ taken from Fedorov et al.
(2023) as a function of the Galactic coordinates.

Figures 1, 2 show the diagonal components of the deformation velocity tensor
$M^{+}_{11}$ and $M^{+}_{22}$, which have been obtained in our previous work
(Fedorov et al., 2023). As can be seen from Fig. 2, in the range of angles
$\theta$ from 210 to 150 degrees and distances $R$ from 4 to 12 kpc, the
magnitude of change in velocity $V_{\theta}$ along the unit vector j
(parameter $M_{22}$), that characterizes contraction - expansion along the
$y^{\prime}$ axis of the local coordinate system is approximately $\pm$2.5$\rm
km\,s^{-1}\,kpc^{-1}$. Compared to the average velocity $V_{\theta}$ (about
240 $\rm km\,s^{-1}$), this is approximately from 0.5 to 12% depending on $R$.
At the same time, in the same region of the Galaxy, the $M^{+}_{11}$
parameter, which contraction - expansion along the $x^{\prime}$ axis of the
system (or along the unit vector i), reaches values of about $\pm$10 $\rm
km\,s^{-1}\,kpc^{-1}$ inside some ring-shaped structures. Compared to the
radial velocity $V_{R}$ (about $\pm$10 $\rm km\,s^{-1}\,kpc^{-1}$), its change
$\partial V_{R}/\partial R$ = $M^{+}_{11}$ is almost 100% depending on $R$ and
is decisive in the formation of velocity field deformations. In this regard,
for the region of the Galaxy under study, we neglect the influence of the
velocity deformation along the $y^{\prime}$ axis of the local coordinate
system on the stability of the stellar structure configuration.

Also, Fig. 1 clearly shows that next to these ring-shaped structures in the
Galactic plane there are regions adjacent, inside which the gradient $\partial
V_{R}/\partial R$ = $M^{+}_{11}$ is close to zero. For example, this is a
ring-shaped structure located at a Galactocentric distance of about 13 kpc.
Similar structures, although less clear, can be seen at distances of about 9
kpc, 6 kpc, and 4 kpc. This means that the velocities of the deformation
velocities along the $R$ direction inside these regions are about zero, as a
result of which there is no expansion or contraction of stellar systems along
the $R$ coordinate. Then in the case when its average value is very close to
zero, there should be no deformations of the structures observed in these
regions, and as a result, their shape should not change with time. The
property of the invariability of the shape of structures containing objects
with the kinematic parameter $\partial V_{R}/\partial R\approx 0$ allows us to
assume that there is a relation between the stability of the spiral pattern in
our Galaxy and the parameter $\partial V_{R}/\partial R\approx 0$.

###  2.1 Sample

To test the assumption mentioned above, we created a special sample of stars
from the third release of the Gaia data – Gaia DR3 (Gaia Collaboration et al.,
2016, 2023), that has been used to determine the kinematic parameters of the
O–M model. The required input data of the model are positions $\alpha,\delta$,
spatial velocities (proper motions $\mu_{\alpha}\,{\rm
cos}\,\delta,\mu_{\delta}$ and radial velocities Vr) of the stars, as well as
their distances (parallaxes $\pi$).

Gaia DR3 provides estimates of 5 astrometric parameters (positions, proper
motions and parallaxes) for 1.8 billion stars, as well as photometry in G, GRP
and GBP bands (Riello et al., 2021). Among them, radial velocities (the sixth
astrometric parameter) have been measured for $\sim$33 million objects. The
presence of 6 astrometric parameters makes it possible to obtain complete
information about stellar motions in space. Of these 33 million objects, we
selected high-luminosity stars with absolute magnitude of M${}_{G}<$ 4 (see
Fig. 3) as the base sample, see details in Akhmetov et al. (2023).

Figure 3: The color (GBP \- GRP) - absolute magnitude MG diagram based on the
Gaia DR3 data.

To provide the best quality of astrometric parameters, we excluded from this
sample the stars for which the following conditions are satisfied:

$\begin{cases}RUWE>1.4,\\\ \pi/\sigma_{\pi}<5,\\\ \end{cases}$ (1)

About 19.3 million stars meet the criteria listed above and we use them for
the kinematic analysis. Since, as is known from Lindegren et al. (2021), the
Gaia DR3 parallaxes are systematically biased relative to the expected
distribution around zero by several tens of microarcseconds, we corrected
them. To correct parallaxes, we used the Parallax bias Z5 computed according
to Table 9 given by Lindegren et al. (2021). Our sample is limited to about
magnitude 17 due to the inclusion of radial velocities. Therefore, the
magnitude range of the sample has been divided into only two parts. In the
range of G magnitudes from 6 to 13, we applied an offset of -30 $\mu$as, and
in the range of G$>$13, we used a value of -40 $\mu$as. Since in Fig. 20 given
by Lindegren et al. (2021) jumps are observed in the range from 11.5 to 13 G
magnitudes, we decided to use an average value of Z5 approximately equal to
-30 $\mu$as in this range.

In addition, we corrected for the proper motions of stars in the magnitude
range from 6 to 13 as well. Brandt (2018); Lindegren et al. (2018) have showed
that in the second release of Gaia (DR2) data, the reference frame of bright
stars rotates relative to faint stars and quasars at a velocity of $\sim$0.15
$\rm mas\,yr^{-1}$. In the early third Gaia release EDR3, this rotation has
been previously excluded (see section 4.5 by Lindegren et al., 2021). However,
an additional correction has been proposed by Cantat-Gaudin & Brandt (2021),
and we use it to align the stellar proper motions in our sample brighter than
G = 13 with the International Celestial Coordinate System.

###  2.2 Distribution of centroids of selected stellar systems in the Galactic
plane

To solve the system of equations of the O–M model, we divided the final sample
into stellar systems, which are regions inside spheres with a radius of
$R_{\rm s}$ = 0.5 kpc. We set the centers of these spheres at the nodes of a
rectangular grid given in the Galactic plane. The grid nodes are separated
from each other along the $x$ and $y$ coordinates at a distance of 100 pc. We
chose the coverage region of the Galactic plane by the rectangular grid of
nodes (centroids) within $140^{\circ}<\theta<220^{\circ}$, 1.5 kpc $<R<$ 13.5
kpc of Galactocentric (GC) cylindrical coordinates with distance form Sun to
GC is $R_{0}=8.15$ kpc(see, for instance, Fig. 1, 2). We solve the system of
equations of the O–M model for 12 unknowns by the least squares method (LSM)
using the astrometric data of stars, the number of which in each sphere is at
least 1000. Figs. 4, 5 and 6 show the dependencies of
$M^{+}_{11},M^{+}_{22},M^{+}_{33}$ on the $X$ and $Y$ coordinates in the range
of angles $\theta$ from 210 to 150 degrees and distances $R$ from 4 to 12 kpc.

Figure 4: The O–M model parameter $M^{+}_{11}$.

Figure 5: The O–M model parameter $M^{+}_{22}$.

Figure 6: The O–M model parameter $M^{+}_{33}$.

After calculating the kinematic parameters inside of each sphere, we create a
subsample of centroids of those spherical regions in which the value of the
parameter $M^{+}_{11}$ turned out to be insignificant, i.e. the calculated
value of the parameter $M^{+}_{11}$ does not exceed twice the value of its
uncertainty. As mentioned above, we do not consider the influence of the
parameters $M^{+}_{22}$ and $M^{+}_{33}$ on stability, but focus only on
establishing a relations between the parameters $M^{+}_{11}$, which
characterize the contraction-expansion of the velocity field in spherical
regions and spiral arms. The numerical value of the radius of the spherical
region equal to 0.5 kpc has been chosen from empirical considerations and is a
compromise between ensuring the reliability of the LSM solution (a sufficient
number of stars within the sphere) and detailing of the mapped kinematic
parameters.

Figure 7: The distribution of the centroids coordinates in the Galactic $XY$
plane, which are formed from stars contained in spheres with a radius of
$R_{\rm s}$ = 0.5 kpc and having insignificant $M^{+}_{11}$.

Figure 8: Same as in Fig. 7 but with superimposed spirals taken from Reid et
al. (2019).

Fig. 7 shows the distribution of the centroid coordinates in the Galactic $XY$
plane, which are formed from stars contained in spheres with a radius of
$R_{\rm s}$ = 0.5 kpc and having insignificant $M^{+}_{11}$. To understand
whether Fig. 7 contains the centroids belonging to some spiral arms, we
superimposed the spiral arms on the centroid distribution with the parameters
obtained by Reid et al. (2019), as shown in Fig. 8. As can be seen from the
latter figure, individual sections of the location of the centroids on the
coordinate plane are in good agreement with the superimposed spirals. In this
regard, we believe that we can consider the positions of centroids with
insignificant $M^{+}_{11}$ as kinematic tracers of spiral arms in the range of
Galactocentric cylindrical coordinates $140^{\circ}<\theta<220^{\circ}$ , 4
kpc $<R<$ 13 kpc . The choice of these coordinate ranges is mainly due to the
accuracy of the parallaxes used and the accuracy of the obtained kinematic
parameters. The total number of centroids from the specified range was 25,110.
It is these centroids that we use further in the work as kinematic tracers of
spiral arms.

##  3 Determination and comparison of spiral arm parameters

To calculate the parameters of the Milky Way spirals, we apply a method that
is widely used to study the spiral structure of the Galaxy from the data of
various objects (Popova & Loktin, 2005; Xu et al., 2013; Bobylev & Bajkova,
2014; Reid et al., 2019; Hao et al., 2021; Hou, 2021). The method is based on
the "position angle - logarithm of the distance" diagram. Considering the
centroid as a fictitious star belonging to the spiral arm, the equation that
determines the position of this fictitious star on the logarithmic spiral can
be written as follows:

$R=a_{0}\,{\rm exp}\\{(\varphi-\varphi_{0}){\rm tg}\,i\\}$ (2)

where $R$ is the Galactocentric distance of the centroid,
$a_{0}=R(\varphi=\varphi_{0}=0)$ is the Galactocentric distance to the point
of intersection of the spiral with the direction of the Galactic center - the
Sun, $\varphi$ is the position angle of the centroid, determined from the
formula ${\rm tg}\varphi=y/(R_{0}-x)$, the value $R_{0}$ is taken equal to
8.15 kpc, $x,y$ are heliocentric Galactic rectangular coordinates of the
centroid, with the $x$ axis directed from the Sun to the Galactic center, and
the direction of the $y$ axis coincides with that of rotation of the Galaxy;
$\varphi_{0}$ is some arbitrarily chosen initial angle, which, for example, is
taken equal to the average value of all centroid azimuths belonging to the
arm. In this paper, we take the angle $\varphi_{0}$ equal to zero. The pitch
angle of the spiral pattern is denoted as $i$ ($i<0$ for a leading spiral).

According to our sample, we calculate the heliocentric distance $r$ to a
particular centroid from the parallaxes corrected by us, as the average
distances to stars included in each spherical region with a radius of 0.5 kpc.
We calculate the Galactocentric distance $R$ to a particular centroid and the
corresponding position angle $\varphi$ using the following formulas

$\displaystyle R^{2}=(R_{0}-x)^{2}+y^{2}$ (3) $\displaystyle{\rm
tg}\,\varphi=y/(R_{0}-x)$ (4)

where $x,y$ are the heliocentric Galactic rectangular coordinates of the
centroid. Obviously, the relations between the position angle $\varphi$ and
the coordinate angle $\theta$ is $\varphi+\theta=180^{\circ}$.

Since we have chosen $\varphi_{0}=0$, equation 2 can be rewritten as:

$ln(R/a_{0})=\varphi\,{\rm tg}\,i$ (5)

or in the form:

$ln(R)=\varphi\,{\rm tg}\,i+ln(a_{0}).$ (6)

The latter can be rewritten as:

$ln(R)=k\varphi\,+b$ (7)

where $k={\rm tg}\,i$, and $b=ln(a_{0})$.

As can be seen, relation 7 is the straight line equation on the plane
"position angle - logarithm of distance". The solution of the system of
conditional equations separately for each spiral arm by the LSM gives two
quantities: $k$ and $b$. Obviously, $i={\rm arctg}k$, where $i$ is the pitch
angle of the spiral, and $a_{0}=e^{b}$ and represents the Galactocentric
distance of the point of intersection of the spiral with the $X$ axis directed
from the Galactic center and passing through the Sun.

The key point in the approach we use is to find those points (centroids) that
belong to a particular spiral arm. Obviously, the choice of such centroids for
constructing segments of straight lines on the plane "position angle $\varphi$
\- logarithm of the distance ln($R$)" is a rather subjective procedure (see
Fig. 7). We objectified the procedure for selecting points (centroids), which
we use later to construct segments of straight lines, by analyzing the
behavior of the centroids‘ radial velocities of our sample depending on the
Galactocentric coordinates $R$ and $\theta$.

Figure 9: Dependence of the radial velocities $V_{R}$ of the centroids of the
original sample as a function of the Galactocentric distance $R$ at the fixed
angle $\theta=180^{\circ}$. Red color shows the centroids for which the
condition that the $M^{+}_{11}$ parameter is not significant has been met.

Figure 10: Same as in Fig. 9 but at other fixed coordinate angles:
$\theta=160^{\circ}$ (top left), $170^{\circ}$ (top right), $190^{\circ}$
(bottom left) and $200^{\circ}$ (bottom right).

For example, Fig. 9 shows the radial velocity of the centroids of the initial
sample depending on $R$ at the fixed angle $\theta=180^{\circ}$. Fig. 10 shows
the same dependencies but at other coordinate angles symmetrically inclined
with respect to $\theta=180^{\circ}$. We built such dependencies of the radial
velocities on $R$ with a step of 1 degree in $\theta$ for the entire part of
the Galaxy studied in this work. The red color on these dependencies shows the
centroids for which the condition that the $M^{+}_{11}$ parameter is not
significant has been met. As expected, the red dots in the figures are grouped
in places of local extrema, as well as in those places where the radial
velocity $V_{R}$ is almost constant.

It is also clearly seen from the figures that when the angle $\theta$ changes,
the colored dots move along the $R$ axis in a certain interval $\Delta\,R$. In
this case, the average value of the distances $R_{m}$ from the intervals
$\Delta\,R$ turned out to correlate with the Galactocentric distances to the
spiral arms $a_{0}$ known from the literature. These empirical facts allowed
us to form an algorithm for the selection of centroids to construct linear
dependencies ${\rm ln}\,R=k\,\varphi+b$.

To form a primary list of centroids belonging to a particular arm, we use the
following conditions:

1\. In a spherical region with a radius of 0.5 kpc, the $M^{+}_{11}$ parameter
is insignificant: $|M^{+}_{11}/\sigma_{M^{+}_{11}}|<2$.

2\. Using the method described above for first approximation values of $a_{0}$
the mean of Galactocentric distance $R_{m}$ of the spiral arms turned out to
be equal 5.5, 6.8, 8.76, 10.4 and 12.3 kpc for the Scutum-Centaurus arm, the
Sagitarius-Carina arm, the Orion local arm, the Perseus arm and the Norma-
Outer arm respectively.

3\. From the analysis of various works (Dambis et al., 2015; Reid et al.,
2019; Veselova & Nikiforov, 2020; Hao et al., 2021; Hou, 2021, and other), in
the first approximation, the pitch angle of all the spiral arms $i$ has been
taken equal to -13 degrees.

4\. For each spiral arm the values $R_{i}$ were calculated using the above
initial values of $a_{0}$, $i$ and position angle $\varphi_{c}$ of the
centroids by the formula 2.

5\. We selected the centroids of the spherical regions with the Galactocentric
distance $R_{c}$ which are in the interval $\pm$300 pc from the value obtained
in the previous step $R_{i}$ for a particular spiral arm, i.e
$|R_{i}-R_{c}|<300$ pc as we can see on Fig.11.

We use the $R_{c}$ and $\varphi_{c}$ coordinates of centroids of spherical
regions satisfying the conditions listed above to solve the equations 7 by the
LSM with iterative excluding centroids according to the criterion: deviation
of an individual point from the solution exceeds $2\sigma$.

The final results of determining the parameters of the spiral arms $i={\rm
arctg}k$ and $a_{0}=e^{b}$ , as well as their uncertainties, are shown in
Table 1 and plotted in Fig. 12.

Table 1: Parameters of spiral arms. I - Scutum–Centaurus arm, II - Sagitarius-Carina arm, III - Perseus arm, IV - Outer arm, Local - Orion arm. N is the number of centroids used. | I | II | Local | III | IV
---|---|---|---|---|---
$i$ | -12.04$\pm$0.34 | -12.07$\pm$0.17 | -12.43$\pm$0.32 | -12.07$\pm$0.18 | -12.43$\pm$0.18
$a_{0}$ | 5.493$\pm$0.014 | 6.878$\pm$0.010 | 8.719$\pm$0.0.022 | 10.470$\pm$0.012 | 12.289$\pm$0.009
N | 151 | 202 | 122 | 215 | 371
Table 2: Parameters of spiral arms. (1) - Bobylev & Bajkova (2014), (2) - Reid et al. (2019), (3) - Veselova & Nikiforov (2020), (4) - Hou (2021), (5) - Xu et al. (2018). I - Scutum–Centaurus arm, II - Sagitarius-Carina arm, III - Perseus arm, IV - Outer arm, Local - Orion arm. | I | II | Local | III | IV | Ref
---|---|---|---|---|---|---
$i$ | -11.2$\pm$4.0 | -9.3$\pm$2.2 | -10.2$\pm$0.3 | -14.8$\pm$0.8 | -11.5$\pm$1.9 | 1
| -13.1$\pm$2.0 | -9.0$\pm$1.9 | -11.4$\pm$1.9 | -9.5$\pm$2.0 | -6.2$\pm$4.2 | 2
| -11.7$\pm$0.9 | -13.1$\pm$1.4 | -9.9$\pm$1.2 | -6.2$\pm$1.6 | -5.2$\pm$2.8 | 3
| -8 | -13.8 | -8.5 | -11.9 |  | 4
| -18.7$\pm$0.8 | -13.5$\pm$0.5 | -11.5$\pm$0.5 | -9.0$\pm$0.1 |  | 5
|  |  |  |  |  |
$a_{0}$ | 4.5$\pm$0.2 | 6.8$\pm$0.3 | 8.1$\pm$0.3 | 9.9$\pm$0.4 | 13.5$\pm$0.5 | 1
| 4.9$\pm$0.1 | 6.0$\pm$0.1 | 8.3$\pm$0.1 | 8.9$\pm$0.1 | 12.2$\pm$0.4 | 2
| 6.07$\pm$0.04 | 6.78$\pm$0.05 | 8.19$\pm$0.05 | 9.74$\pm$0.09 | 12.02$\pm$0.17 | 3
| 6.5 | 6.9 | 8.5 | 9.6 |  | 4
| 5.9$\pm$0.1 | 7.2$\pm$0.1 | 8.3$\pm$0.1 | 10.6$\pm$0.1 |  | 5

The derived parameters of spiral arms at the level of accuracy of their
determination turned out to be in good agreement with the parameters given in
other studies using various tracers (see Table 2).

Below we present plots for linear dependencies on the “position angle -
logarithm of distance” plane and plots of logarithmic spirals constructed in
the Galactic plane according to the data in Table 1. Specific dependencies are
color-coded. Violet color corresponds to the Scutum-Centaurus arm, green – to
the Sagitarius-Carina arm, red – to the Orion local arm, blue – to the Perseus
arm, and black – to the Norma-Outer arm. All figures are built according to
the algorithm described above.

Figure 11: Dependence of ln($R$) on the position angle $\varphi$. Blue dots
and line (–) corresponds to the Perseus arm, green dots and line (–) - to the
Sagitarius-Carina arm, black dots and line (–) - to the Outer arm, violet dots
and line (–) - to the Scutum-Centaurus arm, red dots and line (–) - to the
Orion local arm.

Figure 12: Logarithmic spirals plotted using the corresponding parameters
derived in this work and given in Table 1. Color points show the centroids
belonging to corresponding spiral arms. The gray ellipse shows the Galactic
bar reproduced according to the data obtained by Wegg et al. (2015).

##  4 Results and discussion

Figure 13: Schematic representation of the the spiral pattern of the Milky Way
in the rectangular Galactocentric coordinate system. Blue curve (–)
corresponds to the Perseus arm, green curve (–) - to the Sagitarius-Carina
arm, black curve (–) - to the Outer arm, violet curve (–) - to the Scutum-
Centaurus arm, red curve (–) - to the Orion local arm. The yellow point shows
the position of the Sun. The gray ellipse shows the Galactic bar reproduced
according to the data obtained by Wegg et al. (2015).

Fig. 13 shows the four main spiral arms of the Milky Way in the rectangular
Galactocentric coordinate system. They are built according to the parameters
determined in this work by extrapolation beyond the data we have and
schematically represent the spiral pattern of our Galaxy (Grand design). Also,
Fig. 13 shows the Galactic bar reproduced according to the data obtained by
Wegg et al. (2015), where the major and minor axes of the bar were determined
to be 5 and 1.5 kpc, respectively (see Rattenbury et al., 2007; Cao et al.,
2013, as well). Its major axis inclination to the X-axis is approximately 30∘.
The bar is shown in Fig. 13 schematically as a filled gray ellipse.

Note that, unlike other works, we do not use the concept of "kink" in an arm,
introduced by Reid et al. (2019), but assumed that any arm as a whole can be
characterized by one value of the pitch angle of. It is also worth noting that
we could not determine the width of the spiral arms in this work. However, it
can be recalled that the used tracers of the spiral arms are maximally
separated by 600 pc in the Galactocentric distance, which can be considered as
a lower estimate of their width.

As can be seen from the figure, the spiral arms intersect with the major axis
of the bar at four points. The intersection points of the Norma-Outer and
Scutum-Centaurus arms with the major axis of the bar are located in the second
Galactocentric quadrant Q2 and have approximate coordinates $x_{\rm N}=2.5$,
$y_{\rm N}=1.0$; $x_{\rm Sc}=4.0$, $y_{\rm Sc}=2.5$ kpc, respectively. The
Norma-Outer arm starts at the indicated point and continues into the quadrants
Q3 and Q4, going counterclockwise into the quadrant Q1. The Scutum-Centaurus
arm originates approximately near the point $x_{\rm Sc}=4.0$, $y_{\rm Sc}=2.5$
kpc in the quadrant Q2 and extends counterclockwise into the Q3 quadrant and
then into the quadrant Q1.

The Sagitarius-Carina and Perseus arms start near the far end of the bar and
have approximate coordinates $x_{\rm S}=-3.0$, $y_{\rm S}=-1.5$ kpc; $x_{\rm
P}=-4.0$, $y_{\rm P}=-2.5$ kpc, respectively. The Sagitarius-Carina arm from
the quadrant Q4 passes counterclockwise through the quadrants Q1, Q2 and Q3.
The Perseus arm starts at Q4, passes through Q1, Q2 and Q3, and extends
further counterclockwise.

The behavior of spiral arms described above practically coincides with the
descriptions of their behavior given in various works to which we refer.

##  5 Conclusions

The parameters of the spiral structure obtained in this work indicate their
good agreement with the results obtained by other authors using various
indicators of the spiral structure and different approaches for their
determination. The use in this work of statistically insignificant kinematic
parameters $M^{+}_{11}$ for the selection of centroids, the positions of which
are considered as kinematic tracers of spiral arms, turned out to be a fairly
reliable method.

Note that although the number of centroids used to determine the spiral
parameters is a few hundred, the number of stars in spherical regions reaches
hundreds of thousands. This allows us to determine the kinematic parameters,
as well as the positions of the centroids with high accuracy, which saves us
from arbitrariness when choosing tracers for fixing specific spiral arms.

The proposed approach made it possible to find stable structures in the
Galaxy, in particular, to reveal spiral arms and determine their parameters.
In addition, the method used makes it possible to obtain information not only
about spiral arms, but also some additional information. For example, as Figs.
7 and 11 show, stable structures ($M^{+}_{11}$ \- insignificant) exist not
only in the form of spirals, but also in the form of some formations extending
in the space between them. The interpretation of these facts requires further
accurate study, which is beyond the scope of this article.

## Data availability

The used catalogue data is available in a standardised format for readers via
the CDS (https://cds.u-strasbg.fr). The software code used in this paper can
be made available upon request by emailing the corresponding author.

##  6 Acknowledgements

This work has made use of data from the European Space Agency (ESA) mission
Gaia (https://www.cosmos.esa.int/gaia), processed by the Gaia Data Processing
and Analysis Consortium (DPAC,
https://www.cosmos.esa.int/web/gaia/dpac/consortium). Funding for the DPAC has
been provided by national institutions, in particular the institutions
participating in the Gaia Multilateral Agreement.

We are immensely grateful to the Armed Forces of Ukraine for the fact that in
wartime we still have the opportunity to work and do science.

## References

  * Akhmetov et al. (2023) Akhmetov V. S., Bucciarelli B., Crosta M., Lattanzi M. G., Spagna A., Re Fiorentin P., Bannikova E. Y., submitted in 2023, MNRAS,
  * Bobylev & Bajkova (2014) Bobylev V. V., Bajkova A. T., 2014, MNRAS, 437, 1549
  * Brandt (2018) Brandt T. D., 2018, ApJS, 239, 31
  * Cantat-Gaudin & Brandt (2021) Cantat-Gaudin T., Brandt T. D., 2021, A&A, 649, A124
  * Cao et al. (2013) Cao L., Mao S., Nataf D., Rattenbury N. J., Gould A., 2013, MNRAS, 434, 595
  * Clube (1972) Clube S. V. M., 1972, MNRAS, 159, 289
  * Dambis et al. (2015) Dambis A. K., et al., 2015, Astronomy Letters, 41, 489
  * Dias & Lépine (2005) Dias W. S., Lépine J. R. D., 2005, ApJ, 629, 825
  * Dmytrenko et al. (2023) Dmytrenko A. M., Fedorov P. N., Akhmetov V. S., Velichko A. B., Denyshchenko S. I., 2023, MNRAS, 521, 4247
  * Du Mont (1977) Du Mont B., 1977, A&A, 61, 127
  * Efremov (2011) Efremov Y. N., 2011, Astronomy Reports, 55, 108
  * Fedorov et al. (2021) Fedorov P. N., Akhmetov V. S., Velichko A. B., Dmytrenko A. M., Denischenko S. I., 2021, MNRAS, 508, 3055
  * Fedorov et al. (2023) Fedorov P. N., Akhmetov V. S., Velichko A. B., Dmytrenko A. M., Denyshchenko S. I., 2023, MNRAS, 518, 2761
  * Gaia Collaboration et al. (2016) Gaia Collaboration et al., 2016, A&A, 595, A1
  * Gaia Collaboration et al. (2023) Gaia Collaboration et al., 2023, A&A, 674, A1
  * Georgelin & Georgelin (1976) Georgelin Y. M., Georgelin Y. P., 1976, A&A, 49, 57
  * Gerhard (2011) Gerhard O., 2011, Memorie della Societa Astronomica Italiana Supplementi, 18, 185
  * Hao et al. (2021) Hao C. J., et al., 2021, A&A, 652, A102
  * Hou (2021) Hou L. G., 2021, Frontiers in Astronomy and Space Sciences, 8, 103
  * Hou & Han (2014) Hou L. G., Han J. L., 2014, A&A, 569, A125
  * Hou & Han (2015) Hou L. G., Han J. L., 2015, MNRAS, 454, 626
  * Junqueira et al. (2015) Junqueira T. C., Chiappini C., Lépine J. R. D., Minchev I., Santiago B. X., 2015, MNRAS, 449, 2336
  * Levine et al. (2006) Levine E. S., Blitz L., Heiles C., 2006, Science, 312, 1773
  * Lin & Shu (1964) Lin C. C., Shu F. H., 1964, ApJ, 140, 646
  * Lin et al. (1969) Lin C. C., Yuan C., Shu F. H., 1969, ApJ, 155, 721
  * Lindegren et al. (2018) Lindegren L., et al., 2018, A&A, 616, A2
  * Lindegren et al. (2021) Lindegren L., et al., 2021, A&A, 649, A4
  * Masset & Tagger (1997) Masset F., Tagger M., 1997, A&A, 322, 442
  * Merrifield et al. (2006) Merrifield M. R., Rand R. J., Meidt S. E., 2006, MNRAS, 366, L17
  * Miyamoto & Soma (1993) Miyamoto M., Soma M., 1993, AJ, 105, 691
  * Miyamoto & Zhu (1998) Miyamoto M., Zhu Z., 1998, AJ, 115, 1483
  * Moitinho et al. (2006) Moitinho A., Vázquez R. A., Carraro G., Baume G., Giorgi E. E., Lyra W., 2006, MNRAS, 368, L77
  * Ogorodnikov (1965) Ogorodnikov K. F., 1965, Dynamics of stellar systems
  * Paladini et al. (2004) Paladini R., Davies R., De Zotti G., 2004, Ap&SS, 289, 363
  * Patsis (2006) Patsis P. A., 2006, MNRAS, 369, L56
  * Poggio et al. (2021) Poggio E., et al., 2021, A&A, 651, A104
  * Popova & Loktin (2005) Popova M. E., Loktin A. V., 2005, Astronomy Letters, 31, 171
  * Quillen et al. (2011) Quillen A. C., Dougherty J., Bagley M. B., Minchev I., Comparetta J., 2011, MNRAS, 417, 762
  * Rattenbury et al. (2007) Rattenbury N. J., Mao S., Sumi T., Smith M. C., 2007, MNRAS, 378, 1064
  * Reid et al. (2019) Reid M. J., et al., 2019, ApJ, 885, 131
  * Riello et al. (2021) Riello M., et al., 2021, A&A, 649, A3
  * Russeil (2003) Russeil D., 2003, A&A, 397, 133
  * Sellwood & Binney (2002) Sellwood J. A., Binney J. J., 2002, MNRAS, 336, 785
  * Shu (2016) Shu F. H., 2016, ARA&A, 54, 667
  * Taylor & Cordes (1993) Taylor J. H., Cordes J. M., 1993, ApJ, 411, 674
  * Toomre (1981) Toomre A., 1981, in Fall S. M., Lynden-Bell D., eds, Structure and Evolution of Normal Galaxies. pp 111–136
  * Vázquez et al. (2008) Vázquez R. A., May J., Carraro G., Bronfman L., Moitinho A., Baume G., 2008, ApJ, 672, 930
  * Veselova & Nikiforov (2020) Veselova A. V., Nikiforov I., 2020, Research in Astronomy and Astrophysics, 20, 209
  * Wegg et al. (2015) Wegg C., Gerhard O., Portail M., 2015, MNRAS, 450, 4050
  * Xu et al. (2013) Xu Y., et al., 2013, ApJ, 769, 15
  * Xu et al. (2018) Xu Y., Hou L.-G., Wu Y.-W., 2018, Research in Astronomy and Astrophysics, 18, 146
  * Xu et al. (2021) Xu Y., Hou L. G., Bian S. B., Hao C. J., Liu D. J., Li J. J., Li Y. J., 2021, A&A, 645, L8