# Unifying Simulation and Inference with Normalizing Flows

Haoxing Du<EMAIL_ADDRESS>Department of Physics, University of
California, Berkeley, CA 94720, USA  Claudius Krause
<EMAIL_ADDRESS>Institut für Theoretische Physik, Universität
Heidelberg, Philosophenweg 12, 69120 Heidelberg, Germany  Institute of High
Energy Physics (HEPHY), Austrian Academy of Sciences (OeAW), Georg-Coch-Platz
2, A-1010 Vienna, Austria  Vinicius Mikuni<EMAIL_ADDRESS>National Energy
Research Scientific Computing Center, Berkeley Lab, Berkeley, CA 94720, USA
Benjamin Nachman<EMAIL_ADDRESS>Physics Division, Lawrence Berkeley
National Laboratory, Berkeley, CA 94720, USA  Berkeley Institute for Data
Science, University of California, Berkeley, CA 94720, USA  Ian Pang
<EMAIL_ADDRESS>NHETC, Department of Physics and Astronomy,
Rutgers University, Piscataway, NJ 08854, USA  David Shih
<EMAIL_ADDRESS>NHETC, Department of Physics and Astronomy, Rutgers
University, Piscataway, NJ 08854, USA

###### Abstract

There have been many applications of deep neural networks to detector
calibrations and a growing number of studies that propose deep generative
models as automated fast detector simulators. We show that these two tasks can
be unified by using maximum likelihood estimation (MLE) from conditional
generative models for energy regression. Unlike direct regression techniques,
the MLE approach is prior-independent and non-Gaussian resolutions can be
determined from the shape of the likelihood near the maximum. Using an ATLAS-
like calorimeter simulation, we demonstrate this concept in the context of
calorimeter energy calibration.

††preprint: HEPHY-ML-24-01

##  I Introduction

Detector calibrations are one of the most important and foundational tasks in
experimental physics. In particle and nuclear physics, the largest calibration
step is usually a simulation-based111The residual correction to account for
differences between data and simulation is usually smaller. Methods in this
paper may also be applicable to data-based corrections in certain cases.
correction to ensure that the reported properties of the measured particles
are unbiased. Detector simulations are complex and accurate, but can only run
forward in time and the resolutions are non-trivial. Combined, these
properties mean that we cannot simply invert the simulation to predict true
quantities given measured ones.

This is particularly acute for highly-segmented detectors where many
individual channels are activated for a single particle as the measured phase
space can be high dimensional and complex. A common example of this setting is
calorimeter reconstruction. Unlike tracking detectors that aim to minimally
disrupt the trajectory of a particle, calorimeters are designed to stop
particles and the resulting showers inside dense materials are challenging to
parse and are highly stochastic. Traditional calibration methods are optimized
using relatively low-dimensional summary statistics. Such calibrations have
enabled many science results, but the ultimate precision cannot be reached
until we utilize all of the available low-level information for event
reconstruction.

Machine learning provides a set of tools that can analyze hadronic final
states holistically to achieve the best precision. Recent examples supporting
this claim include the energy calibrations of single hadrons [1, 2, 3, 4, 5,
6, 7, 8, 9, 10, 11], jets [12, 13, 14, 15, 16, 17, 18, 19, 20, 21], and global
event properties [22, 23, 24, 25] at colliders. These techniques can
automatically make use of finer segmentation to improve the resolution of
reconstructed particle energies [4, 10]. While machine learning has broader
utility than reconstructing individual hadron showers within a calorimeter, we
focus on this case because it is particularly challenging and a critical
component of event reconstruction at the Large Hadron Collider and elsewhere.

While existing machine learning methods are promising, most approaches have at
least two undesirable features: they are only point estimates and are prior-
dependent. Firstly, most approaches produce a single estimate without
quantifying ‘uncertainties’. In particle/ nuclear physics, the spread of the
difference between the inferred estimates and the true value is usually called
the resolution, even though it is part of Uncertainty Quantification (UQ) in
the broader machine learning literature. Secondly, techniques that employ
standard machine learning regression tools are not universal: their
performance depends on the distribution of examples (the ‘prior’, e.g. uniform
in energy) used during training. In fact, this prior dependence often results
in a large calibration bias.

Over the last years, some partial solutions to these two central challenges
have been proposed. Generalized numerical inversion [26, 27, 28] is a prior-
independent calibration approach, although it does not scale well to many
output dimensions. The Gaussian Ansatz [17, 18] is a local maximum likelihood
estimator, and thus prior-independent (more in Sec. II). For UQ, loss
functions can be modified to estimate quantiles in addition to just the mean,
median or mode [13, 14] and deep generative models can estimate the exact
resolution function [24, 25]. The Gaussian Ansatz produces a local estimate of
the resolution and so is currently the only method that is both prior-
independent and that goes beyond a point estimate. However, the Gaussian
Ansatz only gives a Gaussian approximation of the resolution. We note that it
is also possible to use Bayesian versions [29, 30, 31, 32, 33] of neural
networks to go beyond point estimates by providing the uncertainty of the
output. However, the uncertainty would still be prior-dependent.

In this paper, we propose a prior-independent method for detector calibrations
that produces per-shower resolution estimates. The approach is based on deep
generative models with access to the explicit likelihood. In particular, we
use normalizing flows222It is also possible to extract the likelihood in many
dimensions from a diffusion model [34] or if the inputs are discretized, from
an autoregressive model [35]. [36, 37, 38, 39, 40, 41]. These machine learning
models are invertible functions with a tractable Jacobian so that one can be
used for both sampling and density estimation. The core idea behind a
normalizing flow is that one starts with a simple random variable with a known
probability density (e.g. a Gaussian) and then applies a series of
transformations with the tractable Jacobian. This results in a complex
probability density that can also be computed via the change of variables
formula. The probability density can be made conditional by letting each
transformation depend on the conditional quantity.

Normalizing flows are state-of-the-art as calorimeter simulation surrogate
models [42, 43, 44, 45, 46, 47, 48, 49]. Additionally, normalizing flows have
shown similar good performance in other tasks in high-energy physics [50, 51,
52, 53, 54, 55, 56, 33, 57, 58, 59, 60, 61, 62, 24, 63, 64, 65, 66, 24, 25,
67, 68, 69, 70, 71, 72]. In this paper, we unify simulation and inference by
showing how models trained for generation can be reused for calibration. As a
demonstration of this approach, we use CaloFlow [42, 43], a normalizing flow-
based calorimeter surrogate model. This machine learning model is trained on
single-pion showers from an extended version [73] of the CaloGAN dataset [74,
75, 76] that now includes both a sampling calorimeter setup [77] and a
hadronic component. The calibration task is to predict the incident pion
energy given the distribution of energies recorded in the cells of the
calorimeter. Importantly, we demonstrate the following key advantages of our
approach:

  1. 1.

Zero-shot calibration:
Normalizing flow-based models trained for surrogate modelling can be reused
for calibration by repurposing the probability density as a likelihood without
any additional model retrainings.

  2. 2.

Access to per-shower resolution:
With access to the complete likelihood, we can compute the exact resolution
function. This allows us to obtain per-shower resolution estimates which give
us additional information about how close the predicted energy is to the true
energy. Furthermore, it enables us to estimate asymmetries in the resolution
function (not possible with the GA).

  3. 3.

Less biased calibration:
We achieve a smaller calibration bias compared to a typical direct regression
approach for the calorimeter setup considered in this work.

This paper is organized as follows. Machine learning-based calibration methods
are introduced in Sec. II. Here we also provide precise definitions for bias
and resolution used in this paper. In Sec. III, we summarize the CaloFlow
algorithm and discuss the results of the calorimeter example which demonstrate
the key advantages listed above. The paper ends with conclusions and outlook
in Sec. IV. We collect some details on the used network architectures in an
appendix.

##  II Methods

Given samples from a forward model (usually a physics-based simulation) $X\sim
p_{X|Z}(x|z)$ for measured values $X$ and true values $Z$, the goal of
calibration is to estimate $Z$ from $X$. In our notation, capital letters
represent random variables and lower case letters represent realizations of
those random variables. Both the measured and true values can be
multidimensional, although we will focus on the common case where $Z$ is one
dimensional and $X\in\mathbb{R}^{N}$.

The usual assumption when deriving a simulation-based calibration is that the
detector response is correct and universal, i.e.
$p_{X|Z}^{\text{train}}(x|z)=p_{X|Z}^{\text{test}}(x|z)$ for different physics
processes (i.e. different $p(z)$) or for actual data.

There are two important quantities related to calibration performance — the
bias and the resolution.

  1. 1.

The bias of a calibration is the deviation between the central tendency of the
inferred estimate $\hat{z}$ and the true reference value $z$. Any measure of
central tendency can be used to measure closure, such as the median or mode.
In this paper, we will focus on the mode333Usually bias refers to the mean,
but experimentalists typically calibrate to the mode. of the inferred
estimates given a fixed $z$ value. To compute the mode in practice, for a
given $z$, we look at $K$ measured values $x_{1},x_{2}\dots x_{K}$
corresponding to that $z$. Each measured value $x_{i}$ undergoes calibration
to derive an inferred estimate $\hat{z}_{i}$. The mode is subsequently
determined based on the collection of inferred estimates
$\\{\hat{z}_{1},\hat{z}_{2},\dots,\hat{z}_{K}\\}$.

  2. 2.

The resolution is the spread of the difference between the inferred estimates
and the true value. In this work, we have two relevant definitions:

     * •

Full resolution: For fixed $z$, it is defined as half of the 68% confidence
interval about the mode of a collection of inferred estimates.

     * •

Per-event resolution: For each $x$, it is defined as half of the 68%
confidence interval of the corresponding maximum likelihood estimate (see Sec.
II.2).

The full resolution is not to be confused with the per-event resolution. The
full resolution is defined for a collection of inferred estimates, whereas the
per-event resolution is defined for each inferred estimate.

###  II.1 Direct Regression

Most proposals for deep learning-based calibration directly regress $Z$ from
$X$ using a loss function like the mean-squared error (MSE):

$\displaystyle L[f]=\sum_{i}(f_{\rm MSE}(x_{i})-z_{i})^{2}\,,$ (1)

for a neural network $f_{\rm MSE}$. Using the calculus of variations, one can
show that with enough training data and sufficiently flexible neural network
architecture and training protocol, the solution to Eq. 1 is the average value
of $Z$ given $X=x$:

$\displaystyle f_{\rm MSE}(x)$ $\displaystyle=\langle Z|X=x\rangle$ (2)
$\displaystyle=\int dz\,z\,p_{Z|X}^{\text{train}}(z|x)$ (3)
$\displaystyle=\int
dz\,z\,p_{X|Z}^{\text{train}}(x|z)\,\frac{p_{Z}^{\text{train}}(z)}{p_{X}^{\text{train}}(x)}\,.$
(4)

Note that $f_{\rm MSE}$ only gives a point estimate of $z$ for each $x$ and
does not have access to per-event resolutions. In other words, there is no
calibration uncertainty quantified by $f_{\rm MSE}$.

For a given $z$, the bias can be computed as the deviation between $z$ and the
mode of a collection of $f_{\rm MSE}(x_{i})$, and the full resolution at $z$
is then computed based on the 68% confidence interval about the mode. Note
that the bias is defined as a function of $z$.444Alternatively, defining the
bias as a function of the inferred values would result in the bias being
dependent on the prior used during testing (i.e. $p^{\rm test}_{Z}(z)$). This
bias also cannot be corrected since we do not know a priori what $p^{\rm
test}_{Z}(z)$ is used during calibration. Hence, we cannot correct for this
bias as we do not have access to $z$ while performing the calibration.

The challenge with Eq. 4 is that it depends on $p_{Z}^{\text{train}}$ even
when we assume that the detector response is universal. Unsurprisingly then,
the bias of direct regression also depends on the training dataset. In Ref.
[18], the prior dependence of MSE-based calibration is shown to be the main
source of large calibration bias.

###  II.2 Maximum Likelihood Inference

Maximum likelihood estimation (MLE) involves finding the $Z$ that maximizes
the likelihood of the data,

$\hat{Z}=\underset{z}{\text{argmax}}\,p_{X|Z}(X|z)\,.$ (5)

Since we assume that $p_{X|Z}$ is universal and Eq. 5 only depends on this
likelihood, then $\hat{Z}$ is universal. We can go beyond the point estimate
and use the likelihood function $p_{X|Z}$ to obtain the per-event resolution
for each maximum likelihood estimate $\hat{z}$. The access to per-event
resolutions is a major advantage of MLE over direct regression in calibration.
Examples showcasing this advantage are discussed in Sec. III.4.

Similar to the direct regression case, the bias and full resolution at a fixed
$z$ value can be defined based on the mode and 68% confidence interval about
the mode.

