# Takin: An Open-Source Software for Experiment Planning, Visualisation, and
Data Analysis

Tobias Weber<EMAIL_ADDRESS>Robert Georgii  Peter Böni  Physik-
Department E21, Technische Universität München, James-Franck-Str. 1, 85748
Garching, Germany  Heinz-Maier-Leibnitz-Zentrum (MLZ), Technische Universität
München, Lichtenbergstr. 1, 85747 Garching, Germany

###### Abstract

Due to the instrument’s non-trivial resolution function, measurements on
triple-axis spectrometers require extra care from the experimenter in order to
obtain optimal results and to avoid unwanted spurious artefacts. We present a
free and open-source software system that aims to ease many of the tasks
encountered during the planning phase, in the execution and in data treatment
of experiments performed on neutron triple-axis spectrometers. The software is
currently in use and has been successfully tested at the MLZ, but can be
configured to work with other triple-axis instruments and instrument control
systems.

This was the 2016 pre-print of our software paper [1], which is now available
at http://dx.doi.org/10.1016/j.softx.2016.06.002 and is licensed under CC BY
4.0 (http://creativecommons.org/licenses/by/4.0/). As pre-print, it is a
modified form of the published paper. Furthermore, small updates have been
included in the present version, e.g. concerning the references.

###### keywords:

triple-axis spectroscopy , instrument control , reciprocal and real space
visualisation , instrument resolution , resolution convolution

††journal: SoftwareX

##  1 Motivation and significance

The invention of the neutron triple-axis spectrometer by Brockhouse [2] has
entailed a huge boost in the understanding of dynamics in solid matter. Today,
triple-axis spectroscopy has become a standard tool in the investigation of
materials in such diverse fields as magnetism, lattice dynamics,
superconductivity, and critical phenomena.

A small selection of very recent examples include an investigation of the
interactions governing the nematic order of iron-based superconductors [3],
helimagnetic band formation in the chiral magnet MnSi [4], and the behaviour
of the helimagnetic bands at the conical-ferromagnetic phase transition [5].
But despite their ubiquity in neutron physics, the triple-axis spectrometer is
often poorly understood by users in terms of its very complicated resolution
function. An understanding of the instrumental resolution is crucial for a
correct planning of measurements and data treatment, as resolution effects
often lead to artefacts in the spectra – called “spurions” – which cannot be
easily distinguished from genuine excitations.

In practice, resolution analysis has led to isolated per-instrument software
applications which are often either very complicated to use, too specialised,
or lacking the required features (e.g. [6, 7, 8]). In this article we present
the free and open-source software system Takin, which is capable of performing
all the required steps in one package and which is general enough to be easily
adaptable to new triple-axis instruments. The software furthermore features a
fully-integrated graphical user interface (GUI) and is therefore much more
accessible than previous command-line tools.

##  2 Software description

Figure 1: Takin main view. Given the sample and instrument parameters, the
software displays the reciprocal lattice (left) and the real lattice (right).
In the reciprocal view, the neutron structure factors $S_{hkl}$ of the Bragg
peaks are shown and the first Brillouin zone is calculated, taking into
account any elastic peaks which are forbidden by the sample space group
symmetry. The scattering triangle can be dragged & dropped and is translated
into the corresponding triple-axis instrument position (not shown) and
resolution (see figure 2). The real space view shows the crystal lattice, the
Wigner-Seitz cell and the atomic positions in the selected unit cell.  Figure
2: View of the corresponding resolution function of the triple-axis
spectrometer. The software performs a realtime update and live display of the
instrument resolution. Supported algorithms include the classical Cooper-
Nathans and Popovici methods as well as the novel Eckold-Sobolev approach to
the resolution function.

Takin is a software system encompassing multiple goals. Its basic feature,
namely a quick visualisation of both the crystal lattice in real and in
reciprocal space and the instrument configuration, is similar to the vTAS
software [8]. But whereas vTAS is mainly focused on simulating realistic
instrumental limits, e.g. walls in the instrument space, Takin does not
include such features, but instead focuses on combining crystallographic
information and instrument visualisation with a special emphasis on resolution
calculation. Moreover, a further emphasis of the software is data treatment
and resolution deconvolution using the Monte-Carlo approach.

###  2.1 Software Functionalities

The main window of the program is divided into two views: a reciprocal and a
real-space representation. The reciprocal space includes the Bragg reflections
that are allowed by the space group’s symmetry operations as well as the
elastic structure factors $S_{hkl}=\left|F_{hkl}\right|^{2}$ for each peak.
Given the allowed peaks, the first Brillouin zone is calculated and displayed.
Furthermore, the scattering triangle is drawn in the reciprocal space and can
be either moved by mouse or imported from scan files (the supported formats
include the ones used by MLZ, PSI, ILL and NIST). Alternatively, the software
can receive the current live instrument position by a network connection to
the triple-axis control software. With respect to network connections, two
instrument control systems are currently supported: NICOS [9], the networked
instrument control system which is used at the MLZ, and SICS [10], the SINQ
instrument control software used at the PSI. Moreover, the object-oriented
network architecture in Takin allows for an easy integration of further
triple-axis control systems. The second part of the screen is reserved for the
real-space view, which shows the real crystal lattice and the Wigner-Seitz
cell. The unit cell and its atom positions, which are calculated by the given
space group symmetry operations, are shown in a projected representation in
this view. As with the reciprocal view, the displayed crystal plane can be
freely chosen.

The software warns if the spectrometer – either the real or the simulated one
– is in a spurious scattering configuration. To that end Takin includes a
module to check for the most important false peak conditions as given by
Shirane et al. [11, Ch. 6]. These include the Currat-Axe spurion which – for
example – appears if $\mathbf{q}$ is parallel to $\mathbf{k_{f}}$ and
$\left|\mathbf{k_{f}}+\mathbf{q}\right|=\left|\mathbf{k_{i}}\right|$. In this
configuration the analyser will Bragg-scatter an elastic incoherent signal
from the sample which subsequently appears as a faux-inelastic peak. Other
spurious signals that are checked for include the higher-order reflections
from the monochromator and the analyser and the Bragg tail. The latter appears
in the $\mathbf{q}$ direction transverse to the lattice vector $\mathbf{G}$
and results from the sloping of the resolution ellipse with respect to the
direction perpendicular to the mean $\mathbf{Q}$ and $E$ position.
Furthermore, a module for calculating powder peaks is included for a quick
assessment of spurious powder scattering from the sample environment, for
instance from aluminium or copper. Please note that more specific spurious
signals – such as thermal diffuse scattering by the analyser – are not
considered in the software.

The entire software package and its associated library [12] are written in
modern C++11 [13] and make extensive usage of the Boost [14] C++ template
library. The graphical user interface is based on Qt [15] and the Qwt [16]
plotting library. Tables with physical information, especially the properties
of the 230 space group types and the scattering lengths, are obtained from the
Clipper [17] library and online from the NIST tables [18], respectively. We
use the magnetic form factors from the ILL tables [19].

####  2.1.1 Resolution

##### Resolution Calculation

The basic resolution calculation module features a C++ re-implementation of
the Cooper-Nathans [20] and Popovici [21] algorithms from Rescal 5 [7].
Furthermore, a novel implementation of the Eckold-Sobolev approach to the
resolution function [22] is available. The latter algorithm is especially
interesting for a realistic treatment of the common feature of monochromator
and analyser focusing, which increases the neutron flux, but also increases
the resolution volume. Each algorithm calculates the covariance matrix for the
wavevector $\mathbf{Q}$ and energy transfer $E$ of the neutron at a specific
instrument position. The inverse of the covariance matrix is the instrumental
resolution matrix $R$ – a quadric describing a four-dimensional ellipsoid. The
angles of the ellipsoid with respect to the $\mathbf{Q}$ and $E$ axes and the
widths with respect to the ellipsoid axes are determined using the principal
axis theorem which involves solving the eigenvector problem. The software uses
a specialised, fast implementation for this reduced eigenvector calculation
involving symmetrical matrices. The result is displayed (see figure 2) using
two- and three-dimensional projections of the ellipsoid in either the
$\left<\mathbf{Q_{\parallel}},\mathbf{Q_{\perp}},\mathbf{Q_{z}},E\right>$ or
the fractional $\left<h,k,l\right>$ crystal coordinate system.

The native C++ implementation of the algorithms is fast enough for a real-time
updating and displaying of the resolution function even on low-end hardware.
On multi-core or multi-processor systems further performance is gained by
taking advantage of the symmetric treatment of the monochromator and the
analyser resolution volumes in the Eckold-Sobolev method which allows for an
efficient parallelisation of the calculations.

##### Resolution Convolution

Measurements using a triple-axis spectrometer do not yield the actual
dynamical structure factor $S\left(\mathbf{Q},\omega\right)$ (where
$\omega=E/\hbar$), but instead give a convolution of $S$ with the instrumental
resolution $R$. Using an externally supplied theoretical
$S\left(\mathbf{Q},\omega\right)$ function, an experiment can thus be
simulated, e.g. for planning future measurements, given the $R$ matrix.

Takin contains a module for Monte-Carlo simulation of the convolution integral
(see figure 3). The Monte-Carlo approach includes generating a predefined
number of neutrons inside the four-dimensional resolution ellipsoid, followed
by sampling $S\left(\mathbf{Q},\omega\right)$ with the $\mathbf{Q}$ and $E$
coordinates of the Monte-Carlo neutrons.

The $S\left(\mathbf{Q},\omega\right)$ function for the excitations to be
simulated can be supplied as either a table or as a Python [23] script or – in
case performance is crucial – as a native C++ class which has to be derived
from a special interface. For tabulated $S\left(\mathbf{Q},\omega\right)$
values, the values from the table are used to generate a four-dimensional
binary search tree (a k-d tree [24]) which can then be very efficiently used
to query nearest positions by the Monte-Carlo simulation. The Monte-Carlo
simulation itself is sped up by making use of all available processor cores
for all model sources except the Python scripts (the latter being due to
threading limitations in the Python interpreter).

##### Resolution Convolution Fit

Free parameters can be defined in the different theoretical
$S\left(\mathbf{Q},\omega\right)$ models. In the Python model – for instance –
all global variables are automatically determined upon loading the script and
offered as fit parameters in the convolution fitter tool. User-selected free
parameters are used to perform a $\chi^{2}$ minimisation in order to get the
optimal fit of the resolution convolution to experimental measurements. For
performing the minimisation, the Minuit software library [25] is used.

Due to its complexity, the fitter program module is currently only accessible
via a command-line interface and takes a special job file format as input. In
order to be in unison with the otherwise fully GUI-based design of the
software a fully graphical convolution fitter is planned for the near future.

##  3 Illustrative Example

Figure 3: Upper panel: For experiment planning, a Monte-Carlo integration of
the resolution convolution integral can be quickly done by selecting a
previously defined crystal, instrument and theoretical
$S\left(\mathbf{Q},\omega\right)$ model. A scan is then simulated along the
defined path. Lower panel: Shown are the results of the scan defined in the
upper panel, i.e. an energy scan at $\mathbf{Q}=\left(1.043\ 0.957\ 0\right)$.
This example plot depicts a 4D convolution of the MIRA instrumental function
with the helimagnons in a MnSi single crystal using the theoretical model
developed by J. Waizner and M. Garst and implemented by G. Brandl [4]. Four
helimagnon bands can be resolved in the focusing direction, which for
$k_{i}=const.$ and the given $\mathbf{Q}$ is at $E<0$. In the non-focusing
direction at $E>0$ the bands cannot be resolved and appear as one broad peak.
The incoherent elastic contribution is not visible.

As Takin contains many different triple-axis features, we restrict ourselves
to showing one specific functionality here: planning an experiment using a
simulated scan, which can be fully and quickly done in the GUI. In the
convolution dialog (figure 3, upper panel), a file describing the crystal and
the instrument geometry have to be given along with a model file containing
the dynamical structure factor $S\left(\mathbf{Q},\omega\right)$. The crystal
and instrument files are defined beforehand in the Takin main view (figure 1).

As described earlier, the $S\left(\mathbf{Q},\omega\right)$ model file can be
either a C++ class, a table or a Python script. In the present example, we use
the chiral magnet MnSi as a sample, MIRA [26, 27] as the instrument and the
helimagnon model from [4] as a theoretical $S\left(\mathbf{Q},\omega\right)$
function. An arbitrary scan path in $\left(\mathbf{Q},E\right)$ space can next
be defined and simulated. Here, we simulate an energy scan at
$\mathbf{Q}=\left(1.043\ 0.957\ 0\right)\mathrm{rlu}$ using a fixed incident
wavenumber, $k_{i}=1.4\,\mathrm{\mbox{\AA}}^{-1}$. The resulting plot (figure
3, lower panel) shows that scans at this $k_{i}$ are viable on MIRA as four
individual helimagnon bands can be resolved for $E<0$, which is the focusing
direction for the defined positions. For the defocused measurement positions
at $E>0$ the simulation shows that the bands cannot be resolved anymore and
instead only one broad, smeared-out peak will be discernible.

Figure 4: An example convolution fit to experimental data [5]. The data points
were measured on MIRA [26, 27], the lines are obtained from the Takin
convolution fitter using the model given in [4]. The momentum transfer $q$ of
this data set differs from the one shown in figure 3.

In a further step, Takin can be used for convolution fits to experimentally
measured data. We can confirm an excellent agreement of the calculated
resolution functions using a focusing monochromator and measurements of the
helimagnon bands in MnSi [5] as depicted in figure 4.

##  4 Impact

With Takin we hope to ease the life of triple-axis users in all aspects of
experiment planning, execution and data analysis. In the planning phase it can
be used to quickly and graphically assess if the experiment is possible on the
given instrument: for example if the angles are in the allowed range; whether
the scattering triangle is closed; or the resolution is good enough.

During the experiment, the software shows the live position of the instrument
in reciprocal space, draws the corresponding resolution function at this
position and warns of possible spurions. After the experiment, Takin can be
used “offline” to do convolution simulations and convolution fits to the
measured spectra providing a choice between different resolution algorithms.

The latest version of the software is currently in regular use at the MLZ
instruments MIRA [26, 27] and PANDA [28]. Moreover, the software has been
successfully tested to work at the SINQ instrument TASP [29] in “live mode”
during a magnon measurement.

The first comprehensive applications for the resolution fitting capabilities
of Takin will be presented in two upcoming papers by the authors. The first
one is concerned with the evolution of the helimagnon band structure in the
vicinity of the conical-ferromagnetic phase transition of MnSi [5], the second
one examines an anomalous phonon behaviour at the cubic-tetragonal transition
in the spinel vanadate MgV2O4 [30].

##  5 Conclusion

We presented a software system for performing calculations and visualisations
of inelastic neutron scattering with a special emphasis on the resolution
function of triple-axis spectrometers.

For future developments, it is planned to offer general support for time-of-
flight spectrometers. Currently, the software already contains a module for
the resolution calculation and convolution using the Violini method [31].
Furthermore the convolution fitter will be equipped with advanced scripting
capabilities and also included in the main GUI.

## Acknowledgements

This work was supported by the German Research Foundation (DFG) and the
Technische Universität München within the funding programme Open Access
Publishing. The work at MIRA was supported by the DFG under GE971/5-1.

We wish to thank Max Kugler and Georg Brandl for many discussions about
triple-axis spectroscopy and instrumentation, and for providing the helimagnon
Python script for the convolution example. We also want to thank Petr Čermák,
Astrid Schneidewind, and Bertrand Roessli for their support and for providing
PANDA and EIGER resolution parameters, respectively.

This work is based on experiments performed at the instrument MIRA operated by
FRM II at the Heinz Maier-Leibnitz Zentrum (MLZ), Garching, Germany.

The software’s DOI is: 10.5281/zenodo.3961491.

## Author Contributions

T.W. designed and developed the software and originally wrote the manuscript
as part of his Ph.D. thesis [32]. P.B. and R.G. contributed to the manuscript
and to the design of the software through many discussions over the years.
T.W. conducted the experiment at MIRA, R.G. was the local contact, and P.B.
the supervisor.

## References

  * [1] T. Weber, R. Georgii, P. Böni, Takin: An open-source software for experiment planning, visualisation, and data analysis, SoftwareX 5 (2016) 121–126.  doi:10.1016/j.softx.2016.06.002.
URL http://www.sciencedirect.com/science/article/pii/S2352711016300152

  * [2] B. N. Brockhouse, in: Inelastic Neutron Scattering in Solids and Liquids, Vol. 487, International Atomic Energy Agency, Vienna, 1961.
  * [3] Q. Wang, Y. Shen, B. Pan, Y. Hao, M. Ma, F. Zhou, P. Steffens, K. Schmalzl, T. R. Forrest, M. Abdel-Hafiez, X. Chen, D. A. Chareev, A. N. Vasiliev, P. Bourges, Y. Sidis, H. Cao, J. Zhao, Strong interplay between stripe spin fluctuations, nematicity and superconductivity in FeSe, Nature Materials 15 (2) (2016) 159–163.  doi:10.1038/nmat4492.
URL http://dx.doi.org/10.1038/nmat4492

  * [4] M. Kugler, G. Brandl, J. Waizner, M. Janoschek, R. Georgii, A. Bauer, K. Seemann, A. Rosch, C. Pfleiderer, P. Böni, M. Garst, Band Structure of Helimagnons in MnSi Resolved by Inelastic Neutron Scattering, Phys. Rev. Lett. 115 (2015) 097203.  doi:10.1103/PhysRevLett.115.097203.
URL http://link.aps.org/doi/10.1103/PhysRevLett.115.097203

  * [5] T. Weber, J. Waizner, G. S. Tucker, R. Georgii, M. Kugler, A. Bauer, C. Pfleiderer, M. Garst, P. Böni, Field dependence of nonreciprocal magnons in chiral MnSi, Phys. Rev. B 97 (2018) 224403.  doi:10.1103/PhysRevB.97.224403.
URL https://link.aps.org/doi/10.1103/PhysRevB.97.224403

  * [6] J. Šaroun, J. Kulda, RESTRAX – a program for TAS resolution calculation and scan profile simulation, Physica B: Condensed Matter 234-236 (1997) 1102–1104, proceedings of the First European Conference on Neutron Scattering.  doi:10.1016/S0921-4526(97)00037-9.
URL http://www.sciencedirect.com/science/article/pii/S0921452697000379

  * [7] D. A. Tennant, D. F. McMorrow, Rescal for Matlab: a computational package for calculating neutron TAS resolution functions, software (1995).
URL https://www.ill.eu/instruments-support/computing-for-science/cs-
software/all-software/matlab-ill/rescal-for-matlab/

  * [8] M. Boehm, A. Filhol, Y. Raoul, J. Kulda, W. Schmidt, K. Schmalzl, E. Farhi, The vTAS suite: A simulator for classical and multiplexed three-axis neutron spectrometers, Nuclear Instruments and Methods in Physics Research Section A: Accelerators, Spectrometers, Detectors and Associated Equipment 697 (2013) 40 – 44.  doi:10.1016/j.nima.2012.08.050.
URL http://www.sciencedirect.com/science/article/pii/S0168900212009217

  * [9] G. Brandl, E. Faulhaber, C. Felder, J. Krüger, A. Lenz, B. Pedersen, NICOS 2 - Networked Instrument Control System (2016).
URL https://forge.frm2.tum.de/cgit/cgit.cgi/frm2/nicos/nicos-core.git

  * [10] H. Heer, M. Könnecke, D. Maden, The SINQ instrument control software system, Physica B: Condensed Matter 241-243 (1997) 124–126, proceedings of the International Conference on Neutron Scattering.  doi:10.1016/S0921-4526(97)00528-0.
URL http://www.sciencedirect.com/science/article/pii/S0921452697005280

  * [11] G. Shirane, S. M. Shapiro, J. M. Tranquada, Neutron Scattering with a Triple-Axis Spectrometer: Basic Techniques, Cambridge University Press, 2002, ISBN 978-0521411264.
  * [12] T. Weber, tlibs - a physical-mathematical C++ template library (2015).
URL https://github.com/t-weber/tlibs

  * [13] B. Stroustrup, The C++ Programming Language, 4th edition, Addison-Wesley, 2013, ISBN 978-0321563842.
  * [14] B. Dawes, D. Abrahams, R. Rivera, Boost C++ Library, software (2016).
URL http://www.boost.org

  * [15] The Qt Company, Qt C++ Library, software (2016).
URL https://www.qt.io

  * [16] U. Rathmann, J. Wilgen, Qwt Plotting Library, software (2016).
URL http://qwt.sourceforge.net

  * [17] K. D. Cowtan, Clipper library (2013).
URL http://www.ysbl.york.ac.uk/~cowtan/clipper

  * [18] A. Munter, Neutron scattering length tables (online) (1999).
URL https://www.ncnr.nist.gov/resources/n-lengths/list.html

  * [19] P. J. Brown, Magnetic form factor tables (online) (1998).
URL https://www.ill.eu/sites/ccsl/ffacts/

  * [20] M. J. Cooper, R. Nathans, The resolution function in neutron diffractometry. I. The resolution function of a neutron diffractometer and its application to phonon measurements, Acta Crystallographica 23 (3) (1967) 357–367.  doi:10.1107/S0365110X67002816.
URL http://dx.doi.org/10.1107/S0365110X67002816

  * [21] M. Popovici, On the resolution of slow-neutron spectrometers. IV. The triple-axis spectrometer resolution function, spatial effects included, Acta Crystallographica Section A 31 (4) (1975) 507–513.  doi:10.1107/S0567739475001088.
URL http://dx.doi.org/10.1107/S0567739475001088

  * [22] G. Eckold, O. Sobolev, Analytical approach to the 4d-resolution function of three axes neutron spectrometers with focussing monochromators and analysers, Nuclear Instruments and Methods in Physics Research Section A: Accelerators, Spectrometers, Detectors and Associated Equipment 752 (2014) 54 – 64.  doi:10.1016/j.nima.2014.03.019.
URL http://www.sciencedirect.com/science/article/pii/S0168900214002940

  * [23] G. Van Rossum, F. L. Drake, The Python Language Reference Manual, Network Theory Ltd., 2011, ISBN 978-1906966140.
  * [24] J. L. Bentley, Multidimensional binary search trees used for associative searching, Commun. ACM 18 (9) (1975) 509–517.  doi:10.1145/361002.361007.
URL http://doi.acm.org/10.1145/361002.361007

  * [25] I. Antcheva, M. Ballintijn, B. Bellenot, M. Biskup, R. Brun, N. Buncic, P. Canal, D. Casadei, O. Couet, V. Fine, L. Franco, G. Ganis, A. Gheata, D. G. Maline, M. Goto, J. Iwaszkiewicz, A. Kreshuk, D. M. Segura, R. Maunder, L. Moneta, A. Naumann, E. Offermann, V. Onuchin, S. Panacek, F. Rademakers, P. Russo, M. Tadel, ROOT – A C++ framework for petabyte data storage, statistical analysis and visualization, Computer Physics Communications 182 (6) (2011) 1384 – 1385.  doi:10.1016/j.cpc.2011.02.008.
URL http://www.sciencedirect.com/science/article/pii/S0010465511000701

  * [26] R. Georgii, K. Seemann, Mira: Dual wavelength band instrument, Journal of large-scale research facilities 1 (2015) A3.  doi:10.17815/jlsrf-1-21.
URL http://nbn-resolving.de/urn:nbn:de:0001-jlsrf-1-21-0

  * [27] R. Georgii, T. Weber, G. Brandl, M. Skoulatos, M. Janoschek, S. Mühlbauer, C. Pfleiderer, P. Böni, The multi-purpose three-axis spectrometer (TAS) MIRA at FRM II, Nuclear Instruments and Methods in Physics Research Section A: Accelerators, Spectrometers, Detectors and Associated Equipmentdoi:10.1016/j.nima.2017.09.063.
URL http://www.sciencedirect.com/science/article/pii/S0168900217310331

  * [28] A. Schneidewind, P. Čermák, Panda: Cold three axes spectrometer, Journal of large-scale research facilities 1 (2015) A12.  doi:10.17815/jlsrf-1-35.
URL http://nbn-resolving.de/urn:nbn:de:0001-jlsrf-1-35-5

  * [29] F. Semadeni, B. Roessli, P. Böni, Three-axis spectroscopy with remanent benders, Physica B: Condensed Matter 297 (1-4) (2001) 152–154, proceeding of the Third International Workshop on Polarised Neutrons.  doi:10.1016/S0921-4526(00)00860-7.
URL http://www.sciencedirect.com/science/article/pii/S0921452600008607

  * [30] T. Weber, B. Roessli, C. Stock, T. Keller, K. Schmalzl, F. Bourdarot, R. Georgii, R. A. Ewings, R. S. Perry, P. Böni, Transverse acoustic phonon anomalies at intermediate wave vectors in ${\mathrm{mgv}}_{2}{\mathrm{o}}_{4}$, Phys. Rev. B 96 (2017) 184301.  doi:10.1103/PhysRevB.96.184301.
URL https://link.aps.org/doi/10.1103/PhysRevB.96.184301

  * [31] N. Violini, J. Voigt, S. Pasini, T. Brückel, A method to compute the covariance matrix of wavevector-energy transfer for neutron time-of-flight spectrometers, Nuclear Instruments and Methods in Physics Research Section A: Accelerators, Spectrometers, Detectors and Associated Equipment 736 (2014) 31 – 39.  doi:10.1016/j.nima.2013.10.042.
URL http://www.sciencedirect.com/science/article/pii/S016890021301423X

  * [32] T. Weber, Ph.D. thesis, Physikdepartment E21, Technische Universität München, Garching, Germany (2016).

## Required Metadata

## Current code version

Nr. |  Code metadata description |
---|---|---
C1 |  Current code version |  1.0
C2 |  Permanent link to code/repository used for this code version |  https://forge.frm2.tum.de/cgit/cgit.cgi/frm2/mira/tastools.git/ and https://github.com/t-weber/takin
C3 |  Legal Code License |  GPL Version 2
C4 |  Code versioning system used |  git
C5 |  Software code languages, tools, and services used |  C++ 11, CMake 3
C6 |  Compilation requirements, operating environments |  Linux, OS X, Unix like, Windows (via MinGW); GCC 4.8 (or later) or Clang, Boost, Qt 4 or 5, Qwt 5 or 6, Clipper, Minuit 2
C7 |  If available Link to developer documentation/manual |  https://forge.frm2.tum.de/cgit/cgit.cgi/frm2/mira/tastools.git/plain/doc/index_help.html
C8 |  Support email for questions |<EMAIL_ADDRESS>
Table 1: Code metadata


# Reliability and limitations of inferring birth radii in the Milky Way disk

Yuxi(Lucy) Lu,1,2, Tobias Buck3, Ivan Minchev3, Melissa K. Ness1,4
1Department of Astronomy, Columbia University, 550 West 120th Street, New
York, NY, USA
2American Museum of Natural History, Central Park West, Manhattan, NY, USA
3Leibniz Institute for Astrophysics Potsdam, An der Sternwarte 16, 14482
Potsdam, Germany
4Center for Computational Astrophysics, Flatiron Institute, 162 5th Avenue,
Manhattan, NY, USA
E-mail<EMAIL_ADDRESS>
(Accepted XXX. Received YYY; in original form ZZZ)

###### Abstract

Recovering the birth radii of observed stars in the Milky Way is one of the
ultimate goals of Galactic Archaeology. One method to infer the birth radius
and the evolution of the ISM metallicity assumes a linear relation between the
ISM metallicity with radius at any given look-back time. Here we test the
reliability of this assumption by using 4 zoom-in cosmological hydrodynamic
simulations from the NIHAO-UHD project. We find that one can infer precise
birth radii only when the stellar disk starts to form, which for our modeled
galaxies happens $\sim$ 10 Gyr ago, in agreement with recent estimates for the
Milky Way. At later times the linear correlation between the ISM metallicity
and radius increases, as stellar motions become more ordered and the azimuthal
variations of the ISM metallicity start to drop. The formation of a central
bar and perturbations from mergers can increase this uncertainty in the inner
and outer disk, respectively.

###### keywords:

Galaxy: evolution – Galaxy: abundances

††pubyear: 2022††pagerange: Reliability and limitations of inferring birth
radii in the Milky Way disk–Reliability and limitations of inferring birth
radii in the Milky Way disk

##  1 Introduction

It is widely accepted that stars move away from their birth locations over
time in the process commonly referred to as radial migration. Radial migration
was first suggested by Grenon (1972, 1989), in which they found old metal-rich
stars in the solar neighborhood that exhibited kinematic and abundance
properties consistent with those of the inner disk. Stars can move away from
their birth radii either because they have high eccentricities or due to a
permanent change in their angular momentum. The latter process is also
referred to as radial migration or mixing. Stars oscillating around their
guiding radii can blur the intrinsic relations, however, this process can be
accounted for by estimating the stellar guiding radius. On the other hand,
migrating stars are largely indistinguishable in their kinematics from stars
born locally (e.g. Sellwood & Binney, 2002; Roškar et al., 2008; Minchev &
Famaey, 2010). Moreover, they can flatten the intrinsic abundance gradient in
the disk significantly over time (e.g. see Figure 5 from Minchev et al.,
2013).

As a result, to truly understand the enrichment history of the Milky Way (MW),
we need to connect the present-day observables to the birth properties of
stars. Without such knowledge, reconstructing MW’s formation history remains
unreachable, and one can easily infer an incorrect enrichment history of the
Galaxy. One example, pointed out in Minchev et al. (2018), is if we attempt to
infer the temporal evolution of the ISM metallicity gradient with stars in
mono-age bins, assuming age is equivalent to look-back time, we would end up
with a flatter ISM metallicity gradient evolution compared to the truth.

The most ambitious approach for recovering common birth environments of stars
is perhaps that of chemical tagging (Freeman & Bland-Hawthorn, 2002). This
method assumes stars that were born in the same cluster are chemically
homogeneous, and that clusters can differentiated from one another given
sufficiently precise and numerous individual abundance measurements from
different nucloesynthetic channels. Then, by grouping stars chemically, one
can trace the stars back to their original birth aggregates in the Galaxy.
However, this method is a challenge for the current data (e.g. Bland-Hawthorn
et al., 2010; Ness et al., 2018; DeSilva et al., 2007; Manea et al., 2022;
Nelson et al., 2021; Ting & Weinberg, 2022; Ratcliffe et al., 2020; Ness et
al., 2022). Nevertheless, even with the currently available abundance
measurements, we can still hope to obtain a subset of birth properties of
stars with sensible assumptions, one of which is the birth radius of a star.

Minchev et al. (2018) presented a largely model-independent approach for
estimating both stellar birth radii, Rbirth, and the Galactic ISM metallicity
evolution with time and radius, [Fe/H](t, R), using metallicity and age
estimates from the local HARPS sample (Adibekyan et al., 2012). The method
relied on the following assumptions: (1) the gas is well mixed azimuthally,
(2) stars are born from a narrow metallicity range at a certain radius at any
given time, (3) the ISM [Fe/H](t, R) evolved smoothly with both radius and
time, and (4) the MW disk formed inside-out. Feltzing et al. (2020) concluded
the method shows promising results by inferring Rbirth with APOGEE DR14 data
and comparing the resulted Rbirth distributions to those in a hydrodynamic
simulation from Agertz et al. (2021). Other derivations of the radial
migration strength can be found in the works by, e.g., Frankel et al. (2019,
2020); Lian et al. (2022).

Although the method described in Minchev et al. (2018) has shown to be
successful, the information about birth radius in the MW is lost, and
therefore, the inferred Rbirth from data cannot be truly validated. Moreover,
inferring a star’s birth radius with its metallicity and age requires a
correlation between the ISM [Fe/H] and R at the time when the star was born
(look-back time = age of the star). In this letter, we will show that precise
Rbirth can only be inferred when a rotationally supported stellar disk is
starting to form. We use MW-like cosmological zoom-in simulations from NIHAO-
UHD (Buck, 2020; Buck et al., 2020), in which we have access to the enrichment
history and Rbirth of stars. We will do this by linking the evolution of ISM
metallicity, the orderliness of the gas-particle orbits, mergers, and the
existence of a stellar bar with the intrinsic uncertainty on inferring Rbirth.

In section 3.1, we investigate at what time does the correlation between the
ISM metallicity and R is strong enough to infer reliable Rbirth. Once the
correlation is established, the dispersion around the metallicity – R linear
relation can still lead to intrinsic uncertainty in inferring Rbirth. This
scatter can be created from the azimuthal variation of ISM metallicity caused
by satellite infall, the spiral arms, the bar, or possibly other perturbations
(e.g. Bellardini et al., 2021; Spitoni et al., 2019; Sánchez-Menguiano et al.,
2020). We examine the connection between the ISM metallicity gradient,
azimuthal ISM metallicity variation, the orderliness of the gas-particle
orbits, and the uncertainty of inferring Rbirth, assuming a linear relation
between the ISM metallicity and R at all times in section 3.2. Finally, we
discuss the implication for inferring Rbirth in the MW in section 4.

##  2 Data & Methods

###  2.1 Simulated Data

The simulated galaxies used in this work were first presented by (Buck et al.,
2018). They are taken from a suite of high-resolution cosmological
hydrodynamical simulations of MW-mass galaxies from the NIHAO-UHD project
(Buck et al., 2020). A modified version of the smoothed particle hydrodynamics
(SPH) solver GASOLINE2 (Wadsley et al., 2017) was used to calculate the
simulations from cosmological initial conditions. Star formation and feedback
(both energetic and chemical elements) are modeled following the prescriptions
in (Stinson et al., 2006; Stinson et al., 2013). The simulations used in this
work trace oxygen and iron separately (see Buck et al., 2021, for a recent
udate to the chemical enrichment) and adopted a metal diffusion algorithm
between particles as described in Wadsley et al. (2008). The parameters for
simulations used in this letter are taken from Buck et al. (2020). The numbers
in the naming represent the total Dark Matter (DM) halo mass in unites of
solar mass. The time separation between every two consecutive output snapshots
is $\sim$ 200 Myr, and all of the galaxies have been integrated until redshift
0. For more information on the simulation details and galaxy properties, we
refer the readers to Buck et al. (2020).

##  3 Results

To infer the birth radius for a star from its [Fe/H] and stellar age alone, a
correlation between the ISM metallicity and R at a look-back time equal to the
stellar age of the star is required. However, as shown in simulations (e.g.
Gurvich et al., 2022; Bird et al., 2013; Bellardini et al., 2022),
observations of high-redshift galaxies (e.g. Ciucă et al., 2021; Wisnioski et
al., 2015; Übler et al., 2019; Elmegreen et al., 2007) and observations of the
MW (e.g. Belokurov & Kravtsov, 2022; Conroy et al., 2022), hot, turbulent gas
is likely to be common in early galaxy formation. This means a correlation
between the ISM metallicity and R at early times may not exist. Similarly, as
pointed out in Minchev et al. (2018), if the MW started with a flat gradient,
their method would fail for stars $\gtrsim$ 10 Gyr.

###  3.1 When can we infer reliable Rbirth from stellar [Fe/H] and age?

Figure 1 summarizes the formation histories of the 4 simulated galaxies of
interest. The metallicity evolution of the stars is shown in the background,
colored by Rbirth. The blue and red histograms show the instantaneous mass
increase ($\Delta$M/M) within 2Rhalf of the galaxies for star and gas
particles, respectively, where Rhalf is the half-light radius of the galaxies
at each look-back time, calculated by (13.8 Gyr - cosmic time). $\Delta$M/M is
calculated by $\Delta$M/M = (Mn+1-Mn)/Mn, where $n$ is the snapshot number and
$n$ = 0 represents the snapshot at a look-back time of 13.8 Gyr. The increase
in stellar mass is a combination of in-situ star formation and ex-situ star
accretion. Similarly, the gas mass increase results from in-situ supernova and
ex-situ gas accretion. The red dashed line shows the time when the stellar
disk is starting to form, which is defined here as most stars are on circular
orbits (see Figure 10 of Buck et al., 2020, or Figure 3 left plot). The black
dashed line in the bottom figure shows the formation of the bar taken from
Buck et al. (2018) Figure 10. The vertical grey/orange lines are the times
when there is a satellite > 1%/10% total mass ratio in the main galaxy’s viral
radius, indicating on-going merger events. The red dashed lines, black dashed
lines, and the grey/orange lines are used throughout the figures in this
letter. The green line shows the absolute value of the Pearson correlation
coefficient (|PCC|) between the metallicity and R of the gas particles at each
look-back time when the snapshot is available. PCC is a statistical quantity
that identifies the linear correlation between two variables, in which -1 and
1 represent a perfectly negative or positive linear correlation, respectively.
Two strongly correlated variables typically have a |PCC| $>$ 0.7, meaning
$\sim$ 50% of the variation has been accounted for in the data by fitting a
line. We have indicated the line of |PCC| = 0.7 in the Figures (cyan dashed
line).

It is clear that all 4 galaxies formed bottom-up, meaning they originated from
a series of galaxy mergers and gas accretion events, as the increase in star
and gas mass are $>$ 100% early on (see blue/red histograms). Most stars are
born in the inner galaxy early on, while later on the birth radius
distribution broadens, indicating inside-out formation. During this early
time, when $\Delta$M/M is large, the |PCC| between the ISM metallicity and R
are mostly $<$ 0.7, meaning no strong linear correlation exists, and thus, we
cannot infer precise Rbirth from a star’s metallicity and age alone. The
correlation between the ISM metallicity and R remains low until a well-ordered
stellar disk begins to form and when a significant amount of stars are formed
in the outer disk ($>\sim$ 5 kpc; indicated by the red dashed line) for most
galaxies. However, once the frequent mergers/accretions have subsided and a
well ordered stellar disk begins to form, the |PCC| stays above 0.7 for all
galaxies111Due to the frequent merger history at low redshift, galaxy g6.96e11
has a lower —PCC— compared to the other galaxies after its stellar disk has
begun to form. . This implies that inferring birth radii with [Fe/H]and age is
most precise and accurate after the turbulent phase of galaxy formation has
ended, which for the MW has been suspected to be $\sim$ 10 Gyr ago based on
comparing simulations with observations (Belokurov & Kravtsov, 2022)222Conroy
et al. (2022) also pointed out the thick disk formed $\sim$ 13 Gyr ago defined
by star formation efficiency. However, for the purpose of this study, we find
the definition of disk formation from Belokurov & Kravtsov (2022) fits best
with this study..

Figure 1: Formation history of 4 galaxies used in this work. The metallicity
evolution of the stars is shown in the background, colored by Rbirth. The blue
and red histograms show the instantaneous mass increase ($\Delta$M/M) within
2*Rhalf of the galaxies for star and gas, respectively, where Rhalf is the
half-light radius of the galaxy at each look-back time. These values are
calculated using (Mn+1-Mn)/Mn, where $n$ is the snapshot number and $n$ = 0
represents the snapshot at a look-back time of 13.8 Gyr. The green lines show
the absolute Pearson correlation coefficient (|PCC|) between [Fe/H] and R for
the ISM at each look-back time when the snapshot is available. The cyan dashed
lines indicate PCC = 0.7, which is when the ISM [Fe/H] and R are strongly
correlated. The red dashed lines show the approximate times when the stellar
disks started to from. The black dashed line in the bottom figure shows the
formation of the bar taken from Buck et al. (2018) Figure 10. The vertical
grey/orange lines are the times when there is a satellite > 1%/10% total mass
ratio in the main galaxy’s viral radius, indicating on-going merger events.
This figure implies that inferring birth radii with [Fe/H]and age is most
accurate after the turbulent phase of galaxy formation has ended, and when the
stellar disk has begun to form, which for the MW is suspected to be $\sim$ 10
Gyr ago based on combining observations and simulations (Belokurov & Kravtsov,
2022).

###  3.2 ISM metallicity gradient and azimuthal variation, circularisation of
the orbits, and the uncertainty on Rbirth

How well does a straight line fit the metallicity gradient at all times, and
what uncertainty will that bring in inferring Rbirth? The answer lies within
Figure 2. The top plot shows the metallicity gradient evolution over time, and
the solid lines display the smoothed version of the points, using a 1D
Gaussian filter after the disk has started to form for each galaxy. The slopes
are calculated by fitting a line between the gas particles’ metallicity and R
within 25 kpc of the galaxy. The red dashed line shows a natural log function
to guide the eye. Minchev et al. (2018) derived a natural-log functional form
for the metallicity gradient evolution, and it fits well with most of the
galaxies soon after their disks have begun to form. As pointed out earlier,
the g6.96e11 galaxy has frequent mergers/infalls at low redshift, causing a
more chaotic ISM evolution. The metallicity dispersion ([Fe/H]disp) around the
linear relation between the ISM metallicity and R at each snapshot is shown in
the middle plot. Similar to the conclusion from section 3.1, the ISM is more
turbulent at earlier times and thus, the |PCC| between ISM [Fe/H] and R is
low, and [Fe/H]disp is high. However, as the stellar disk started to form,
[Fe/H]disp settled to $<$ 0.1 dex.

We can also estimate the intrinsic uncertainty in inferring Rbirth,
$\sigma_{Rb}$, assuming 1) we know the linear fits of the ISM metallicity at a
given look-back time, and that 2) there is no error on the observables ([Fe/H]
and age; see the bottom plot of Figure 2). The intrinsic uncertainty is
estimated using $\sigma_{Rb}$ = [Fe/H]disp/(|$d$[Fe/H]ISM/$d$R|). It is clear
that even with perfectly measured observables, the intrinsic scattering around
the ISM metallicity evolution, assuming a linear relation between the ISM
[Fe/H] and R at all times, yields an uncertainty on the Rbirth measurement to
be $\gtrsim$ 2 kpc for all 4 simulated galaxies. For the MW, based on the
estimation of the current metallicity gradient of $\sim$ -0.05 dex/kpc with a
dispersion of 0.1 dex derived from Cepheids (Luck, 2018), the intrinsic
uncertainty would be $\sigma_{\rm Rbirth}$$\sim$ 2 kpc, which is comparable
with the simulated galaxies. However, as pointed out by (Luck, 2018), the
intrinsic dispersion of 0.1 dex is smaller than the measurement uncertainty of
0.14 dex, implying the actual dispersion can be smaller. We also note that the
metallicity gradient in the Galaxy is still not well constrained (e.g. Spina
et al., 2022; Eilers et al., 2022)