Although the MLE-based calibration is always prior-independent, the prior
dependence does not guarantee that the method is always unbiased. Only in
certain scenarios (e.g. 1D Gaussian noise model such that $p_{X|Z}(x|z)\sim
N(z,\sigma)$) can the MLE calibration be shown to be unbiased. Certain
detector responses $p_{X|Z}$ may result in biased calibration. In general,
prior independence is a necessary but insufficient condition for unbiased
calibration. Even so, we show in Sec. III.3 that MLE-based calibration results
in a smaller bias compared to direct regression for an ATLAS-like calorimeter
setup. This suggests that in our application, the dependence on training prior
is a bigger contributing factor to the bias than the detector response.

Also, the challenge with MLE is that we usually do not know $p_{X|Z}$
explicitly. Nevertheless, we are able to sample from this conditional density
by running a simulation. The Gaussian Ansatz [17, 18] approaches this problem
by estimating the likelihood locally around the mode with a clever neural
network parameterization and learning rate schedule. In this work, we instead
aim to learn the entire likelihood using a neural network — an approach
sometimes referred to as neural likelihood estimation in the simulation-based
inference literature [78, 79, 80, 81]. This is substantially more work, but we
observe that the work may already be done: given a fast simulation based on
neural networks with access to the likelihood, we can use it for calibration
in addition to generation without requiring any additional retraining.

Our tool of choice is the normalizing flow (NF). NFs are neural networks that
are optimized using maximum likelihood estimation. With access to an estimate
of the full likelihood, we can study both Gaussian and non-Gaussian aspects of
the resolution.

##  III Regressing particle incident energy

In this section, we compare the performance of direct regression and NF in
regressing the particle incident energy from the resulting calorimeter shower
information.

###  III.1 Dataset

The events are $\pi^{+}$ calorimeter showers from a new version [73] of the
CaloGAN dataset which we generated with Geant4 [82, 83, 84] for this study. We
generated 100k showers with $\pi^{+}$ incident energy uniformly distributed in
the range [1,100] GeV. Additionally, we generated 100k showers with incident
energy log-uniformly distributed in the same range. The calorimeter setup used
in the original dataset only included the electromagnetic calorimeter (ECAL).
In this new dataset, we also include the hadronic calorimeter (HCAL) which is
positioned behind the ECAL. Also, the original dataset included energy
contributions from both active and inactive calorimeter layers, whereas this
new sampling calorimeter dataset only includes energy contributions from the
active layers as would be available in practice.

The ECAL is a three-layer sampling calorimeter cube with 480 mm side-length
that is inspired by the ATLAS liquid argon (LAr) electromagnetic calorimeter
[85]. For the ECAL, the active material is LAr and the absorber material is
lead (Pb). Note that the ECAL used here is identical to the one used in Ref.
[77]. However, photon (instead of $\pi^{+}$) showers were considered in that
study. The HCAL is a three-layer sampling calorimeter cube with 2000 mm side-
length located behind the ECAL. For the HCAL, the active material is LAr and
the absorber material is tungsten (W). The sampling fractions of the ECAL and
HCAL are $\sim$ 20% and $\sim$ 1.3%, respectively. The calorimeter showers are
represented as three-dimensional images that are binned in position space.

In this representation, the calorimeter shower geometry is made up of voxels
(volumetric pixels) and the details of the calorimeter voxel dimensions are
included in Table 1. The energy distribution in the calorimeter for uniformly
distributed $\pi^{+}$ incident energies is shown in Fig. 1. Note that the
actual energy deposited in the sampling calorimeter is usually much smaller
than the incident energy. Hence, it is common practice that the energy
deposited is rescaled555In this work, the rescaling is done by dividing the
actual energy deposited in the ECAL or HCAL by their respective sampling
fractions. based on the sampling fraction to ensure that the dataset contains
deposited energy $E_{\rm dep}$ that is close to the incident energy $E_{\rm
inc}$. For this reason, it is possible for some showers to have $E_{\rm
dep}/E_{\rm inc}>1$ as shown in Fig. 1.

Figure 1: Generated versus reference energy distributions in our new ECAL+HCAL calorimeter setup for uniformly distributed incident energies $E_{\rm inc}$. The energy deposited in layer $i$ is denoted by $E_{i}$ and the energy deposited in the calorimeter is denoted by $E_{\rm dep}=\sum_{i=0}^{6}E_{i}$. The ratio of energy deposited in the calorimeter and the incident energy of the particle is denoted by $E_{\text{dep}}/E_{\text{inc}}$. |  | Layer
---
index
|  $z$ length
---
(mm)
|  $\eta$ length
---
(mm)
|  $\phi$ length
---
(mm)
| Number
---
of voxels
ECAL | 0 | 90 | 5 | 160 | $3\times 96$
1 | 347 | 40 | 40 | $12\times 12$
2 | 43 | 80 | 40 | $12\times 6$
HCAL | 3 | 375 | 20.83 | 666.67 | $3\times 96$
4 | 667 | 166.67 | 166.67 | $12\times 12$
5 | 958 | 333.33 | 166.67 | $12\times 6$
Table 1: Dimensions of a calorimeter voxel. The positive $z-$axis (radial
direction in full detector) is the direction of particle propagation, the
$\eta$ direction is along the proton beam axis, and $\phi$ is perpendicular to
$z$ and $\eta$. For the number of voxels, the first (second) number is the
number of bins in the $\phi$ ($\eta$) direction (e.g., $12\times 6$ refers to
12 $\phi$ bins and 6 $\eta$ bins).

###  III.2 Calorimeter Shower simulation with Normalizing Flows

CaloFlow [42, 43] is an approach to fast calorimeter simulation based on
conditional normalizing flows. In the context of fast calorimeter simulation,
CaloFlow is designed to generate the voxel level shower energies
$\vec{\mathcal{I}}$ conditioned on the corresponding incident energies of the
showers $E_{\text{inc}}$ denoted by $p(\vec{\mathcal{I}}|E_{\text{inc}})$. In
particular, it uses a two-flow approach to perform this task: Flow-I is
formulated to learn the probability density $p_{1}(\vec{E}|E_{\text{inc}})$ of
calorimeter layer energies666The layer energy of a given calorimeter layer is
the sum of all the voxel energies in that layer. $\vec{E}$ conditioned on the
incident energy $E_{\text{inc}}$, while Flow-II is designed to learn the
probability density of the normalized voxel level shower energies conditioned
on the calorimeter layer energies
$p_{2}\left(\vec{\mathcal{I}}|\vec{E},E_{\text{inc}}\right)$.

For the regression task discussed in this section, we found that having layer
energy information was sufficient and adding voxel level information did not
significantly increase the regression performance. This is consistent with
Ref. [10], which found that increasing longitudinal segmentation is more
effective than increasing transverse segmentation (equivalent to adding voxel
level information) in improving the energy resolution. As we solely used
$\vec{E}$ as input, only Flow-I of CaloFlow was required to perform the
calibration task.

In our present dataset, $\vec{E}=(E_{0},E_{1},E_{2},E_{3},E_{4},E_{5})$ is
six-dimensional, twice the dimensionality found in previous versions of the
dataset. Hence, there were some minor modifications made to the original
CaloFlow primarily to deal with different dimensionality of the new dataset.
Nevertheless, the main CaloFlow algorithm remains the same. The details of the
architecture and training are outlined in Appendix A. It is known that using
an ensemble of generative models generally helps to improve generation
performance at the distribution level (e.g. Ref. [86]) and would ideally be
done for the task of fast calorimeter surrogate modelling. Hence, we opted to
train two ensembles of 10 independent Flow-I models each. In particular, 10
models are trained on uniformly distributed $E_{\rm inc}$, and the other 10
models are trained on log-uniformly distributed $E_{\rm inc}$.

When generating samples with an ensemble of models, each model in the ensemble
is used to generate an equal fraction of the total number of samples. We show
in Fig. 1 that an ensemble of Flow-I models is able to generate energy
distributions that agree reasonably well with the reference distributions with
some slight mismatch in tails and at places with sharp changes in the
reference distributions.

Once CaloFlow is trained, it can be repurposed for the calibration of particle
incident energy without any additional retraining. For each shower, the
predicted particle incident energy corresponds to the incident energy which
maximizes the likelihood function estimated by CaloFlow. The “full” likelihood
function of the ensemble of Flow-I models is taken to be mean of the
likelihood functions from each of the models in the ensemble. For the
remainder of this paper, this “full” likelihood is what we use to perform the
MLE-based calibration.

###  III.3 Prior-independent and less biased calibration

Using the dataset described in Sec. III.1 to train our models, we predict the
incident energy of the incoming pion $E_{\text{inc}}$ given the energy
deposited in the six calorimeter layers $\vec{E}$. To test the prior
dependence of these methods, we also compared the performance of these methods
when trained on showers with uniformly distributed $E_{\text{inc}}$ versus
log-uniformly distributed $E_{\text{inc}}$ in the range [1,100] GeV.

Our chosen direct regression model is a dense neural network (DNN) with two
fully connected hidden layers, each comprising 256 nodes with ReLU activation
functions. The direct regression DNN is implemented in PyTorch and optimized
with Adam [87] using a batch size of 200 and 200 epochs. We found that using
an ensemble of DNNs helped to reduce the uncertainty of the calibration. Like
for CaloFlow , we opted to have two ensembles, each with 10 independent DNNs.
In particular, 10 models are trained on uniformly distributed $E_{\rm inc}$,
and the other 10 models are trained on log-uniformly distributed $E_{\rm
inc}$. For each shower, the predicted incident energy of the ensemble is then
taken as the mean predicted incident energy of all the individual DNNs in the
ensemble. Each DNN has $\sim$134k model parameters, while each NF has
$\sim$107k model parameters.

In Fig. 2, we show examples of the likelihood evaluated by CaloFlow for three
given calorimeter showers that originate from a 30 GeV pion, 60 GeV pion and
90 GeV pion, respectively. Specifically, the difference between $-2\log
p(\vec{E}|E_{\rm inc})$ and $-2\log p_{\rm max}\equiv-2\log p(\vec{E}|E_{\rm
inc})|_{E_{\rm inc}=E_{\rm pred}}$ is plotted for a range of $E_{\text{inc}}$,
where $E_{\rm pred}$ is the MLE prediction from CaloFlow. Note that the value
of $E_{\rm pred}$ corresponds to the minimum of the blue curve in each of the
three plots in Fig. 2.

To evaluate the bias, the mode of the $E_{\rm pred}/E_{\rm true}$ distribution
has to be computed for each fixed $E_{\rm true}$. In this study, the mode is
estimated using kernel density estimation (KDE)777We used the
scipy.stats.gaussian_kde [88] function to perform the KDE. and the uncertainty
of the estimate is determined by bootstrapping [89]. To obtain the mode of the
$E_{\rm pred}/E_{\rm true}$ distribution at each fixed $E_{\rm true}$, the
steps we took are as follows:

  1. 1.

Draw with replacement $N$ samples from $N$ values of $E_{\rm pred}$, where $N$
is the number of showers in the evaluation dataset for a given fixed $E_{\rm
true}$.

  2. 2.

Perform kernel density estimation of the drawn samples with kernel bandwith
determined using Scott’s rule [90].

  3. 3.

Identify the position of the mode of the estimated density.

  4. 4.

Repeat steps 1-3 for a total of 20 times.

  5. 5.

Compute the mean and standard deviation of the 20 estimated values of the
mode.

Figure 2: Plot of likelihood evaluated using CaloFlow for a given calorimeter
shower originating from a 30 GeV pion (left), 60 GeV pion (middle) and 90 GeV
pion (right), respectively. In each plot, the true incident energy $E_{\rm
true}$ is shown by the vertical black line, and the evaluated likelihood by
the NF is shown in blue. The lower (upper) bounds of the 68% confidence
interval about the MLE is shown by the red vertical dashed lines.

The mean of the multiple mode estimates is taken to be the final mode estimate
at each fixed $E_{\rm true}$. The corresponding standard deviation is taken to
be the uncertainty of the final mode estimate. The mode estimation is
performed for both the uniform and log-uniform cases. In Fig. 3, we show the
the mode of $E_{\rm pred}/E_{\rm true}$ distribution at $E_{\rm
true}\in\\{10,20,30,40,50,60,70,80,90\\}$ GeV for the different regression
methods. Here we used nine evaluation datasets, each with 100k $\pi^{+}$
showers at one of the nine fixed $E_{\rm true}$ values.

Figure 3: Plot of the mode of $E_{\rm pred}/E_{\rm true}$ distribution,
denoted by “Mode$(E_{\rm pred}/E_{\rm true})$”, and the corresponding
uncertainty at different fixed values of $E_{\rm
true}\in\\{10,20,30,40,50,60,70,80,90\\}$ GeV. Results for models trained on
(log-)uniformly distributed $E_{\text{inc}}$ are shown in solid (dashed)
lines.

We observe that the flow-based calibration is generally less biased (closer to
unity) compared to the MSE-based calibration. The larger bias of the MSE-based
calibration is clearly seen at the edges of $E_{\rm true}$ range. Furthermore,
we see that the MSE-based calibration corresponding to the log-uniform prior
is generally more biased than the calibration corresponding the uniform prior
at intermediate energies.

The prior independence of the flow-based calibration is confirmed by observing
consistent biases of CaloFlow trained on uniform and log-uniform priors. For
the MSE-based calibration, there is prior dependence at many of the $E_{\rm
true}$ values used in the evaluation as seen from the difference in the biases
corresponding to the uniform and log-uniform priors. For CaloFlow , the origin
of the dip at $E_{\rm true}=60$ GeV is unclear. However, we observe that the
biases of CaloFlow trained on uniform and log-uniform priors agree with each
other.

###  III.4 Resolution estimation

MLE-based calibration methods, such as the NF and Gaussian Ansatz, are able to
provide estimates of the per-shower resolution and this allows these methods
to go beyond a point estimate by providing per-event resolutions (i.e. per-
shower resolutions in calorimeter context). This gives MLE-based calibration
an important advantage over direct regression.

Using CaloFlow, we perform per-shower resolution estimation by computing the
two $E_{\rm inc}$ values corresponding to the lower (upper) bound of 68%
confidence interval of the MLE prediction. Equivalently, these are the two
$E_{\rm inc}$ values for which $-2\log p+2\log p_{\rm max}=1$ as shown by the
vertical dashed lines in Fig. 2. The per-shower resolution prediction is then
the average difference of the two $E_{\rm inc}$ values from $E_{\rm pred}$.
The full resolution is computed as half of the 68% confidence interval about
the mode of the $E_{\rm pred}/E_{\rm true}$ distribution. In Fig. 4, we
compare the mean predicted per-shower resolution and the full resolution at
each of the nine fixed $E_{\rm true}$ from the NF. The results in this section
were produced using the NF trained on showers with uniformly distributed
incident energies. We observe that the predicted resolution from the NF
closely matches the full resolution. Also, we found that the uncertainty in
the determination of the full resolution to be small and almost unnoticeable
in Fig. 4.

Although not shown in Fig. 4, we found that the MSE-based calibration obtains
a full resolution that closely matches that from the NF. However, it is unable
to provide per-shower resolutions.

Figure 4: Plot of full resolution (solid) and mean predicted per-shower
resolution (dashed) at $E_{\rm true}\in\\{10,20,30,40,50,60,70,80,90\\}$ GeV.

Next we show that the per-shower resolution is a reliable guide to the
calibration accuracy. In Fig. 5, we plot 2D histograms of
$\absolutevalue{E_{\rm pred}/E_{\rm true}-1}$ versus the per-shower resolution
at three different $E_{\rm true}$ values. We observe that showers with small
per-shower resolution tend to have small $\absolutevalue{E_{\rm pred}/E_{\rm
true}-1}$, while those with large per-shower resolutions appear to be evenly
distributed across different values of $\absolutevalue{E_{\rm pred}/E_{\rm
true}-1}$. In other words, a smaller per-shower resolution is indicative of a
more accurately calibrated shower. For $E_{\rm true}=20$ GeV, there is a drop
in the per-shower resolution for $\absolutevalue{E_{\rm pred}/E_{\rm
true}-1}\gtrsim 0.75$ where the NF is relatively confident that the shower
originates from an $E_{\rm inc}>E_{\rm true}$. Nevertheless, we note that the
per-shower resolution in this range of $\absolutevalue{E_{\rm pred}/E_{\rm
true}-1}$ is still larger than that at low values of $\absolutevalue{E_{\rm
pred}/E_{\rm true}-1}$.

Figure 5: 2D histograms of $\absolutevalue{E_{\rm pred}/E_{\rm true}-1}$
versus per-shower resolution at $E_{\rm true}=20$ GeV (left), $E_{\rm
true}=50$ GeV (middle) and $E_{\rm true}=80$ GeV (right).

Since the NF is able to access the complete likelihood, this enables it to
capture non-Gaussian components of the resolution. This is one advantage it
has over the the Gaussian Ansatz approach [17, 18] — the only other method
that is both prior-independent and can provide a per-shower resolution
estimate — which is, by construction, unable to predict any asymmetry. We
quantify this by defining an asymmetry observable based on the upper and lower
bounds of the 68% confidence interval of the MLE prediction. In particular,
the asymmetry is computed as the ratio $\frac{\Delta_{L}}{\Delta_{U}}$, where
$\Delta_{L}$ ($\Delta_{U}$) is the difference between the lower (upper) bound
of the 68% confidence interval and the MLE.

An example is included in Fig. 6 which shows the per-shower resolution
asymmetry predicted by the NF for 100k showers with $E_{\rm true}=50$ GeV. We
observe that $\frac{\Delta_{L}}{\Delta_{U}}$ is largely greater than unity.
This implies that for most showers, lower $E_{\rm inc}$ values tend to have
larger likelihoods. While one might suspect that this resolution asymmetry is
merely an artifact of the flow, we observe in Fig. 7 that the distribution of
$E_{\rm inc}$ at fixed vertical slices of $E_{\rm dep}$ is also asymmetric
towards lower values of $E_{\rm inc}$. The particular shape of the 2D
histogram in Fig. 7 is due to the nature of the $\pi^{+}$ showers. At fixed
$E_{\rm inc}$, fully electromagnetic events will have larger $E_{\rm dep}$,
since all the energy is deposited via electromagnetic interactions. On the
other hand, fully hadronic events will have significantly lower $E_{\rm dep}$
due to energy carried away invisibly by particles like neutrons and neutrinos
produced in nuclear interactions. The $\pi^{+}$ showers considered in this
work contain both electromagnetic and hadronic components, so on a shower-by-
shower basis, the $E_{\rm dep}$ can range anywhere between these two extremes
of a fully electromagnetic shower and a fully hadronic shower [91]. This
results in a tail at high energies in the $E_{\rm dep}$ distribution for fixed
$E_{\rm inc}$, presumably because the typical $\pi^{+}$ shower is mostly
hadronic in nature. Thinking of $p(E_{\rm dep}|E_{\rm inc})$ as a proxy for
the likelihood suggests that the resolution asymmetry is not merely an
artifact of the flow but a genuine characteristic to be expected. We
hypothesize that the ability to predict the exact resolution would allow the
NF to outperform the Gaussian Ansatz in resolution estimation in calibration
tasks where the asymmetries in the resolution are important.

Figure 6: Histogram of asymmetry $\frac{\Delta_{L}}{\Delta_{U}}$ predicted by
CaloFlow for 100k showers with $E_{\rm true}=50$ GeV. Figure 7: 2D histogram
of the deposited energy $E_{\rm dep}$ versus $E_{\rm inc}$ for 100k showers
with uniformly distributed $E_{\rm inc}$. A solid line is drawn across the
mode of $E_{\rm inc}$ at each $E_{\rm dep}$ to guide the eye.

##  IV Conclusions

Using CaloFlow as an example, we show that normalizing flow-based deep
generative models can be repurposed for prior-independent detector
calibration.

We compared the calibration performance of CaloFlow and a DNN direct
regression model on $\pi^{+}$ calorimeter showers from our new sampling
calorimeter dataset [73], and found that the direct regression method is
clearly more biased relative to the flow-based calibration method. This
highlights the advantage of utilizing prior-independent MLE-based calibration
for tasks such as particle energy regression.

We demonstrate a second advantage of using a MLE-based calibration methods
over direct regression by estimating per-shower resolutions. The average
estimated per-shower resolution obtained from CaloFlow closely aligns with the
full resolution. In contrast, the direct regression approach yields only point
estimates and lacks the capacity to make per-shower resolution predictions.
For CaloFlow, we found that smaller per-shower resolutions tend to coincide
with more accurately calibrated showers. This is evidence that the per-shower
resolutions from CaloFlow are reliable and have the potential to be leveraged
for improved calibration. Notably, utilizing a NF for calibration grants
access to the complete resolution function, a feature inaccessible in previous
methods like the Gaussian Ansatz. This characteristic enables the NF to
capture asymmetries in the resolution, thereby expected to confer an advantage
over the Gaussian Ansatz in specific calibration tasks. In future work, it
would be interesting to study how the per-shower resolution information
obtained by the NF can be utilized to further improve the calibration.

## Data and code availability

The datasets used in this study can be found at Ref. [73] and the software to
generate these datasets are located at https://github.com/hep-
lbdl/CaloGAN/tree/two_layer. The machine learning software is at
https://github.com/Ian-Pang/regression_with_CF.

## Acknowledgements

We would like to thank Rikab Gambhir for the helpful discussions related to
the Gaussian Ansatz model. CK would like to thank the Baden-Württemberg-
Stiftung for financing through the program Internationale Spitzenforschung,
project Uncertainties – Teaching AI its Limits (BWST_IF2020-010). IP and DS
are supported by the U.S. Department of Energy (DOE), Office of Science grant
DOE-SC0010008 and HD, VM, and BN are supported by the DOE under contract DE-
AC02-05CH11231.

##  Appendix A Architecture and training

Here we briefly describe the architecture and training procedure used for
CaloFlow (see Refs. [42, 43] for more details). There are some differences
compared to the implementation in the original CaloFlow papers [42, 43], but
most of main algorithm remains the same. One main difference is that only
Flow-I is used in this study.

The flows used in this work are Masked Autoregressive Flows (MAFs) [92] with
compositions of Rational Quadratic Splines (RQS) [93] as transformations. The
RQS transformations are parameterized using neural networks known as MADE
blocks [94]. Identical flow architectures are used in each of the two cases
with uniformly and log-uniformly distributed $E_{\rm inc}$. Each flow consists
of six MADE blocks, each with two hidden layers of 64 nodes. The RQS
transformations are defined with 8 bins and a tail bound of 14.

The incident energy of the incoming photon is preprocessed as

$E_{\text{inc}}\to\log_{10}(E_{\text{inc}}/10\ \text{GeV})\,.$ (6)

The layer energies are preprocessed as

$E_{i}\to 2\left(\log_{10}(E_{i}+1\ \text{keV})-1\right)\,.$ (7)

The index $i$ denotes the layer number. In the original CaloFlow, a different
preprocessing was used for the layer energies $E_{i}$ in Flow-I where $E_{i}$
were transformed to unit-space (see [42]).

Uniform noise in the range [0,0.1] keV was applied to the voxel energies
during training and evaluation. The same range was used for the pion dataset
in [44]. The range of uniform noise used in the original CaloFlow [42, 43] was
[0,1] keV. The addition of noise was found to prevent the flow from fitting
unimportant features. The training of flows in this work is optimized using
independent Adam optimizers [87]. The flows were trained by minimizing $-\log
p_{1}(\vec{E}|E_{\text{inc}})$ for 150 epochs with a batch size of 200. The
initial learning of $10^{-4}$ was chosen and a multi-step learning schedule
was used when training the flow which halves the learning rate after each
selected epoch milestone during the training.