Figure 2: Metallicity gradient evolution (top), metallicity dispersion around
the linear fit, [Fe/H]disp (middle), and uncertainty in inferring
Rbirth($\sigma_{Rb}$) assuming the metallicity gradient evolution is known and
that there is no error on [Fe/H] and age measurements (bottom). The
uncertainty, $\sigma_{Rb}$, is calculated by [Fe/H]disp/(|$d$[Fe/H]ISM/$d$R|).
The solid lines are the smoothed version of the points, using a 1D Gaussian
filter, after the stellar disk begins to form. The red line in the top plot
shows a natural log function to guide the eye.

We can also connect the ISM evolution and $\sigma_{\rm Rbirth}$ with the
“orderliness” of the stars orbit (a proxy for circularity of the star’s
orbits), $<v_{\phi}>/\sigma_{v\phi}$, for the ISM, where $<v_{\phi}>$ is the
average azimuthal velocity, and $\sigma_{v\phi}$ is the azimuthal velocity
dispersion. Figure 3 left column plots show the time evolution of
$<v_{\phi}>/\sigma_{v\phi}$ for the 4 simulated galaxies. We look at
$<v_{\phi}>/\sigma_{v\phi}$ in 6 different radial bins between 0-18 kpc. The
solid lines are the smoothed version of the points, using a 1D Gaussian filter
after the stellar disk begins to form. The black dashed line in the bottom
figure shows the formation of the bar taken from Buck et al. (2018) Figure 10.
Darker vertical lines mean more satellites exist in the viral radius at that
time. Figure 3 right plot shows the same as the left but with the ISM
metallicty variation in a small radial bin ($\sigma_{\rm[Fe/H]}$) at 6
different radius between 0-18 kpc. The metallicity variation increases towards
the outer disk, and can also be affected by the bar, agreeing with findings
from other simulations and chemical evolution models (e.g. Bellardini et al.,
2021; Spitoni et al., 2019).