## References

  * [1] L. de Oliveira, B. Nachman and M. Paganini, _Electromagnetic Showers Beyond Shower Shapes_ , _Nucl. Instrum. Meth. A_ 951 (2020) 162879 [1806.05667].
  * [2] CMS collaboration, _The Phase-2 Upgrade of the CMS Endcap Calorimeter_ , .
  * [3] ATLAS Collaboration, _Deep Learning for Pion Identification and Energy Calibration with the ATLAS Detector_ , _ATL-PHYS-PUB-2020-018_ (2020) .
  * [4] C. Neubüser, J. Kieseler and P. Lujan, _Optimising longitudinal and lateral calorimeter granularity for software compensation in hadronic showers using deep neural networks_ , _Eur. Phys. J. C_ 82 (2022) 92 [2101.08150].
  * [5] N. Akchurin, C. Cowden, J. Damgov, A. Hussain and S. Kunori, _On the Use of Neural Networks for Energy Reconstruction in High-granularity Calorimeters_ , _JINST_ 16 (2021) P12036 [2107.10207].
  * [6] J. Kieseler, G.C. Strong, F. Chiandotto, T. Dorigo and L. Layer, _Calorimetric Measurement of Multi-TeV Muons via Deep Regression_ , _Eur.Phys.J.C_ 82 (2021) 79 [2107.02119].
  * [7] N. Akchurin, C. Cowden, J. Damgov, A. Hussain and S. Kunori, _Perspectives on the Calibration of CNN Energy Reconstruction in Highly Granular Calorimeters_ , 2108.10963.
  * [8] ATLAS collaboration, _Point Cloud Deep Learning Methods for Pion Reconstruction in the ATLAS Experiment_ , .
  * [9] S.R. Qasim, N. Chernyavskaya, J. Kieseler, K. Long, O. Viazlo, M. Pierini et al., _End-to-end multi-particle reconstruction in high occupancy imaging calorimeters with graph neural networks_ , _Eur.Phys.J.C_ 82 (2022) 753 [2204.01681].
  * [10] F.T. Acosta, B. Karki, P. Karande, A. Angerami, M. Arratia, K. Barish et al., _The Optimal use of Segmentation for Sampling Calorimeters_ , 2310.04442.
  * [11] ATLAS collaboration, _The application of neural networks for the calibration of topological cell clusters in the ATLAS calorimeters_ , .
  * [12] R. Haake and C. Loizides, _Machine Learning based jet momentum reconstruction in heavy-ion collisions_ , _Phys. Rev. C_ 99 (2019) 064904 [1810.06324].
  * [13] CMS collaboration, _A Deep Neural Network for Simultaneous Estimation of b Jet Energy and Resolution_ , _Comput. Softw. Big Sci._ 4 (2020) 10 [1912.06046].
  * [14] S. Cheong, A. Cukierman, B. Nachman, M. Safdari and A. Schwartzman, _Parametrizing the Detector Response with Neural Networks_ , _JINST_ 15 (2020) P01030 [1910.03773].
  * [15] ALICE collaboration, _Machine Learning based jet momentum reconstruction in Pb-Pb collisions measured with the ALICE detector_ , _PoS_ EPS-HEP2019 (2020) 312 [1909.01639].
  * [16] CMS collaboration, _Mass regression of highly-boosted jets using graph neural networks_ , .
  * [17] R. Gambhir, B. Nachman and J. Thaler, _Learning Uncertainties the Frequentist Way: Calibration and Correlation in High Energy Physics_ , _Phys.Rev.Lett._ 129 (2022) 082001 [2205.03413].
  * [18] R. Gambhir, B. Nachman and J. Thaler, _Bias and Priors in Machine Learning Calibrations for High Energy Physics_ , _Phys.Rev.D_ 106 (2022) 036011 [2205.05084].
  * [19] ATLAS Collaboration, _New techniques for jet calibration with the ATLAS detector_ , _Eur.Phys.J.C_ 83 (2023) 761 [2303.17312].
  * [20] ATLAS collaboration, _Simultaneous energy and mass calibration of large-radius jets with the ATLAS detector using a deep neural network_ , 2311.08885.
  * [21] ALICE collaboration, _Measurement of the radius dependence of charged-particle jet suppression in Pb-Pb collisions at $\sqrt{s_{\rm NN}}$ = 5.02 TeV_, 2303.00592.
  * [22] M. Diefenthaler, A. Farhat, A. Verbytskyi and Y. Xu, _Deeply Learning Deep Inelastic Scattering Kinematics_ , _Eur.Phys.J.C_ 82 (2021) 1064 [2108.11638].
  * [23] M. Arratia, D. Britzger, O. Long and B. Nachman, _Reconstructing the Kinematics of Deep Inelastic Scattering with Deep Learning_ , _Nucl.Instrum.Meth.A_ 1025 (2021) 166164 [2110.05505].
  * [24] M. Leigh, J.A. Raine and T. Golling, _$\nu$ -Flows: conditional neutrino regression_, _SciPost Phys._ 14 (2022) 159 [2207.00664].
  * [25] J.A. Raine, M. Leigh, K. Zoch and T. Golling, _$\nu^{2}$ -Flows: Fast and improved neutrino reconstruction in multi-neutrino final states with conditional normalizing flows_, 2307.02405.
  * [26] A. Cukierman and B. Nachman, _Mathematical Properties of Numerical Inversion for Jet Calibrations_ , _Nucl. Instrum. Meth. A_ 858 (2017) 1 [1609.05195].
  * [27] ATLAS Collaboration, _Simultaneous Jet Energy and Mass Calibrations with Neural Networks_ , Tech. Rep. ATL-PHYS-PUB-2020-001, CERN, Geneva (Jan, 2020).
  * [28] ATLAS Collaboration, _Generalized Numerical Inversion: A Neural Network Approach to Jet Calibration_ , Tech. Rep. ATL-PHYS-PUB-2018-013, CERN, Geneva (Jul, 2018).
  * [29] D.J. MacKay, _Probable networks and plausible predictions-a review of practical bayesian methods for supervised neural networks_ , _Network: computation in neural systems_ 6 (1995) 469.
  * [30] R.M. Neal, _Bayesian learning for neural networks_ , vol. 118, Springer Science & Business Media (2012).
  * [31] Y. Gal et al., _Uncertainty in deep learning_ , .
  * [32] G. Kasieczka, M. Luchmann, F. Otterpohl and T. Plehn, _Per-Object Systematics using Deep-Learned Calibration_ , 2003.11099.
  * [33] M. Bellagente, M. Haußmann, M. Luchmann and T. Plehn, _Understanding Event-Generation Networks via Uncertainties_ , _SciPost Phys._ 13 (2021) 003 [2104.04543].
  * [34] V. Mikuni and B. Nachman, _High-dimensional and Permutation Invariant Anomaly Detection_ , 2306.03933.
  * [35] T. Finke, M. Krämer, A. Mück and J. Tönshoff, _Learning the language of QCD jets with transformers_ , _JHEP_ 06 (2023) 184 [2303.07364].
  * [36] E.G. Tabak and C.V. Turner, _A family of nonparametric density estimation algorithms_ , _Communications on Pure and Applied Mathematics_ 66 (2013) 145.
  * [37] L. Dinh, D. Krueger and Y. Bengio, _Nice: Non-linear independent components estimation_ , _arXiv preprint arXiv:1410.8516_ (2014) .
  * [38] D.J. Rezende and S. Mohamed, _Variational inference with normalizing flows_ , 2015.  10.48550/ARXIV.1505.05770.
  * [39] L. Dinh, J. Sohl-Dickstein and S. Bengio, _Density estimation using real nvp_ , _arXiv preprint arXiv:1605.08803_ (2016) .
  * [40] G. Papamakarios, E. Nalisnick, D.J. Rezende, S. Mohamed and B. Lakshminarayanan, _Normalizing flows for probabilistic modeling and inference_ , .
  * [41] C. Winkler, D.E. Worrall, E. Hoogeboom and M. Welling, _Learning likelihoods with conditional normalizing flows_ , _CoRR_ abs/1912.00042 (2019) [1912.00042].
  * [42] C. Krause and D. Shih, _CaloFlow: Fast and Accurate Generation of Calorimeter Showers with Normalizing Flows_ , _Phys.Rev.D_ 107 (2021) 113003 [2106.05285].
  * [43] C. Krause and D. Shih, _CaloFlow II: Even Faster and Still Accurate Generation of Calorimeter Showers with Normalizing Flows_ , _Phys.Rev.D_ 107 (2021) 113004 [2110.11377].
  * [44] C. Krause, I. Pang and D. Shih, _CaloFlow for CaloChallenge Dataset 1_ , 2210.14245.
  * [45] M.R. Buckley, C. Krause, I. Pang and D. Shih, _Inductive CaloFlow_ , 2305.11934.
  * [46] S. Diefenbacher, E. Eren, F. Gaede, G. Kasieczka, C. Krause, I. Shekhzadeh et al., _L2LFlows: Generating High-Fidelity 3D Calorimeter Images_ , _JINST_ 18 (2023) P10017 [2302.11594].
  * [47] A. Xu, S. Han, X. Ju and H. Wang, _Generative Machine Learning for Detector Response Modeling with a Conditional Normalizing Flow_ , 2303.10148.
  * [48] I. Pang, J.A. Raine and D. Shih, _SuperCalo: Calorimeter shower super-resolution_ , 2308.11700.
  * [49] F. Ernst, L. Favaro, C. Krause, T. Plehn and D. Shih, _Normalizing Flows for High-Dimensional Detector Simulations_ , 2312.09290.
  * [50] B. Nachman and D. Shih, _Anomaly Detection with Density Estimation_ , _Phys. Rev. D_ 101 (2020) 075042 [2001.04990].
  * [51] C. Gao, J. Isaacson and C. Krause, _i-flow: High-Dimensional Integration and Sampling with Normalizing Flows_ , 2001.05486.
  * [52] E. Bothmann, T. Janßen, M. Knobbe, T. Schmale and S. Schumann, _Exploring phase space with Neural Importance Sampling_ , 2001.05478.
  * [53] C. Gao, S. Höche, J. Isaacson, C. Krause and H. Schulz, _Event Generation with Normalizing Flows_ , _Phys. Rev. D_ 101 (2020) 076002 [2001.10028].
  * [54] M. Bellagente, A. Butter, G. Kasieczka, T. Plehn, A. Rousselot, R. Winterhalder et al., _Invertible Networks or Partons to Detector and Back Again_ , _SciPost Phys._ 9 (2020) 074 [2006.06685].
  * [55] B. Stienen and R. Verheyen, _Phase space sampling and inference from weighted events with autoregressive flows_ , _SciPost Phys._ 10 (2021) 038 [2011.13445].
  * [56] S. Bieringer, A. Butter, T. Heimel, S. Höche, U. Köthe, T. Plehn et al., _Measuring QCD Splittings with Invertible Networks_ , _SciPost Phys._ 10 (2020) 126 [2012.09873].
  * [57] A. Hallin, J. Isaacson, G. Kasieczka, C. Krause, B. Nachman, T. Quadfasel et al., _Classifying Anomalies THrough Outer Density Estimation (CATHODE)_ , _Phys.Rev.D_ 106 (2021) 055006 [2109.00546].
  * [58] T. Bister, M. Erdmann, U. Köthe and J. Schulte, _Inference of cosmic-ray source properties by conditional invertible neural networks_ , _Eur.Phys.J.C_ 82 (2021) 171 [2110.09493].
  * [59] A. Butter, T. Heimel, S. Hummerich, T. Krebs, T. Plehn, A. Rousselot et al., _Generative Networks for Precision Enthusiasts_ , _SciPost Phys._ 14 (2021) 078 [2110.13632].
  * [60] R. Winterhalder, V. Magerya, E. Villa, S.P. Jones, M. Kerner, A. Butter et al., _Targeting Multi-Loop Integrals with Neural Networks_ , _SciPost Phys._ 12 (2022) 129 [2112.09145].
  * [61] A. Butter, S. Diefenbacher, G. Kasieczka, B. Nachman, T. Plehn, D. Shih et al., _Ephemeral Learning – Augmenting Triggers with Online-Trained Normalizing Flows_ , _SciPost Phys._ 13 (2022) 087 [2202.09375].
  * [62] R. Verheyen, _Event Generation and Density Estimation with Surjective Normalizing Flows_ , _SciPost Phys._ 13 (2022) 047 [2205.01697].
  * [63] A. Butter, T. Heimel, T. Martini, S. Peitzsch and T. Plehn, _Two Invertible Networks for the Matrix Element Method_ , _SciPost Phys._ 15 (2022) 094 [2210.00019].
  * [64] A. Hallin, G. Kasieczka, T. Quadfasel, D. Shih and M. Sommerhalder, _Resonant anomaly detection without background sculpting_ , _Phys.Rev.D_ 107 (2022) 114012 [2210.14924].
  * [65] T. Heimel, R. Winterhalder, A. Butter, J. Isaacson, C. Krause, F. Maltoni et al., _MadNIS – Neural Multi-Channel Importance Sampling_ , _SciPost Phys._ 15 (2022) 141 [2212.06172].
  * [66] M. Backes, A. Butter, M. Dunford and B. Malaescu, _An unfolding method based on conditional Invertible Neural Networks (cINN) using iterative training_ , 2212.08674.
  * [67] D. Sengupta, S. Klein, J.A. Raine and T. Golling, _CURTAINs Flows For Flows: Constructing Unobserved Regions with Maximum Likelihood Estimation_ , 2305.04646.
  * [68] J. Ackerschott, R.K. Barman, D. Gonçalves, T. Heimel and T. Plehn, _Returning CP-Observables to The Frames They Belong_ , 2308.00027.
  * [69] T. Heimel, N. Huetsch, R. Winterhalder, T. Plehn and A. Butter, _Precision-Machine Learning for the Matrix Element Method_ , 2310.07752.
  * [70] T. Heimel, N. Huetsch, F. Maltoni, O. Mattelaer, T. Plehn and R. Winterhalder, _The MadNIS Reloaded_ , 2311.01548.
  * [71] C. Bierlich, P. Ilten, T. Menzo, S. Mrenna, M. Szewc, M.K. Wilkinson et al., _Towards a data-driven model of hadronization using normalizing flows_ , 2311.09296.
  * [72] R. Das, G. Kasieczka and D. Shih, _Residual ANODE_ , 2312.11629.
  * [73] H. Du, C. Krause, V. Mikuni, B. Nachman, I. Pang and D. Shih, _Electromagnetic + Hadronic Sampling Calorimeter Shower Images_ , Apr., 2024.  10.5281/zenodo.11073232.
  * [74] B. Nachman, L. de Oliveira and M. Paganini, _Electromagnetic calorimeter shower images_ , _Mendeley Data_ (2017) .
  * [75] M. Paganini, L. de Oliveira and B. Nachman, _Accelerating Science with Generative Adversarial Networks: An Application to 3D Particle Showers in Multilayer Calorimeters_ , _Phys. Rev. Lett._ 120 (2018) 042003 [1705.02355].
  * [76] M. Paganini, L. de Oliveira and B. Nachman, _CaloGAN : Simulating 3D high energy particle showers in multilayer electromagnetic calorimeters with generative adversarial networks_ , _Phys. Rev._ D97 (2018) 014021 [1712.10321].
  * [77] C. Krause, B. Nachman, I. Pang, D. Shih and Y. Zhu, _Anomaly detection with flow-based fast calorimeter simulators_ , 2312.11618.
  * [78] G. Papamakarios, D. Sterratt and I. Murray, _Sequential neural likelihood: Fast likelihood-free inference with autoregressive flows_ , in _The 22nd international conference on artificial intelligence and statistics_ , pp. 837–848, PMLR, 2019.
  * [79] G. Papamakarios, _Neural density estimation and likelihood-free inference_ , _arXiv preprint arXiv:1910.13233_ (2019) .
  * [80] K. Cranmer, J. Brehmer and G. Louppe, _The frontier of simulation-based inference_ , _Proceedings of the National Academy of Sciences_ 117 (2020) 30055.
  * [81] S. Dirmeier, C. Albert and F. Perez-Cruz, _Simulation-based inference using surjective sequential neural likelihood estimation_ , _arXiv preprint arXiv:2308.01054_ (2023) .
  * [82] GEANT4 collaboration, _GEANT4–a simulation toolkit_ , _Nucl. Instrum. Meth. A_ 506 (2003) 250.
  * [83] J. Allison, K. Amako, J. Apostolakis, H. Araujo, P. Arce Dubois, M. Asai et al., _Geant4 developments and applications_ , _IEEE Transactions on Nuclear Science_ 53 (2006) 270.
  * [84] J. Allison, K. Amako, J. Apostolakis, P. Arce, M. Asai, T. Aso et al., _Recent developments in geant4_ , _Nuclear Instruments and Methods in Physics Research Section A: Accelerators, Spectrometers, Detectors and Associated Equipment_ 835 (2016) 186.
  * [85] ATLAS collaboration, _ATLAS liquid-argon calorimeter: Technical Design Report_ , Technical design report. ATLAS, CERN, Geneva (1996), 10.17181/CERN.FWRW.FOOQ.
  * [86] A. Butter, S. Diefenbacher, G. Kasieczka, B. Nachman and T. Plehn, _Amplifying statistics with ensembles of generative models_ , in _ICLR 2021 SimDL Workshop https://simdl. github. io/files/18. pdf_ , 2021.
  * [87] D.P. Kingma and J. Ba, _Adam: A method for stochastic optimization_ , 2014.
  * [88] P. Virtanen, R. Gommers, T.E. Oliphant, M. Haberland, T. Reddy, D. Cournapeau et al., _SciPy 1.0: Fundamental Algorithms for Scientific Computing in Python_ , _Nature Methods_ 17 (2020) 261.
  * [89] B. Efron, _Bootstrap Methods: Another Look at the Jackknife_ , _The Annals of Statistics_ 7 (1979) 1 .
  * [90] D. Scott, _Multivariate Density Estimation: Theory, Practice, and Visualization_ , John Wiley & Sons, New York, Chicester (1992).
  * [91] C.W. Fabjan and F. Gianotti, _Calorimetry for particle physics_ , _Rev. Mod. Phys._ 75 (2003) 1243.
  * [92] G. Papamakarios, T. Pavlakou and I. Murray, _Masked autoregressive flow for density estimation_ , _Advances in neural information processing systems_ 30 (2017) .
  * [93] C. Durkan, A. Bekasov, I. Murray and G. Papamakarios, _Neural spline flows_ , _Advances in Neural Information Processing Systems_ 32 (2019) 7511 [1906.04032].
  * [94] M. Germain, K. Gregor, I. Murray and H. Larochelle, _Made: Masked autoencoder for distribution estimation_ , in _International conference on machine learning_ , pp. 881–889, PMLR, 2015.


# Kinematics of the Local Spiral Structure Revealed by Young Stars in _Gaia_
DR3

Dejian Liu  Purple Mountain Observatory, Chinese Academy of Sciences, Nanjing
210023, People’s Republic of China  School of Astronomy and Space Science,
University of Science and Technology of China, Hefei 230026, People’s Republic
of China  Ye Xu  Purple Mountain Observatory, Chinese Academy of Sciences,
Nanjing 210023, People’s Republic of China  School of Astronomy and Space
Science, University of Science and Technology of China, Hefei 230026, People’s
Republic of China  Chaojie Hao  Purple Mountain Observatory, Chinese Academy
of Sciences, Nanjing 210023, People’s Republic of China  School of Astronomy
and Space Science, University of Science and Technology of China, Hefei
230026, People’s Republic of China  Shuaibo Bian  Purple Mountain Observatory,
Chinese Academy of Sciences, Nanjing 210023, People’s Republic of China  Zehao
Lin  Purple Mountain Observatory, Chinese Academy of Sciences, Nanjing 210023,
People’s Republic of China  School of Astronomy and Space Science, University
of Science and Technology of China, Hefei 230026, People’s Republic of China
Yingjie Li  Purple Mountain Observatory, Chinese Academy of Sciences, Nanjing
210023, People’s Republic of China  Jingjing Li  Purple Mountain Observatory,
Chinese Academy of Sciences, Nanjing 210023, People’s Republic of China

###### Abstract

Using young open clusters and O–B2-type stars in  _Gaia_ DR3, we investigate
the kinematics of the local spiral structure. In general, the young sources in
the outer spiral arms may present larger peculiar motions than those in the
inner spiral arms. The young open clusters appear to have smaller peculiar
motions than the O–B2-type stars, and the sources in both the Perseus and
Local Arms may show an inward motion toward the Galactic center and rotate
slower than Galactic rotation. Meanwhile, the sources in the Carina Arm may
move in the opposite direction from the Sun to the Galactic center and rotate
marginally faster than Galactic rotation. In addition, using young open
clusters and O–B2-type stars, we have improved the distance estimations of
kinematic methods for several regions near the Sun.

Galaxy kinematics (602); Open star clusters (1160); OB stars (1141)

,

##  1 Introduction

The kinematics of the spiral arms of the Milky Way are of great interests to
astronomers, as they can help unveil the formation and evolution of our home
Galaxy. As key elements of the spiral arms, the measurements and analysis of
the spatial distributions, velocities, and rotational orbits of stars and gas
therein can provide valuable insights into Galactic structure and dynamics
(Binney & Tremaine, 2008), as well as the gravitational potential of the Milky
Way (Sakai et al., 2019; Immer et al., 2019).

Numerous young objects have been used to trace the spiral-arm structure of our
Galaxy, e.g., high-mass star-formation region (HMSFR) masers yielded by the
Bar and Spiral Structure Legacy (BeSSeL) Survey (Brunthaler et al., 2011; Reid
et al., 2014, 2019) and the Very Long Baseline Interferometry (VLBI)
Exploration of Radio Astrometry (VERA) array (VERA Collaboration et al.,
2020), as well as observations of young massive O–B-type stars (Xu et al.,
2018, 2021), Hii regions (Georgelin & Georgelin, 1976; Hou & Han, 2014), young
open clusters (YOCs, Dias et al., 2019; Hao et al., 2021; Gaia Collaboration
et al., 2022a), etc. Although the global Galactic spiral structure is still a
topic of debate, the spiral structure within about 5 kpc of the Sun is better
understood. There are three spiral segments in the solar neighborhood, namely
the Sagittarius–Carina, Local, and Perseus Arms (Reid et al., 2019). Recently,
Xu et al. (2023) reconsidered the spiral structure of the Milky Way and
proposed that the Sagittarius Arm may be not connected to the Carina Arm. Due
to a lack of available sources in the Sagittarius Arm, in this work, we mainly
concentrate on the three distinct spiral-arm segments in the solar
neighborhood: the Carina, Local, and Perseus Arms.

The kinematics of the Milky Way are closely related to its structure. In
recent years, substantial progress in studying the local kinematics using
stars provided by _Gaia_ and VLBI masers has been achieved. Using Gaia data
release 2 (Gaia Collaboration et al., 2018a), Gaia Collaboration et al.
(2018b) looked at the kinematics of the Milky Way disk, within a radius of
several kiloparsecs around the Sun, and found rich and complex kinematics of
the Galactic disk and streaming motions in all three dimensional velocity
components. Gaia Collaboration et al. (2022a) studied the Galactic spiral
structure and kinematics using young stellar populations in _Gaia_ data
release 3 (DR3, Gaia Collaboration et al., 2022b), and revealed that the
velocity field of the OB stars shows streaming motions whose characteristic
length is similar to the spiral arm density. As reviewed by Immer & Rygl
(2022), based on VLBI masers, researchers have studied the kinematics of
sources located in different spiral arm segments. However, many of these
studies focused only on the peculiar motion of masers in individual spiral arm
segments, while systematic kinematic studies of the local spiral structure,
especially using young sources located in the spiral arms, are still very
lacking. In this work, based on the spiral-arm model proposed by Xu et al.
(2023), we aim to study the kinematics of the local spiral structure using
young stars.

Accurate trigonometric parallaxes and proper motions of HMSFR masers in the
northern hemisphere sky have been determined (e.g., Reid et al., 2019; Xu et
al., 2021; Bian et al., 2022). The high-precision six-dimensional (6D; 3D
spatial and 3D velocity) astrometric parameters of OB stars and YOCs, provided
by  _Gaia_ DR3, enable us to investigate the kinematics of young stars in the
local spiral structure, and compare them with those of HMSFR masers.
Additionally, based on the astrometric parameters of young stars in the solar
neighborhood, we try to provide valuable indications for future distance
estimation using kinematic methods.

##  2 Data

###  2.1 YOCs

Xu et al. (2023) compiled a sample that contains more than 5 000 open clusters
(OCs) from three catalogs (Hao et al., 2021, 2022b; Castro-Ginard et al.,
2022). In this work, we increase this catalog by considering the 38 newly
reported OCs by Hao et al. (2022a) and the 101 newly reported OCs by Qin et
al. (2023). We also identify potentially duplicate OCs. If the difference
between the mean parameters of two OCs are within 3$\sigma_{i}$
($i=\alpha,\delta,\varpi,\mu_{\alpha}\cos\delta,\mu_{\delta}$) in five-
dimensional parametric space, they can be regarded as duplicate objects, and
we confirm them via visual inspection. After this step, 4 871 OCs remain (see
Table 1), of which 3 794 have _Gaia_ radial velocity measurement. In the
following analysis, we only use young OCs (age $<$ 20 Myr) with parallax
accuracies better than 20% and radial velocity uncertainties smaller than 10
${\rm km~{}s}^{-1}$. Determination of the radial velocity uncertainties of OCs
is conducted following the same method as Soubiran et al. (2018). Finally, we
obtain a sample containing 285 YOCs, and their astrometric parameters are
taken from _Gaia_ DR3. We perform a parallax zero-point correction for the
member stars of each OC (Lindegren et al., 2021), and then use the average
parallax of all members as the OC’s parallax.

Table 1: Summary of the Mean Parameters of the Cluster Sample

cluster | $N$ | $\alpha$ | $\delta$ | $\varpi$ | $\sigma_{\varpi}$ | $\mu_{\alpha}\cos\delta$ | $\sigma_{\mu_{\alpha}\cos\delta}$ | $\mu_{\delta}$ | $\sigma_{\mu_{\delta}}$ | $v_{\rm los}$ | $\sigma_{v_{\rm los}}$ | $N_{v_{\rm los}}$ | $\log t$ | Ref
---|---|---|---|---|---|---|---|---|---|---|---|---|---|---
|  | (deg) | (deg) | (mas) | (mas) | (mas yr-1) | (mas yr-1) | (mas yr-1) | (mas yr-1) | (${\rm km~{}s}^{-1}$) | (${\rm km~{}s}^{-1}$) |  |  |
UBC1194 | 64 | 0.020 | 63.900 | 0.330 | 0.023 | -2.496 | 0.024 | -1.261 | 0.023 | -77.9 | 10.7 | 5 | 7.78 | Castro2022
Berkeley_58 | 115 | 0.063 | 60.937 | 0.331 | 0.084 | -3.361 | 0.091 | -1.705 | 0.089 | -4.1 | 10.0 | 1 | 8.47 | Scholz2015
NGC_7801 | 33 | 0.077 | 50.711 | 0.444 | 0.064 | -3.642 | 0.056 | -2.439 | 0.075 | -12.5 | 14.4 | 4 | 9.26 | Scholz2015
… |  |  |  |  |  |  |  |  |  |  |  |  |  |

Note. — This table is available in its entirety in machine-readable form.

###  2.2 O–B2-type Stars