Figure 3 suggest when the stellar disk begins to form, the orbits transition
rapidly into a more ordered configuration and the ISM metallicity azimuthal
variation also started to decrease. The rapid increase in
$<v_{\phi}>/\sigma_{v\phi}$ and the decrease of ISM metallicity azimuthal
variation shows a fast transition from the turbulent gas phase to the quiesent
disk phase (see also Gurvich et al., 2022). This rapid phase transition also
leads to a well-established metallicity gradient with a high correlation
(|PCC| $>$ 0.7). This means, after the rapid increase in the orderliness of
the gas-particle orbits, a star formed from the ISM has a well-constraint
metallicity at a certain radius, meaning we can infer birth radii with smaller
$\sigma_{\rm Rbirth}$ (Figure 2). However, late mergers can induce turbulence
in gas, and therefore, increase $\sigma_{\rm[Fe/H]}$ and $\sigma_{\rm Rbirth}$
in mostly the outer galaxy, c.f. dispersion in $<v_{\phi}>/\sigma_{v\phi}$.
The bar can also induce variation in metallicity in the inner galaxy. As a
result, the uncertainty in inferred Rbirth for the inner Galaxy after the
formation of the bar, and the outer Galaxy during a time of merger could be
higher.

Figure 3: Time evolution of $<v_{\phi}>/\sigma_{v\phi}$ (left) and the
azimuthal ISM metallicity variation, $\sigma_{\rm[Fe/H]}$, (right) for the
four simulated galaxies, in which $<v_{\phi}>$ is the average azimuthal
velocity, and $\sigma_{v\phi}$ is the azimuthal velocity dispersion.
$<v_{\phi}>/\sigma_{v\phi}$ can be taken as the “orderliness” of the gas-
particle orbits, or a proxy for the circularity of the star’s orbits. We look
at $<v_{\phi}>/\sigma_{v\phi}$ in 6 different radial bins and
$\sigma_{\rm[Fe/H]}$ in 6 narrow radial bins ranging from 0-18 kpc. The solid
lines are the smoothed version of the points, using a 1D Gaussian filter after
the stellar disk begins to form.

##  4 Conclusions and Discussion of the implications for the MW

By using MW-like simulations from the NIHAO-UHD project, we are able to
connect the phases of galaxy formation with the possibility of inferring birth
radii of stars. In particular, this letter shows:

  * •

Precise birth radii can be inferred from stellar ages, guiding center radii,
and metallicities only once the stellar disk starts to form - $\sim$ 10 Gyr
ago for the MW (Belokurov & Kravtsov, 2022).

  * •

At later times the linear correlation between the ISM metallicity and R
increases, as stellar motions become more ordered and the azimuthal variations
of the ISM metallicity decrease by an order of magnitude (see Figure 1 and
Figure 3).

  * •

A central bar and perturbations from mergers can increase the intrinsic
uncertainty of inferred Rbirth in the inner and outer disk by up to 50%,
respectively (see Figure 3).

  * •

Considering the observed present-day ISM metallicity gradient, derived from
Cepheids, of -0.05 dex/kpc with a scatter of 0.1 dex (Luck, 2018), the minimum
intrinsic uncertainty on inferring Rbirth in the MW is $\sim$ 2 kpc. However,
as pointed out in Luck (2018), this scatter is smaller than the measurement
uncertainty (0.14 dex), and therefor, the true scatter and the intrinsic
uncertainty is likely to be smaller than 2 kpc.