We adopt the method of Xu et al. (2023) to identify O–B2-type stars in _Gaia_
DR3 and a renormalized unit weight error (ruwe) $<1.4$ (Lindegren et al.,
2021) is used to remove objects with unreliable astrometric solutions. After
that, 1 135 O–B2-type stars with _Gaia_ radial velocities are selected (see
Table 2). The values of the parameter ipd_gof_harmonic_amplitude above 0.1 in
combination with ruwe being larger than 1.4 are indicative of resolved
binaries, which are still not correctly handled in the astrometric processing,
and may cause spurious solutions (Gaia Collaboration et al., 2021). Hence,
these two parameters are adopted to reduce the number of potential binaries,
and then reduce the effect due to them, i.e., ipd_gof_harmonic_amplitude
$<0.1$ and ruwe $<1.4$ (Fabricius et al., 2021). Applying the same criteria of
a parallax precision better than 20% and a radial velocity error smaller than
10 ${\rm km~{}s}^{-1}$, we obtained a final sample of 300 O–B2-type stars.
Same as the OCs, we also perform a parallax zero-point correction for each
O–B2-type star in the sample (Lindegren et al., 2021).

Table 2: Summary of the Parameters of the O–B2-type Stars

_Gaia_ DR3 ID | $\alpha$ | $\delta$ | $\varpi$ | $\sigma_{\varpi}$ | $\mu_{\alpha}\cos\delta$ | $\sigma_{\mu_{\alpha}\cos\delta}$ | $\mu_{\delta}$ | $\sigma_{\mu_{\delta}}$ | $v_{\rm los}$ | $\sigma_{v_{\rm los}}$ | $T_{\rm eff}$
---|---|---|---|---|---|---|---|---|---|---|---
| (deg) | (deg) | (mas) | (mas) | (mas yr-1) | (mas yr-1) | (mas yr-1) | (mas yr-1) | (${\rm km~{}s}^{-1}$) | (${\rm km~{}s}^{-1}$) | (K)
528594342521399168 | 0.4453 | 67.5070 | 1.010 | 0.012 | -1.566 | 0.011 | -1.795 | 0.012 | -42.7 | 19.3 | 42190
429366957178085376 | 1.0123 | 60.8458 | 0.303 | 0.023 | -2.544 | 0.018 | -1.969 | 0.021 | -51.0 | 34.9 | 24480
431544127639952384 | 1.9834 | 63.3354 | 0.427 | 0.018 | -2.823 | 0.021 | -0.925 | 0.020 | -55.5 | 6.9 | 20000
… |  |  |  |  |  |  |  |  |  |  |

Note. — This table is available in its entirety in machine-readable form.

##  3 Kinematics of Young Objects

Using the parallaxes to convert proper motions to linear velocities and
combining them with the radial velocities, we derive the peculiar motions
($U_{s},V_{s},{\rm and}~{}W_{s}$) of the YOCs and O–B2-type stars, using the
same method as Reid et al. (2009). $U_{s},V_{s}$, and $W_{s}$ are the velocity
components towards the Galactic center (GC, radial), in the direction of
Galactic rotation ( azimuthal), and towards the north Galactic pole,
respectively. In the above calculation, the fundamental parameters of the
Milky Way, $R_{0}=8.15$ kpc, $\Theta_{0}=236$ ${\rm km~{}s}^{-1}$, and solar
motions of $U_{\odot}$ = 10.6 ${\rm km~{}s}^{-1}$, $V_{\odot}$ = 10.7 ${\rm
km~{}s}^{-1}$, and $W_{\odot}$ = 7.6 ${\rm km~{}s}^{-1}$ (Reid et al., 2009)
are adopted. Meanwhile, a universal rotation curve (Persic et al., 1996)
consisting of two parameters, a2 = 0.96 and a3 = 1.62 (Reid et al., 2019), is
assumed.

Figure 1 displays the peculiar motions of the YOCs and O–B2-type stars in the
Galactic plane, which are almost concordant with each other. There are two
significant kinematic features of sources in the Perseus and Local Arms. Near
the direction of $l$ $\sim$ $100^{\circ}$–$135^{\circ}$, the YOCs and
O–B2-type stars in the Perseus Arm exhibit large peculiar motions, which is
consistent with the result shown using HMSFR masers (Reid et al., 2019).
Around the direction of $l~{}\sim 70^{\circ}$–$90^{\circ}$, both the YOCs and
O–B2-type stars in the Local Arm exhibit opposite motions with respect to
Galactic rotation, similar to the results of HMSFR masers (Xu et al., 2013).
Therefore, these results suggest that the peculiar motions of the different
young populations in the solar neighborhood probably are comparable.

To address the issue of whether there are kinematic differences between young
stars in different evolutionary stages, we divide the O–B2-type stars into the
O-type stars and B0–B2-type stars based on their effective temperature, where
the O-type stars are hotter than 30 000 K and B0–B2-type stars are hotter than
20 000 K but colder than 30 000 K (Chen et al., 2013). As shown in Figure 2,
we present the peculiar motions of the O- and B0–B2-type stars in the Galactic
plane. It appears that the B0–B2-type stars, overall, possess larger peculiar
motions than those of the O-type stars.

Next, we investigate the kinematic characteristics of the YOCs and O–B2-type
stars in each spiral arm near the Sun, and analyze the average peculiar motion
($\overline{U_{s}}$, $\overline{V_{s}}$, and $\overline{W_{s}}$) of those
sources with respect to the Galactocentric radius $R$. Following Reid et al.
(2019), we define an outlier as having greater than a 3$\sigma$ residual in
any peculiar motion component and, based on this, we remove 20 YOCs and 44
O–B2-type stars from the following analysis.

(a) YOCs (b) O–B2-type stars

Figure 1: Peculiar motions of the YOCs (left) and O–B2-type stars (right) in
the Galactic plane. The parallax uncertainties are smaller than 20% and the
radial velocity uncertainties are smaller than 10 ${\rm km~{}s}^{-1}$. The
solid curved lines trace the centers of the spiral arms: from top to bottom
are the Perseus, Local, and Carina Arms, respectively (Xu et al., 2023). The
shadows and dashed lines are the $1\sigma$ and $2\sigma$ widths from the
center of the spiral arms, respectively. The Sun (black Sun symbol) is at (0,
8.15) kpc (Reid et al., 2019). A 50 ${\rm km~{}s}^{-1}$ scale vector is shown
in the lower-left corner.

(a) O-type stars (b) B0–B2-type stars

Figure 2: Peculiar motions of the O (left) and B0–B2 (right) stars in Galactic
plane. The others features are the same as in Figure 1.

###  3.1 Kinematics of Young Objects in the Spiral Arms

In this work, we focus on analyzing the kinematic characteristics of sources
in the three spiral arms near the Sun, including the Perseus, Local, and
Carina Arms. Besides, by comparing the kinematics of sources at different
distances (i.e., $1\sigma$ and $2\sigma$) from the center of each spiral arms,
we can investigate the areas most influenced by the gravitational force of
each spiral arms. Adopting the spiral-arm model proposed by Xu et al. (2023),
the average peculiar motions and std values are similar for sources located
within $1\sigma$ and $2\sigma$ widths from the center of each spiral arm, so
we just use the $2\sigma$ for this study. In Figures 3 and 4, we show
histograms of peculiar motions of YOCs and O–B2-type stars. Then, we calculate
the average peculiar motions and corresponding errors, and the standard
deviation (std) of the peculiar motions of sources in each spiral arm, which
are listed in Tables 3 and 4. For the 199 HMSFR masers summarized by Reid et
al. (2019), and the five newly reported by Xu et al. (2021) and Bian et al.
(2022), Xu et al. (2023) re-assigned them to the spiral arms in their new
proposed spiral-arm model. Here, the masers belonging to the Perseus Arm in
the second and third quadrants, the Local Arm, and the Carina Arm are selected
and their average peculiar motions are calculated and compared with those of
the YOCs and O–B2-type stars, as listed in Tables 3 and 4. Due to the limited
number of samples, the average peculiar motions from Table 3 are not
significant at a $3\sigma$ level, but they still reflect the potential
kinematic characteristics of young objects located in the spiral arms.

As shown in Figure 1, in the Perseus Arm, 18 YOCs and 40 O–B2-type stars are
located in the second quadrant, and 13 YOCs and 11 O–B2-type stars are located
in the third quadrant. This distribution is consistent with the fact that
there are fewer and generally faint clouds in the third quadrant of the
Perseus Arm (Dame et al., 2001). 75% YOCs and 73% O–B2-type stars have
positive values for the $U_{s}$ component, and 81% YOCs and 76% O–B2-type
stars have negative values for the $V_{s}$ component. As listed in Table 3,
the YOCs and O–B2-type stars in the Perseus Arm have positive average values
of the $U_{s}$ component and negative average values of the $V_{s}$ component,
indicating that the young objects are likely to have a tendency to move
towards the GC and move slower than Galactic rotation. It is in line with the
previous result obtained from masers (Sakai et al., 2019; Zhang et al., 2019).
Besides, on average, the O–B2-type stars may exhibit larger peculiar motions
than the masers and YOCs. When separating the O- and B0–B2-type stars, we
discover that the $U_{s}$ and $V_{s}$ components of the peculiar motions of
the O-type stars may be smaller than those of the B0–B2-type stars, while the
uncertainties of O-type stars are large.

In the Local arm, 64 YOCs and 44 O–B2-type stars are located in the second
quadrant, and 76 YOCs and 8 O–B2-type stars in the third quadrant. 57% YOCs
and 75% O–B2-type stars have positive values for the $U_{s}$ component, 53%
YOCs and 71% O–B2-type stars have negative values for the $V_{s}$ component.
Combining with the average peculiar motions of young objects (listed in Table
3), it is likely that O–B2-type stars, on average, present a movement toward
the GC and a lagging motion with respect to the Galactic rotation. YOCs may
have the same kinematic features as O–B2-type stars, but they are not
significant. The sources in the Local Arm may have smaller motions toward the
GC than those of sources in the Perseus Arm. The average peculiar motions of
the masers in the Local Arm may be larger than those of the YOCs, but smaller
than those of the O–B2-type stars in the $V_{s}$ component. The B0–B2-type
stars seem to present larger peculiar motions than the O-type stars.

In the Carina arm, 51% YOCs and 61% O–B2-type stars have negative values of
the $U_{s}$ component, 51% YOCs and 59% O–B2-type stars have positive values
of the $V_{s}$ component. The O–B2-type stars in the Carina Arm may have
different kinematic features, in that they appear to be moving towards the
Galactic anticenter. The O–B2-type stars appear to have no lagging motion and
may even be ahead of Galactic rotation. The average values for the $U_{s}$ and
$V_{s}$ components are $0.3\pm 2.0$ ${\rm km~{}s}^{-1}$ and $-2.3\pm 1.9$
${\rm km~{}s}^{-1}$. Thus, the YOCs appear to have no peculiar motion in the
radial direction and co-rotate with the Galactic orbits. The average $V_{s}$
component of the masers in the Carina Arm is close to that of the young stars,
and they may not exhibit lagging motion in relative to Galactic rotation.
Contrary to the features in the Perseus and Local Arms, the peculiar motions
of the O-type stars in the Carina Arm may be larger than those of the
B0–B2-type stars.

The average peculiar motions and std values of the O–B2-types may be larger
than those of YOCs or masers in the spiral arms near the Sun. Additionally, in
each spiral arm, the B0–B2-type stars perhaps show larger average peculiar
motions than the O-type stars. Interestingly, the peculiar motions of sources
located in the Perseus and Local Arms may be larger than those in the Carina
Arm, implying that the peculiar motions of the young objects at different
Galactocentric radii may be different. This issue will be investigated in the
next subsection.

Figure 3: Histograms of peculiar motions of the YOCs.