However, in this study, we only investigated a linear relation between the ISM
metallicity and R. Future work should investigate the possibility that other
functional forms exist, that better describe this correlation. We also only
used the ISM gradient derived from Cepheids and other studies have found
slightly different gradients ranging from around -0.05 to -0.07 dex/kpc (e.g.
Anders et al., 2017; Maciel & Andrievsky, 2019). Studies derive local
metallicity gradient using neutral gas and HII regions have found a slightly
shallower gradient of $\sim$ -0.04 dex/kpc (e.g. De Cia et al., 2021; Esteban
et al., 2022). Most of the studies either do not report the metallicity
intrinsic dispersion around the relations or have a scatter that is on the
order of the measurement uncertainty. As a result, the intrinsic uncertainty
on Rbirth for the MW is not constrained. We also want to point out that the
evolution of the simulated galaxies may not match exactly with that of the MW.
One discrepancy is that the disk seems to form earlier for the MW compared to
simulations (see also Belokurov & Kravtsov, 2022), which could indicate
missing/incorrect physics in the simulations. Future work is needed to further
investigate this discrepancy.

## Acknowledgements

Melissa Ness is supported in part by a Sloan Fellowship. Tobias Buck
acknowledges support from the European Research Council under ERC-CoG grant
CRAGSMAN-646955. This research made use of pynbody Pontzen et al. (2013). We
gratefully acknowledge the Gauss Centre for Supercomputing e.V. (www.gauss-
centre.eu) for funding this project by providing computing time on the GCS
Supercomputer SuperMUC at Leibniz Supercomputing Centre (www.lrz.de). This
research was carried out on the High Performance Computing resources at New
York University Abu Dhabi. We greatly appreciate the contributions of all
these computing allocations.