Figure 4: Histograms of peculiar motions of the O–B2 type stars. Table 3: Average Peculiar Motions of the Sources in the Individual Spiral Arms  Arm | Component | YOCs | O–B2 | O | B0–B2 | Maser
---|---|---|---|---|---|---
Perseus |  $\overline{U_{s}}$ (${\rm km~{}s}^{-1}$) | $6.0\pm 2.6$ | $8.8\pm 3.1$ | $-3.8\pm 6.7$ | $13.1\pm 3.1$ | $8.9\pm 2.4$
$\overline{V_{s}}$ (${\rm km~{}s}^{-1}$) | $-6.5\pm 1.6$ | $-8.5\pm 2.3$ | $-4.0\pm 6.2$ | $-10.0\pm 2.1$ | $-8.0\pm 2.5$
$\overline{W_{s}}$ (${\rm km~{}s}^{-1}$) | $0.3\pm 0.7$ | $-0.7\pm 0.7$ | $-1.2\pm 2.0$ | $-0.5\pm 0.6$ | $0.7\pm 1.7$
Num | 31 | 51 | 13 | 38 | 28
Local |  $\overline{U_{s}}$ (${\rm km~{}s}^{-1}$) | $1.8\pm 0.7$ | $6.5\pm 1.5$ | $3.1\pm 4.0$ | $7.3\pm 1.5$ | $-0.4\pm 1.8$
$\overline{V_{s}}$ (${\rm km~{}s}^{-1}$) | $-1.8\pm 0.7$ | $-12.2\pm 2.7$ | $-6.8\pm 6.9$ | $-13.3\pm 2.9$ | $-7.9\pm 1.0$
$\overline{W_{s}}$ (${\rm km~{}s}^{-1}$) | $0.3\pm 0.3$ | $-0.7\pm 0.7$ | $0.8\pm 1.8$ | $-1.0\pm 0.7$ | $1.8\pm 1.6$
Num | 140 | 52 | 9 | 43 | 28
Carina |  $\overline{U_{s}}$ (${\rm km~{}s}^{-1}$) | $0.3\pm 2.0$ | $-4.8\pm 2.0$ | $-9.0\pm 3.7$ | $-2.6\pm 2.4$ | $3.7\pm 1.7$
$\overline{V_{s}}$ (${\rm km~{}s}^{-1}$) | $-2.3\pm 1.9$ | $4.5\pm 1.8$ | $-3.3\pm 3.2$ | $8.5\pm 1.9$ | $0.4\pm 1.4$
$\overline{W_{s}}$ (${\rm km~{}s}^{-1}$) | $0.0\pm 0.8$ | $-1.7\pm 0.5$ | $-1.7\pm 1.1$ | $-1.6\pm 0.6$ | $-4.2\pm 1.2$
Num | 35 | 85 | 29 | 56 | 17
Table 4: Std Values of the Peculiar Motions of Sources in the Individual Spiral Arms  Arm | Component | YOCs | O–B2 | O | B0–B2 | Maser
---|---|---|---|---|---|---
Perseus |  $\sigma_{U_{s}}$ (${\rm km~{}s}^{-1}$) | 14.7 | 21.8 | 24.2 | 19.1 | 13.0
$\sigma_{V_{s}}$ (${\rm km~{}s}^{-1}$) | 8.8 | 16.3 | 22.5 | 13.2 | 13.2
$\sigma_{W_{s}}$ (${\rm km~{}s}^{-1}$) | 3.7 | 4.7 | 7.0 | 3.5 | 9.0
Local |  $\sigma_{U_{s}}$ (${\rm km~{}s}^{-1}$) | 8.3 | 10.6 | 11.9 | 10.1 | 9.4
$\sigma_{V_{s}}$ (${\rm km~{}s}^{-1}$) | 8.7 | 19.5 | 20.6 | 19.1 | 5.5
$\sigma_{W_{s}}$ (${\rm km~{}s}^{-1}$) | 4.1 | 4.7 | 5.3 | 4.5 | 8.4
Carina |  $\sigma_{U_{s}}$ (${\rm km~{}s}^{-1}$) | 11.7 | 18.8 | 19.9 | 17.9 | 6.9
$\sigma_{V_{s}}$ (${\rm km~{}s}^{-1}$) | 11.2 | 16.4 | 17.4 | 14.3 | 5.8
$\sigma_{W_{s}}$ (${\rm km~{}s}^{-1}$) | 4.9 | 4.7 | 5.7 | 4.1 | 5.1

###  3.2 Peculiar Motions Variations along Galactocentric Radius

To study the peculiar motions of the YOCs and O–B2-type stars with respect to
Galactocentric radius, we analyze the average peculiar motion as a function of
the Galactocentric radius in the range 6–11 kpc, as presented in Figures 5–6.
Since sources located in the two spiral arms may have the same Galactocentric
distance, we show the average peculiar motions variations along the
Galactocentric radius for all sources and sources located in different spiral
arms, respectively. Here, the average peculiar motion of each type of sources
is calculated within a window of 1 kpc width in steps of 0.25 kpc.

In these figures, the undulations of the peculiar motions of the O–B2-type
stars as a function of Galactocentric radius are stronger than those displayed
by the YOCs. The $U_{s}$ component of YOCs may have two “peak”, which are
located at $R\sim 7.0$ kpc and $R\sim 9.0$ kpc. However, considering the
uncertainties of the values, the peaks are not significant. Thus, more YOCs
are needed to further verify them. The $U_{s}$ component of the O–B2-type
stars increases from $R\sim$ 7.0 kpc and reaches a maximum of $\sim$ 14 ${\rm
km~{}s}^{-1}$ at a distance of $R\sim 9.0$ kpc. For the $V_{s}$ components of
the YOCs, they show a flat distribution from 6.0 to 9.0 kpc, and decrease
beyond 9.0 kpc. For the $V_{s}$ components of the O–B2-type stars, they
decline from 7.0 to 8.5 kpc and reach minimum of $\sim-12$ ${\rm km~{}s}^{-1}$
at a distance of $R\sim 8.5$ kpc. Outside 7.5 kpc, the $V_{s}$ component of
the YOCs and O–B2-type stars remain smaller than zero, indicating the young
objects may rotate slower than Galactic rotation. The YOCs are more massive
than single O-B2-type stars and would therefore better represent the average
motion. Their declining $V_{s}$ velocity as a function of radius beyond 9.0
kpc could also suggest that the assumed rotation curve may overestimete the
rotation velocity of the Milky Way. As for the $W_{s}$ components of the YOCs
and O–B2-type stars, their undulations are both weak and do not exceed
$\pm~{}2$ ${\rm km~{}s}^{-1}$ in the Galactocentric range of 6–11 kpc.

In the Perseus Arm, the $U_{s}$ components appear to decrease with the
Galactocentric radius for both YOCs and O–B2-type stars. The $V_{s}$
components of YOCs appear to decrease with Galactocentric radius, but increase
with the radius for the O–B2-type stars. In the Local Arm, both the $U_{s}$
and $V_{s}$ components appear to show a “U” structure along the Galactocentric
radius for YOCs and O–B2-type stars. The lowest point of “U” structure is at
about 8.0 kpc. In the Carina Arm, the $U_{s}$ components of YOCs and O–B2-type
stars appear to increase as the Galactocentric radius and would pass through
the zero point at 7.0–7.5 kpc. The $V_{s}$ components of YOCs and O–B2-type
stars do not appear to fluctuate significantly with the Galactocentric radius.

Between 7.0 and 8.0 kpc, the sources located in the Carina Arm and the Local
Arm may have the same Galactocentric radius. We further investigate the
kinematic features of sources located in different spiral arms while having
the same Galactocentric radius. The $U_{s}$ components decrease in this range
for sources located in the Local Arm, while they increase for sources located
in the Carina Arm. They would intersect at around 7.5 kpc. The $V_{s}$
components of the YOCs located in the Local Arm are larger than those in the
Carina Arm. O–B2-type stars seem to exhibit the opposite phenomenon of YOCs.
Both young objects indicate that the kinematic features may be different for
the sources located in the two arms, while it needs to be confirmed with more
young stars.

Figure 5: The average peculiar motions of the YOCs as a function of
Galactocentric distance, $R$. The black solid lines indicate all sources, and
shaded area represent the uncertainties. The magenta, blue, and green lines
represent the sources located in the Perseus, Local, and Carina Arms. The
average is calculated within a window of 1 kpc width in steps of 0.25 kpc, and
$R_{0}$ is 8.15 kpc. Each interval must have more than five sources.

Figure 6: The average peculiar motions of the O–B2-type stars as a function of
Galactocentric distance, $R$. The others are the same as Figure 5.

##  4 Kinematic Distance

Traditionally, one can construct a simple model of the rotation speeds of
sources as a function of distance from the center of the Milky Way. Then, if
one measures the line-of-sight velocity (radial velocity) of an object, one
can determine its distance by matching the observation with the model
prediction, which yields a kinematic distance. However, the distances of
sources derived from kinematic methods may be unreliable due to noncircular
(peculiar) motions. The significant peculiar motions of the sources are
expected to affect their estimated kinematic distance. Xu et al. (2006) found
a parallax distance of $1.95\pm 0.04$ kpc of the high-mass star-forming region
W3(OH). The kinematic distance (i.e., $>$ 4 kpc) was more than twice as far as
the parallax distance. Hence, it is crucial to correct the kinematic model in
regions where the peculiar motions of sources are prominent to obtain
relatively precise kinematic distance estimations.

First, we convert the heliocentric radial velocities of the YOCs and O–B2-type
stars to the local standard of rest (LSR) velocities. The solar motion values
derived from stars in different evolutionary stages are different (Binney &
Tremaine, 2008). Since the objects discussed here are mainly young stars, we
do not use the IAU solar motion. Instead, we use a “revised LSR velocity”
($V_{\rm LSR}$) by applying the updated solar motion values obtained from
masers (Reid et al., 2019) to estimate the kinematic distance of each object,
which can be calculated by

$V_{\rm LSR}=V_{\rm Helio}+(V_{\odot}\sin l+U_{\odot}\cos l)\cos
b+W_{\odot}\sin b,$ (1)

where $V_{\rm Helio}$ is the heliocentric radial velocity, $l$ and $b$ are the
Galactic longitude and latitude, and $U_{\odot}$, $V_{\odot}$, and $W_{\odot}$
are the solar motions. For a source with a given Galactic longitude and
distance from the Sun, we can calculate its circular motion, $\Theta$, in the
Galactic plane using the universal rotation curve (Persic et al., 1996; Reid
et al., 2019).

The radial velocity of a source is the key to determining its kinematic
distance. Here, we define an assumed LSR radial velocity in a pure Galactic
circular orbit, $V_{\rm LSR}^{\rm mod}$, which can be calculated by

$V_{\rm LSR}^{\rm mod}=(\Theta(r)\cos\gamma-\Theta(R_{0})\sin l)*\cos b,$ (2)

where $\Theta(r)$ is the theoretical rotation velocity at Galactocentric
radius $r$ calculated by using the rotation curve. $\gamma$ is the angle
between the line-of-sight direction and the rotation direction of the source,
and $R_{0}=8.15$ kpc (Reid et al., 2019). A schematic depiction of these
parameters is given in Figure 7. In contrast, if we know the LSR velocity of a
source, its kinematic distance, $d$, can be derived by Equation (2).

Figure 7: Schematic depiction of the source and Galactic parameters.

We adopt the probability density function (PDF, $P$) to estimate the kinematic
distance. Following Reid et al. (2016), we calculate it via

$P=\frac{1-\exp[-(\Delta V_{\rm LSR}/\sigma_{V_{\rm LSR}})^{2}/2]}{(\Delta
V_{\rm LSR}/\sigma_{V_{\rm LSR}})^{2}/2},$ (3)

where $\Delta V_{\rm LSR}$ is the difference between the LSR velocity and that
expected from the rotation model of the Galaxy, $\sigma_{V_{\rm LSR}}$ is the
uncertainty of the LSR velocity. Uncertainties in the kinematic distances
include the 1$\sigma$ width of the PDF and uncertainties caused by the
parameters in the rotation curve, i.e., a2 ($0.96\pm 0.05$) and a3 ($1.62\pm
0.02$, Reid et al., 2019).

Figure 8 (a) presents the distributions of YOCs, O–B2-type stars, and masers
in a longitude–velocity ($l$–$v$) diagram. Since YOCs and O–B2-type stars are
mainly located within 6 kpc from the Sun, only masers in this range are
selected for comparison. The distributions of the YOCs and O–B2-type stars are
roughly consistent with that of the masers in the $l$–$v$ diagram, although
they all have large dispersions. We calculate the average LSR radial
velocities of the masers, YOCs, and O–B2-type stars with a Galactic longitude
step every 20∘, as shown Figure 8 (b). Overall, the LSR velocities of the
three young populations are concordant in most locations. Near the Galactic
longitude of $30^{\circ}$, the large differences in the LSR velocities among
the three groups may be due to a lack of YOCs and O–B2-type stars (see Figure
8 (a)). Therefore, we can combine the YOCs and O–B2-type stars to revise the
model of kinematic distance estimation.

In Figure 9, we display the distributions of the three spiral arms near the
Sun in the $l$–$v$ diagram and of the YOCs and O–B2-type stars located within
$1\sigma$ from the centers of the spiral arms in the Galactic plane for
comparison. Here, the spiral arms are from the model presented by Xu et al.
(2023), and $V_{\rm LSR}^{\rm mod}$ of each spiral arm is calculated using
Equation (2). In the $l$–$v$ diagram, the selected YOCs and O–B2-type stars
are expected located near the spiral arms, while there are significant and
systematic discrepancies between their measured $V_{\rm LSR}$ (using Equation
(1), shown as the colored points in Figure 9) and $V_{\rm LSR}^{\rm mod}$
values (the colored lines in Figure 9). It is caused by the peculiar motion of
sources. This suggests that the kinematic distances of these objects directly
obtained from radial velocities deviate from their parallax distances.
Therefore, in order to obtain more reliable kinematic distances for these
sources, one can consider deducting the peculiar motions by adding offsets to
radial velocities during the estimation.

Considering that the amount of YOCs and O–B2-type stars is limited and that
they are mainly concentrated in spiral arms, we focus on correcting the
kinematic distance model toward the regions where the LSR velocity
discrepancies between the sources and spiral arms are significant. For
example, near the direction of $l\sim 60^{\circ}-90^{\circ}$ in the Local Arm,
sources with LSR velocities less than $-10$ ${\rm km~{}s}^{-1}$ have
significant discrepancies between $V_{\rm LSR}$ and $V_{\rm LSR}^{\rm mod}$,
as shown in Figure 9. After an offset of approximately 28 ${\rm km~{}s}^{-1}$
is applied to those sources, their distributions are close to the Local Arm in
the $l$–$v$ diagram. Hence, this offset can yield more reliable kinematic
distance estimations for these sources. In total, six regions are found to
need compensation via LSR velocity offsets for accurate kinematic distance
estimations (see the expected values of these offsets in Table 5). Although
the offsets are determined using sources in the spiral arms, in practice, it
is not necessary to consider whether the source is located in a spiral arm,
but only the conditions listed in Table 5, i.e., Galactic longitude and radial
velocity.

The offsets are derived from the YOCs and O–B2-type stars located within
1$\sigma$ from the center of the spiral arms in the Galactic plane. Then we
apply them to all young objects that meet the criteria in Table 5. The
parallax distances, standard distances, and revised distances of the YOCs and
O–B2-type stars are presented in Figure 10(a). Here, the standard kinematic
distances are estimated without considering peculiar motion, while the revised
kinematic distances are obtained by adding the offsets listed in Table 5 to
counteract the peculiar motion. For the YOCs and O–B2-type stars, their
standard kinematic distances are generally greater than their parallax
distances, while their revised kinematic distances are closer to their
parallax distances. Moreover, we divide the YOCs and O–B2-type stars randomly
into two groups, and calculate the standard kinematic distances for one group
and the revised kinematic distances for the other group. After 100 re-samples
and repetitions of the experiment, the average median value of the offset
between standard kinematic distance and the parallax distance is 2.2 kpc, and
average dispersion is 1.6 kpc. The corresponding values for the revised
kinematic distance and parallax distance are 0.2 and 1.9 kpc, respectively.
The typical relative uncertainty of the kinematic distance for these young
objects is about 0.3. Furthermore, for the offsets in correction that are
listed in Table 5, we check them using the masers, since their distribution
are consistent with those of the YOCs and O–B2-type stars in the $l$–$v$
diagram (Figure 8). As shown in Figure 10(b), compared with the standard
kinematic distances, the revised kinematic distances of the masers are in
agreement with their parallax distances, and the median offset is reduced from
2.0 to 0.3 kpc. The offsets, which are obtained from a fraction of YOCs and
O–B2-type stars, are successfully applied not only all YOCs and O–B2-type
stars, but also masers. It suggests that the method adopted here is effective.

The above results indicates that, within the range of 6 kpc from the Sun, if a
source satisfies the conditions listed in Table 5, one can add the
corresponding offset to its observed radial velocity to determine a more
reliable kinematic distance than the standard kinematic distance. However,
using the present-day data, we can make kinematic distance correction in only
six regions in the solar neighborhood. With an increase of available data in
the future, we are expected to modify the kinematic distance model in all
regions near the Sun, and even in the distant regions of the Milky Way.

##  5 Summary

We investigate the kinematics of YOCs and O–B2-type stars in the local spiral
structure. The YOCs and O–B2-type stars in the different spiral arms show
different kinematic features. The young populations located in the outer
spiral arms may exhibit larger peculiar motions than those in the inner spiral
arms. In the Perseus and Local Arms, the YOCs and O–B2-type stars perhaps move
towards the GC and rotate slower than the rotation of the Milky Way. In
contrast, the young populations in the Carina Arm may move toward the Galactic
anticenter and rotate slightly faster of Galactic rotation. Using the
astrometric parameters of the young stellar populations, we also perform
kinematic distance corrections for several regions in the solar neighborhood.
At present, the corrections to the kinematic distance estimates are only valid
up to 6 kpc from the Sun, yet it provides a useful reference for future
kinematic distance estimates.

Figure 8: (a): $l$–$v$ diagram. Colored dots are the sources of masers (orange
dots), YOCs (red triangles), and O–B2-type stars (cyan diamonds). (b):
variations in the average $V_{\rm LSR}$ per 20∘ for the masers (orange), YOCs
(red), O–B2-type stars (cyan) along Galactic longitude. Each interval must
have more than five sources. Figure 9: $V_{\rm LSR}$ values of the YOCs
(triangles) and O–B2-type stars (diamonds) being assigned to the corresponding
spiral arms (lines) as a function of Galactic longitude. The magenta, blue,
and green lines represent the Perseus Arm, the Local Arm, and the Carina Arm,
respectively, whose $V_{\rm LSR}^{\rm mod}$ values are calculated by using
Equation (2).

(a) O–B2-type stars and YOCs (b) Masers

Figure 10: Comparison between the parallax and kinematic distances. The young stellar population is shown in the left panel, and the masers are shown in the right panel. $D_{\pi}$ is the distance converted from the measured parallax. $D_{k}$ is the kinematic distance. The standard kinematic distance is obtained without considering the peculiar motion. The revised kinematic distance is obtained by adding the offsets to radial velocity to counteract the peculiar motion. Table 5: Offsets Used for the Kinematic Distance Estimation  $l~{}(^{\circ})$ |  condition |  offset (${\rm km~{}s}^{-1}$)
---|---|---
60–90 |  $V_{\rm LSR}<-10$ ${\rm km~{}s}^{-1}$ |  $~{}~{}~{}28$
90–130 |  $V_{\rm LSR}<-45$ ${\rm km~{}s}^{-1}$ |  $~{}~{}~{}20$
130–150 |  $V_{\rm LSR}<-50$ ${\rm km~{}s}^{-1}$ |  $~{}~{}~{}27$
230–250 |  $V_{\rm LSR}>~{}~{}50$ ${\rm km~{}s}^{-1}$ |  $-10$
250–280 |  $V_{\rm LSR}<~{}~{}~{}~{}0$ ${\rm km~{}s}^{-1}$ |  $~{}~{}~{}30$
290–330 |  $V_{\rm LSR}<-40$ ${\rm km~{}s}^{-1}$ |  $~{}~{}~{}25$

We appreciate the anonymous referee for the instructive comments that helped
us to improve the paper. This work was funded by the NSFC Grands 11933011,
National SKA Program of China (Grant No. 2022SKA0120103) and the Key
Laboratory for Radio Astronomy. L.Y.J. thanks the support of the NSFC grant
No. 12203104, the Natural Science Foundation of Jiangsu Province (grant No.
BK20210999). 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.

## References

  * Bian et al. (2022) Bian, S. B., Xu, Y., Li, J. J., et al. 2022, AJ, 163, 54, doi: 10.3847/1538-3881/ac3d90
  * Binney & Tremaine (2008) Binney, J., & Tremaine, S. 2008, Galactic Dynamics: Second Edition (Princeton University Press)
  * Brunthaler et al. (2011) Brunthaler, A., Reid, M. J., Menten, K. M., et al. 2011, AN, 332, 461, doi: 10.1002/asna.201111560
  * Castro-Ginard et al. (2022) Castro-Ginard, A., Jordi, C., Luri, X., et al. 2022, A&A, 661, A118, doi: 10.1051/0004-6361/202142568
  * Chen et al. (2013) Chen, Y., Zhou, P., & Chu, Y.-H. 2013, ApJ, 769, L16, doi: 10.1088/2041-8205/769/1/L16
  * Dame et al. (2001) Dame, T. M., Hartmann, D., & Thaddeus, P. 2001, ApJ, 547, 792, doi: 10.1086/318388
  * Dias et al. (2019) Dias, W. S., Monteiro, H., Lépine, J. R. D., & Barros, D. A. 2019, MNRAS, 486, 5726, doi: 10.1093/mnras/stz1196
  * Fabricius et al. (2021) Fabricius, C., Luri, X., Arenou, F., et al. 2021, A&A, 649, A5, doi: 10.1051/0004-6361/202039834
  * Gaia Collaboration et al. (2018a) Gaia Collaboration, Brown, A. G. A., Vallenari, A., et al. 2018a, A&A, 616, A1, doi: 10.1051/0004-6361/201833051
  * Gaia Collaboration et al. (2018b) Gaia Collaboration, Katz, D., Antoja, T., et al. 2018b, A&A, 616, A11, doi: 10.1051/0004-6361/201832865
  * Gaia Collaboration et al. (2021) Gaia Collaboration, Brown, A. G. A., Vallenari, A., et al. 2021, A&A, 649, A1, doi: 10.1051/0004-6361/202039657
  * Gaia Collaboration et al. (2022a) Gaia Collaboration, Drimmel, R., Romero-Gomez, M., et al. 2022a, arXiv e-prints, arXiv:2206.06207, doi: 10.48550/arXiv.2206.06207
  * Gaia Collaboration et al. (2022b) Gaia Collaboration, Vallenari, A., Brown, A. G. A., et al. 2022b, arXiv e-prints, arXiv:2208.00211, doi: 10.48550/arXiv.2208.00211
  * Georgelin & Georgelin (1976) Georgelin, Y. M., & Georgelin, Y. P. 1976, A&A, 49, 57
  * Hao et al. (2022a) Hao, C. J., Xu, Y., Wu, Z. Y., et al. 2022a, A&A, 668, A13, doi: 10.1051/0004-6361/202244570
  * Hao et al. (2022b) —. 2022b, A&A, 660, A4, doi: 10.1051/0004-6361/202243091
  * Hao et al. (2021) Hao, C. J., Xu, Y., Hou, L. G., et al. 2021, A&A, 652, A102, doi: 10.1051/0004-6361/202140608
  * Hou & Han (2014) Hou, L. G., & Han, J. L. 2014, A&A, 569, A125, doi: 10.1051/0004-6361/201424039
  * Immer et al. (2019) Immer, K., Li, J., Quiroga-Nuñez, L. H., et al. 2019, A&A, 632, A123, doi: 10.1051/0004-6361/201834208
  * Immer & Rygl (2022) Immer, K., & Rygl, K. L. J. 2022, Universe, 8, 390, doi: 10.3390/universe8080390
  * Lindegren et al. (2021) Lindegren, L., Bastian, U., Biermann, M., et al. 2021, A&A, 649, A4, doi: 10.1051/0004-6361/202039653
  * Persic et al. (1996) Persic, M., Salucci, P., & Stel, F. 1996, MNRAS, 281, 27, doi: 10.1093/mnras/278.1.27
  * Qin et al. (2023) Qin, S., Zhong, J., Tang, T., & Chen, L. 2023, ApJS, 265, 12, doi: 10.3847/1538-4365/acadd6
  * Reid et al. (2016) Reid, M. J., Dame, T. M., Menten, K. M., & Brunthaler, A. 2016, ApJ, 823, 77, doi: 10.3847/0004-637X/823/2/77
  * Reid et al. (2009) Reid, M. J., Menten, K. M., Zheng, X. W., et al. 2009, ApJ, 700, 137, doi: 10.1088/0004-637X/700/1/137
  * Reid et al. (2014) Reid, M. J., Menten, K. M., Brunthaler, A., et al. 2014, ApJ, 783, 130, doi: 10.1088/0004-637X/783/2/130
  * Reid et al. (2019) —. 2019, ApJ, 885, 131, doi: 10.3847/1538-4357/ab4a11
  * Sakai et al. (2019) Sakai, N., Reid, M. J., Menten, K. M., Brunthaler, A., & Dame, T. M. 2019, ApJ, 876, 30, doi: 10.3847/1538-4357/ab12e0
  * Soubiran et al. (2018) Soubiran, C., Cantat-Gaudin, T., Romero-Gómez, M., et al. 2018, A&A, 619, A155, doi: 10.1051/0004-6361/201834020
  * VERA Collaboration et al. (2020) VERA Collaboration, Hirota, T., Nagayama, T., et al. 2020, PASJ, 72, 50, doi: 10.1093/pasj/psaa018
  * Xu et al. (2023) Xu, Y., Hao, C. J., Liu, D. J., et al. 2023, ApJ, 947, 54, doi: 10.3847/1538-4357/acc45c
  * Xu et al. (2021) Xu, Y., Hou, L. G., Bian, S. B., et al. 2021, A&A, 645, L8, doi: 10.1051/0004-6361/202040103
  * Xu et al. (2006) Xu, Y., Reid, M. J., Zheng, X. W., & Menten, K. M. 2006, Science, 311, 54, doi: 10.1126/science.1120914
  * Xu et al. (2013) Xu, Y., Li, J. J., Reid, M. J., et al. 2013, ApJ, 769, 15, doi: 10.1088/0004-637X/769/1/15
  * Xu et al. (2018) Xu, Y., Bian, S. B., Reid, M. J., et al. 2018, A&A, 616, L15, doi: 10.1051/0004-6361/201833407
  * Zhang et al. (2019) Zhang, B., Reid, M. J., Zhang, L., et al. 2019, AJ, 157, 200, doi: 10.3847/1538-3881/ab141d