## Data Availability

The data underlying this article will be shared on reasonable request to the
corresponding author.

## References

  * Adibekyan et al. (2012) Adibekyan V. Z., Sousa S. G., Santos N. C., Delgado Mena E., González Hernández J. I., Israelian G., Mayor M., Khachatryan G., 2012, A&A, 545, A32
  * Agertz et al. (2021) Agertz O., et al., 2021, MNRAS, 503, 5826
  * Anders et al. (2017) Anders F., et al., 2017, A&A, 600, A70
  * Bellardini et al. (2021) Bellardini M. A., Wetzel A., Loebman S. R., Faucher-Giguère C.-A., Ma X., Feldmann R., 2021, MNRAS, 505, 4586
  * Bellardini et al. (2022) Bellardini M. A., Wetzel A., Loebman S. R., Bailin J., 2022, arXiv e-prints, p. arXiv:2203.03653
  * Belokurov & Kravtsov (2022) Belokurov V., Kravtsov A., 2022, arXiv e-prints, p. arXiv:2203.04980
  * Bird et al. (2013) Bird J. C., Kazantzidis S., Weinberg D. H., Guedes J., Callegari S., Mayer L., Madau P., 2013, ApJ, 773, 43
  * Bland-Hawthorn et al. (2010) Bland-Hawthorn J., Krumholz M. R., Freeman K., 2010, ApJ, 713, 166
  * Buck (2020) Buck T., 2020, MNRAS, 491, 5435
  * Buck et al. (2018) Buck T., Ness M. K., Macciò A. V., Obreja A., Dutton A. A., 2018, ApJ, 861, 88
  * Buck et al. (2020) Buck T., Obreja A., Macciò A. V., Minchev I., Dutton A. A., Ostriker J. P., 2020, MNRAS, 491, 3461
  * Buck et al. (2021) Buck T., Rybizki J., Buder S., Obreja A., Macciò A. V., Pfrommer C., Steinmetz M., Ness M., 2021, MNRAS, 508, 3365
  * Ciucă et al. (2021) Ciucă I., Kawata D., Miglio A., Davies G. R., Grand R. J. J., 2021, MNRAS, 503, 2814
  * Conroy et al. (2022) Conroy C., et al., 2022, arXiv e-prints, p. arXiv:2204.02989
  * De Cia et al. (2021) De Cia A., Jenkins E. B., Fox A. J., Ledoux C., Ramburuth-Hurt T., Konstantopoulou C., Petitjean P., Krogager J.-K., 2021, Nature, 597, 206
  * DeSilva et al. (2007) DeSilva G. M., Freeman K. C., Asplund M., Bland-Hawthorn J., Bessell M. S., Collet R., 2007, AJ, 133, 1161
  * Eilers et al. (2022) Eilers A.-C., Hogg D. W., Rix H.-W., Ness M. K., Price-Whelan A. M., Mészáros S., Nitschelm C., 2022, ApJ, 928, 23
  * Elmegreen et al. (2007) Elmegreen D. M., Elmegreen B. G., Ravindranath S., Coe D. A., 2007, ApJ, 658, 763
  * Esteban et al. (2022) Esteban C., E. M. J., Garćıa-Rojas J., Arellano-Córdov K. Z., 2022, arXiv e-prints, p. arXiv:2204.12258
  * Feltzing et al. (2020) Feltzing S., Bowers J. B., Agertz O., 2020, MNRAS, 493, 1419
  * Frankel et al. (2019) Frankel N., Sanders J., Rix H.-W., Ting Y.-S., Ness M., 2019, ApJ, 884, 99
  * Frankel et al. (2020) Frankel N., Sanders J., Ting Y.-S., Rix H.-W., 2020, ApJ, 896, 15
  * Freeman & Bland-Hawthorn (2002) Freeman K., Bland-Hawthorn J., 2002, ARA&A, 40, 487
  * Grenon (1972) Grenon M., 1972, in Cayrel de Strobel G., Delplace A. M., eds, IAU Colloq. 17: Age des Etoiles. p. 55
  * Grenon (1989) Grenon M., 1989, Ap&SS, 156, 29
  * Gurvich et al. (2022) Gurvich A. B., et al., 2022, arXiv e-prints, p. arXiv:2203.04321
  * Lian et al. (2022) Lian J., et al., 2022, MNRAS, 511, 5639
  * Luck (2018) Luck R. E., 2018, AJ, 156, 171
  * Maciel & Andrievsky (2019) Maciel W. J., Andrievsky S., 2019, arXiv e-prints, p. arXiv:1906.01686
  * Manea et al. (2022) Manea C., Hawkins K., Maas Z. G., 2022, arXiv e-prints, p. arXiv:2201.07809
  * Minchev & Famaey (2010) Minchev I., Famaey B., 2010, ApJ, 722, 112
  * Minchev et al. (2013) Minchev I., Chiappini C., Martig M., 2013, A&A, 558, A9
  * Minchev et al. (2018) Minchev I., et al., 2018, MNRAS, 481, 1645
  * Nelson et al. (2021) Nelson T., Ting Y.-S., Hawkins K., Ji A., Kamdar H., El-Badry K., 2021, ApJ, 921, 118
  * Ness et al. (2018) Ness M., et al., 2018, ApJ, 853, 198
  * Ness et al. (2022) Ness M. K., Wheeler A. J., McKinnon K., Horta D., Casey A. R., Cunningham E. C., Price-Whelan A. M., 2022, ApJ, 926, 144
  * Pontzen et al. (2013) Pontzen A., Roškar R., Stinson G. S., Woods R., Reed D. M., Coles J., Quinn T. R., 2013, pynbody: Astrophysics Simulation Analysis for Python
  * Ratcliffe et al. (2020) Ratcliffe B. L., Ness M. K., Johnston K. V., Sen B., 2020, ApJ, 900, 165
  * Roškar et al. (2008) Roškar R., Debattista V. P., Stinson G. S., Quinn T. R., Kaufmann T., Wadsley J., 2008, ApJ, 675, L65
  * Sánchez-Menguiano et al. (2020) Sánchez-Menguiano L., Sánchez S. F., Pérez I., Ruiz-Lara T., Galbany L., Anderson J. P., Kuncarayakti H., 2020, MNRAS, 492, 4149
  * Sellwood & Binney (2002) Sellwood J. A., Binney J. J., 2002, MNRAS, 336, 785
  * Spina et al. (2022) Spina L., Magrini L., Cunha K., 2022, Universe, 8, 87
  * Spitoni et al. (2019) Spitoni E., Cescutti G., Minchev I., Matteucci F., Silva Aguirre V., Martig M., Bono G., Chiappini C., 2019, A&A, 628, A38
  * Stinson et al. (2006) Stinson G., Seth A., Katz N., Wadsley J., Governato F., Quinn T., 2006, MNRAS, 373, 1074
  * Stinson et al. (2013) Stinson G. S., et al., 2013, MNRAS, 436, 625
  * Ting & Weinberg (2022) Ting Y.-S., Weinberg D. H., 2022, ApJ, 927, 209
  * Übler et al. (2019) Übler H., et al., 2019, ApJ, 880, 48
  * Wadsley et al. (2008) Wadsley J. W., Veeravalli G., Couchman H. M. P., 2008, MNRAS, 387, 427
  * Wadsley et al. (2017) Wadsley J. W., Keller B. W., Quinn T. R., 2017, MNRAS, 471, 2357
  * Wisnioski et al. (2015) Wisnioski E., et al., 2015, ApJ, 799, 209