00footnotetext: Received

# Classify the Higgs decays with the PFN and ParticleNet at electron-positron
colliders

Gang Li1,1<EMAIL_ADDRESS>Libo Liao2,2<EMAIL_ADDRESS>Xinchou
Lou1 Peixun Shen1 Weimin Song3
Shudong Wang1,4,3<EMAIL_ADDRESS>and Zhaoling Zhang3 1 Institute of High
Energy Physics, Chinese Academy of Sciences, 19B Yuquan Road, Shijingshan
District, Beijing, China
2 Wuzhou University, 82 Fumin Third Road, Wanxiu District, Wuzhou, China,
3 College of Physics, Jilin University, 2699 Qianjin Street, Changchun, China,
4 School of Physical Sciences, University of Chinese Academy of Sciences,
Beijing, China

###### Abstract

Various Higgs factories are proposed to study the Higgs boson precisely and
systematically in a model-independent way. In this study, the Particle Flow
Network and ParticleNet techniques are used to classify the Higgs decays into
multi-categories and the ultimate goal is to realize an “end-to-end” analysis.
A Monte Carlo simulation study is performed to demonstrate the feasibility,
and the performance looks rather promising. This result could be the basis of
a “one-shop” analysis to measure all the branching fractions of the Higgs
decays simultaneously.

###### keywords:

the Higgs Boson, event classification, Particle Flow Network, ParticleNet

##  1 Introduction

The historic observation of the Higgs boson in 2012 at the Large Hadron
Collider (LHC) [1, 2] declared the discovery of the last missing piece of the
most fundamental building blocks in the Standard Model (SM). Although the SM
has been remarkably successful in describing experimental phenomena, a
precision Higgs physics program would be critically important given that the
SM does not predict the parameters in the Higgs potential, nor does it involve
particle candidates for dark matter. The precision determination of the Higgs
couplings to the SM particles, gauge bosons and leptons/quarks, are the agents
probing the Higgs mechanism for generating the masses [3]. In particular,
potential observable deviations of the Higgs couplings from the SM
expectations would indicate new physics [4]. Therefore, the Higgs discovery
marks the beginning of a new era of both theoretical and experimental
explorations. Various $e^{+}e^{-}$ colliders as Higgs factories were proposed
by the high energy physics community, such as ILC [5], CLIC [6], FCC-ee [7],
and CEPC [8, 9].

The most important advantages of a Higgs factory are that the center of mass
(CM) energy is precisely defined and that they could perform absolute
measurements to the Higgs boson. Neglecting the $Z$ fusion production, in an
$e^{+}e^{-}\to ZH$ event, where the $Z$ decays to a pair of visible fermions (
$Z\rightarrow
e^{+}e^{-},~{}\mu^{+}\mu^{-},~{}\tau^{+}\tau^{-},~{}\mbox{or}~{}q\bar{q}$ ),
the Higgs boson can be identified with the kinematics of these fermion pairs
independent of its decays. The production cross-section and most of the decay
branching fractions of the Higgs could be measured model-independently by the
counting method. For example, the CEPC can measure the cross-section of
$e^{+}e^{-}\rightarrow ZH$, $\sigma(ZH)$, at 240 GeV, to a precision of 0.5%
and the branching fractions to a few percent, respectively, by combining the
four decay modes of the $Z$ boson [9, 10].

The physics goal of a Higgs factory must be accomplished by optimizing the
detector design and making use of the latest developments in data science.
Recently, various Machine Learning (ML) techniques have already shown very
promising performances in the data analysis of high energy physics because of
their strong inductive biases[11], in particular for jet studies. For
instance, jets are treated as images [12, 13, 14, 15, 16, 17] or as sequences
[18, 19, 20, 21], trees [22, 23], graphs [24], or sets [25, 26] of particles,
and ML techniques, most notably deep neural networks (DNNs), are used to build
new jet tagging algorithms automatically from (labeled) simulated samples and
even data [27, 28, 29, 30]. While the above ML techniques are used at jet-
level for case studies, they naturally can apply for the event level in
$e^{+}e^{-}$ collision, which has much simpler topology and pile-up free.

In this article, two ML approaches are used studying the classification
problem of Higgs events. The classification results can serve as the basis of
an “end-to-end” (E2E) analysis, which means that it starts from particle level
information and is efficient and balanced because it can study almost all the
Higgs decays modes simultaneously through state-of-the-art ML techniques. The
approach also is a ”one-shop” analysis to support extracting all Higgs
couplings and taking into account the correlations and commonalities of the
same detector for the experiment.

The rest of this paper is organized as follows. The ML methods used in this
study are introduced in Sec. 2, followed by the implementation of the ML with
a Monte Carlo (MC) simulation in Sec. 3. Finally, a summary is presented.

##  2 Machine Learning methods

Recently, various ML techniques were proposed for jet tagging studies. Among
them, the PFN [25] and ParticleNet [26] got relatively superior performance.

In [25] the authors applied the Deep Sets concept [34] to the jet tagging
problem. They propose two elegant model architectures, named EnergyFlow
Network (EFN) and ParticleFlow Network (PFN), with provable physics
properties, such as infrared and colinear safety. In these two architectures,
the features of each particle are encoded into a latent space of $\Phi$ [34]
and the category, $F$, is extracted from the summed representation in that
latent space. Both $\Phi$ and $F$ are approximated by neural networks. The key
mathematical fact is that a generic function of a set of particles can be
decomposed into an arbitrarily good approximation according to the Deep Set
Theorem [34]. The performance of these models in classification problems is
comparable with other more complicated models. The authors also tried to
interpret and visualize what the model has learned [25].

Motivated by the success of CNNs, the ParticleNet [26] approach based on the
Dynamic Graph Convolutional Neural Network (DGCNN) is proposed for learning on
particle cloud data. The edge convolution (“EdgeConv”) operation, a
convolution-like operation for point clouds, is used instead of the regular
convolution operation. One important feature of the EdgeConv operation is that
it can be easily stacked, just as regular convolutions. Therefore, another
EdgeConv operation can be applied subsequently, which make it possible to
learn features of point clouds hierarchically. Another important feature is
that the proximity of points can be dynamically learned with EdgeConv
operations. The study shows the graph describing the point clouds are
dynamically updated to reflect the changes in the edges, i.e., the neighbors
of each point. Reference [26] shows that this leads to better performance than
keeping the graph static.

As suggested by the authors [25] and according to the performances of EFN and
PFN, we choose the PFN and ParticleNet to classify the Higgs decays. This ML
attempt contains some distinct features in contrast to conventional data
analysis. First, the ML approach is used to classify many physics processes at
the same time. If some tiny decays are neglected, there are about 9 branching
fractions of the Higgs decay to be measured. The number of classes is greater
than 9 when the SM backgrounds are included. In addition, the classification
results could be the basis of an E2E analysis, which means that all the
particle level information, such as four momenta, PID, and impact parameters
of charged particles, is used as input directly, and the network calculates
the scores of each event. In this case, the analysis no longer needs some
dedicated and complicated reconstruction tools, such as lepton/photon
isolation, jet-clustering, and $\tau$ finder, and so on.

##  3 Classify the Higgs decays

In this section, a 9-classification of the all accessible Higgs decay final
states is realized with the PFN method, and the confusion matrices of them are
determined. And as an attempt, a more ambitious 39-classification is tried
with the ParticleNet, and promising and consistent results are achieved.

###  3.1 ML model setup

For the two ML models used to classify the Higgs decays, kinematic information
of energy, polar and azimuthal angles are always given for each particle. We
should note that energies and polar angles are used instead of the transverse
momenta and rapidities in the original study [25, 26], respectively, since the
models are utilized for $e^{+}e^{-}$ collider experiments in this study. All
features are pre-processed: the energies are divided by the scalar $\sum E$
and the polar and azimuthal angles are centered based on those of the event.
The inputs also include the PID and impact parameters of charged particles.

The PFN architecture [25] is designed to parameterize the functions $\Phi$ and
$F$ in a sufficiently general way, several dense neural network layers are
used as universal approximators. For $\Phi$, three dense layers are employed
with 100, 100, and $l$ nodes, respectively, where $l$ is the latent dimension
that takes 256 after comparing the performances of 128 and 256. For $F$, we
use the same configuration as the original paper, which has three dense
layers, each with 100 nodes. Each dense layer uses the ReLU activation
function and He-uniform parameter initialization [38]. A nine-unit layer
(depends on the number of classes) with a SoftMax activation function is the
output layer.

The ParticleNet[26] architecture consists of three EdgeConv blocks, one
aggregation layer, and two fully-connected layers. The first EdgeConv block
uses the spatial coordinates of the particles in the $\theta-\phi$ space to
compute the distances, while the subsequent blocks use the learned feature
vectors as coordinates. The number of nearest neighbors $k$ is 16 for all
three blocks, and the number of channels $C$ for each EdgeConv block is (64,
64, 64), (128, 128, 128), and (256, 256, 256), respectively. After the
EdgeConv blocks, a channel-wise global average pooling operation is applied to
aggregate the learned features over all particles in the cloud. This is
followed by a fully connected layer with 256 units and the ReLU activation. A
dropout layer with a drop probability of 0.1 is included to prevent
overfitting. A fully connected layer with 39 units, followed by a softmax
function, is used to generate the output for the 39 classification task.

###  3.2 Data samples

In this study, there are 4 production modes of the Higgs boson at 240 GeV to
be analyzed, i.e., $e^{+}e^{-}\to e^{+}e^{-}H$, $\mu^{+}\mu^{-}H$,
$\tau^{+}\tau^{-}H$, and $q\bar{q}H$. And in each production mode, the same 9
decay modes are measured, which are $H\to c\bar{c}$, $b\bar{b}$,
$\mu^{+}\mu^{-}$, $\tau^{+}\tau^{-}$, $gg$, $\gamma\gamma$, $ZZ$,
$W^{+}W^{-}$, and $\gamma Z$, respectively. So there are 36 processes in
total. For each process, 400,000 events are generated with the WHIZARD 1.9.5
[35] and fed to Pythia6 [36] for hadronization. It should be noticed that the
sequential decays of $W$ and $Z$ are not dealt specifically to avoid the
complication, though it can enhance the classification performance if more
decay knowledge is used.

All the generated samples are simulated by a simplified way to model detector
responses. In detail, all particles are simulated according to the performance
of the baseline detector in the CEPC CDR [9]. The momentum resolution of
charged tracks is $\frac{\sigma(p_{t})}{p_{t}}=2\times
10^{-5}\oplus\frac{0.001}{p\sin^{3/2}\theta}[\mathrm{GeV}^{-1}]$, the energy
resolution of photons is
$\frac{\sigma(E)}{E}=0.01\oplus\frac{0.20}{\sqrt{E/\mathrm{(GeV)}}}$, and this
of neural hadrons is
$\frac{\sigma(E)}{E}=0.03\oplus\frac{0.50}{\sqrt{E/\mathrm{(GeV)}}}$, which
means all the particles are reconstructed with a perfect particle flow
detector. In case of impact parameters and particle identification, they are
taken directly from the truth of generation. While the simulation is a bit too
ideal, it is sufficient for a feasibility study.

###  3.3 9-category classification: training and evaluation

During model training, the common properties of the neural network include
categorical cross-entropy loss function, the Adam optimization algorithm [37],
batch size of 1000, and a learning rate of 0.001. 400,000 events are used for
each production mode, and the total number of events of 9 decays is 3,600,000.
The full data set is split into training, validation, and test samples
according to fractions of 8:1:1. The monitoring of loss and accuracy on
training and validation samples shows that the models converge well and there
is no obvious over-training after the models are trained for 100 epochs, see
Figure 1 as an example.

Figure 1:  Accuracy and loss versus the number of epochs of $e^{+}e^{-}\to
e^{+}e^{-}H$ process during the training.

The outputs of the classifier, which are from a nine-unit layer with a SoftMax
activation function, are visualized in various ways. The SoftMax is essential
because it helps to produce scores comprising 9 probabilities proportional to
the exponential of the input information for each event, which is input for a
cut-based data analysis. Figure 2 presents the 9 scores for each category.
Taking the bottom left panel as an example, these events are of $H\to ZZ$, and
the curves in different colors represent the probability distributions
assuming them to be $H\to ZZ$. The blue curve peaks when the score approaches
1, which means the classifier can identify $H\to ZZ$ signals. There are two
small peaks in the blue and brown curves around 0.8, which shows that $H\to
ZZ$ and $H\to\gamma Z$ can contaminate each other because of the similarity of
their cascade decays. From Figure 2, it can be seen that high-dimensional data
is difficult to visualize intuitively. A better way is that data in lower
dimensions are plotted to show the inherent structures. To aid visualization
of the structure of 9 outputs, the $t$-SNE[39] method is used. Figure 3 shows
the distribution of the two largest components after the dimensionality
reduction, where 1-9 represents the 9 decay modes of the Higgs boson from
$c\bar{c}$ to $\gamma Z$ in the same order as the above, respectively. The
patterns in Figure 3 are consistent with those in Figure 2 but much clearer.
It can be seen that $\mu^{+}\mu^{-}$(3), $\gamma\gamma$(6),
$\tau^{+}\tau^{-}$(4) and $\gamma Z$(9) modes are almost isolated clusters and
background free. The clusters of the others can also be seen and the overlaps
are also significant.

Figure 2:  The distributions of 9 outputs for each true category, taking
$e^{+}e^{-}H$ as an example. Each score is calculated by assuming that the
event belongs to that category.  Figure 3:  The classification performance on
the test set is visualized with $t$-SNE , where the two largest components are
used, taking 10,000 events of $e^{+}e^{-}\to e^{+}e^{-}H$ for illustration.

Some standard quantities can measure the performances of classifiers. For
instance, accuracy (ACC) measures the fraction of correctly classified
observations, ROC curve (Receiver Operating Characteristic curve) visualizes
the True Positive Rate (TPR) versus the False Positive Rate (FPR), and AUC
(Area Under the Curve) is the area under the ROC curve. If we have a better
classification for each threshold value, the area grows, and a perfect
classification leads to an AUC of 1.0. The average ACC and AUC for all 36
processes in the 4 tagging modes are summarized in Table 1. Several
conclusions can be drawn from the table. First, the average accuracy of each
tagging mode reaches about 87%, which is good and adequate for further
analysis. The decays of $H\to\mu^{+}\mu^{-}$, $\tau^{+}\tau^{-}$, and
$\gamma\gamma$ have the best accuracy and largest AUCs as expected. Last but
not least, the accuracy of $H\to ZZ$ or $W^{+}W^{-}$ is not good as the
others, which leaves room for further improvement.

Finally, the confusion matrices are used to evaluate the performance of the ML
model and to be used as an important gradient for the further data analysis.
Confusion matrices are calculated by comparing the prediction of the model and
the true labels. Figure 4 shows the confusion matrices of the four
classifiers, respectively. In terms of the confusion matrix, the accuracy
appears as the diagonal elements of the corresponding confusion matrices, and
the off-diagonal ones represent misclassification rates. So confusion matrices
contain the complete information of both the correct and incorrect
classifications, which could help to unfold the generated numbers of signals,
$N_{i}$.

Decay Mode | $e^{+}e^{-}H$ | $\mu^{+}\mu^{-}H$ | $\tau^{+}\tau^{-}H$ | $q\bar{q}H$
---|---|---|---|---
ACC | AUC | ACC | AUC | ACC | AUC | ACC | AUC
$H\to c\bar{c}$ | 0.880 | 0.991 | 0.882 | 0.991 | 0.857 | 0.987 | 0.755 | 0.966
$H\to b\bar{b}$ | 0.908 | 0.994 | 0.893 | 0.994 | 0.877 | 0.991 | 0.733 | 0.972
$H\to\mu^{+}\mu^{-}$ | 0.997 | 1.000 | 0.986 | 1.000 | 0.981 | 1.000 | 0.983 | 1.000
$H\to\tau^{+}\tau^{-}$ | 0.993 | 0.999 | 0.985 | 0.999 | 0.985 | 0.999 | 0.982 | 0.999
$H\to gg$ | 0.810 | 0.985 | 0.830 | 0.986 | 0.816 | 0.982 | 0.736 | 0.954
$H\to\gamma\gamma$ | 0.997 | 1.000 | 0.999 | 1.000 | 1.000 | 1.000 | 0.997 | 1.000
$H\to ZZ$ | 0.650 | 0.958 | 0.667 | 0.960 | 0.585 | 0.947 | 0.535 | 0.926
$H\to W^{+}W^{-}$ | 0.806 | 0.981 | 0.801 | 0.981 | 0.771 | 0.974 | 0.632 | 0.952
$H\to\gamma Z$ | 0.921 | 0.996 | 0.936 | 0.996 | 0.910 | 0.993 | 0.896 | 0.993
Table 1:  Accuracy (left) and AUCs (right) of four classifiers.

Figure 4:  Confusion matrices of (a): $e^{+}e^{-}\to e^{+}e^{-}H$, (b):
$e^{+}e^{-}\to\mu^{+}\mu^{-}H$, (c): $e^{+}e^{-}\to\tau^{+}\tau^{-}H$, and
(d): $e^{+}e^{-}\to q\bar{q}H$, respectively.

###  3.4 An attempt of 39-classification

The above study shows that the multi-classification is very promising in data
analysis, so here a more ambitious case of 39-classification will be tried.
For the signal processes, considering that $Z$ decay to 4 categories (neglect
neutrino decay and $W$ fusion processes up to now), i.e, $e^{+}e^{-}$,
$\mu^{+}\mu^{-}$, $\tau^{+}\tau^{-}$, and $q\bar{q}$, and that the Higgs to
the same 9 decay modes as the above, so there are 36 signals. For a realistic
analysis, the backgrounds must be taken into account, especially the
irreducible ones. In the analysis of $e^{+}e^{-}\to ZH$ study, the irreducible
backgrounds mainly come from the SM process of $e^{+}e^{-}\to ZZ$. The
background can be categorized into 3 classes depending on the decays of $Z$
bosons, i.e, pure leptonic ($ZZ_{l}$), semi-leptonic($ZZ_{sl}$), and hadronic
($ZZ_{h}$) decays. In conclusion it is a 39-classification problem.

Same data sets of the signal and extra 3 background processes are pre-
processed with the same procedure, which has 39$\times$400 000= 15 600 000
events, which is very challenging because of memory usage. So we switch to
anther deep learning framework, ParticleNet/Weaver[26, 40], which has more
flexible memory strategy.

The confusion matrix of the 39-classification is presented in Fig. 5, which
shows very good separation power among all 39 processes. For the signal, four
blocks of the $e^{+}e^{-}H$, $\mu^{+}\mu^{-}H$, $\tau^{+}\tau^{-}H$, and
$q\bar{q}H$ processes can be seen clearly, which demonstrate similar patterns
as the corresponding ones in Fig. 4. In each sub-matrix, the
$H\to\gamma\gamma$, $\mu^{+}\mu^{-}$, and $\tau^{+}\tau^{-}$ decays achieve
the best performances. Among the four “blocks”, the mis-classification rates
are rather small. And the $H\to ZZ$ decay doesn’t achieve sufficiently good
performance as the other decays, which is also consistent with those of the
9-classification. For the irreducible backgrounds of $e^{+}e^{\to}ZZ$, all of
three processes are labeled correctly with very high accuracy, greater than
90%, which indicates that the kinematics of different events can be learnt to
discriminate the irreducible backgrounds by the ParticleNet.

In this 39-classification, all 9 Higgs decays in 4 tagging modes with the
irreducible backgrounds together can be classified with rather good accuracy.
It is different with only one single tagging mode, which indicates the the
Higgs decays can be determined with a combined way and much more information.
So better precision could be expected naturally.

Figure 5:  Confusion matrix of the 39-classification

##  4 Summary and discussion

In this paper, we presented a study of the classification of the Higgs decays
with the state-of-the-art ML approaches at electron-positron colliders. It
deploys the ML techniques and try to classify both the signal and background
events with only particle level information and obtain the confusion matrices,
which can be used in the further data analysis. This approach is the basis of
an efficient and balanced “one-shop” analysis, which is possible to measure
all Higgs couplings using all detector information and taking all the
commonalities and correlations into account. For the analyses of tens or
hundreds of channels, they can be repeated using this technique in a few days
if all data samples are ready. In contrast, the time could be considerably
longer using the conventional analysis method.

To be frank, this work is only a feasibility study. There are various
possibilities to improve and further validate it. One is to enhance the
performance by taking the sequential decays of $W$ and $Z$ bosons into
account, which can adopt more information and enhance the classification
performance. Another endeavor with more physical significance is incorporating
some physics processes beyond the SM in the analysis, such as the invisible
and semi-invisible decays of the Higgs boson, which can enhance the
sensitivity of an experiment to new physics. It is also very constructive to
take the full SM backgrounds and main systematic uncertainties into account.

## Acknowledge

The authors present special thanks to Yunxuan Song, Congqiao Li, Dr. Yu Bai,
and Dr. Huilin Qu for useful discussion and advice. The authors thank the IHEP
Computing Center for its firm support. This work is supported by the National
Natural Science Foundation of China (NSFC) under Contracts Nos. 12075271 and
12047569.

## References

  * [1] ATLAS Collaboration (G. Aad et al.), _Observation of a new particle in the search for the Standard Model Higgs boson with the ATLAS detector at the LHC_ , _Phys. Lett. B_ 716 1-29 (2012).
  * [2] CMS Collaboration (S. Chatrchyan et al.), _Observation of a new boson at a mass of 125 GeV with the CMS experiment at the LHC_ , _Phys. Lett. B_ 716 30-61 (2012).
  * [3] R. Lafaye, T. Plehn, M. Rauch, D. Zerwas, and M. Duhrssen, _Measuring the Higgs Sector_ , _JHEP_ 08, 009 (2009), doi:10.1088/1126-6708/2009/08/009, arXiv:0904.3866 (2009).
  * [4] C. Englert, A. Freitas, M. M. Mühlleitner, T. Plehn, M. Rauch, M. Spira, and K. Walz, _Precision Measurements of Higgs Couplings: Implications for New Physics Scales_ , _J. Phys. G_ 41, 113001, doi:10.1088/0954-3899/41/11/113001, arXiv:1403.7191 (2014).
  * [5] H. Baer, T. Barklow, K. Fujii, Y. Gao, A. Hoang, S. Kanemura, J. List, H. E. Logan, A. Nomerotski and M. Perelstein, et al., _The International Linear Collider Technical Design Report - Volume 2: Physics_ , arXiv:1306.6352 (2013).
  * [6] M. Aicheler, P. Burrows, M. Draper, T. Garvey, P. Lebrun, K. Peach, N. Phinney, H. Schmickler, D. Schulte and N. Toge, _A Multi-TeV Linear Collider Based on CLIC Technology: CLIC Conceptual Design Report_ , doi:10.5170/CERN-2012-007 (2012).
  * [7] A. Abada et al., _FCC-ee: The Lepton Collider: Future Circular Collider Conceptual Design Report Volume 2_ , _Eur. Phys. J. ST_ 228, no.2, 261-623 doi:10.1140/epjst/e2019-900045-4 (2019).
  * [8] CEPC Study Group, _CEPC Conceptual Design Report: Volume 1 - Accelerator_ , arXiv:1809.00285 (2018).
  * [9] CEPC Study Group, _CEPC Conceptual Design Report: Volume 2 - Physics & Detector_, arXiv:1811.10545 (2018).
  * [10] F. An, et al., _Precision Higgs physics at the CEPC_ , _Chin. Phys. C_ 43, no.4, 043002, doi:10.1088/1674-1137/43/4/043002, arXiv:1810.09037 (2019).
  * [11] M. D. Schwartz, _Modern Machine Learning and Particle Physics_ , arXiv:2103.12226 (2021).
  * [12] A. J. Larkoski, I. Moult, and B. Nachman, _Jet Sub-structure at the Large Hadron Collider: A Review of Recent Advances in Theory and Machine Learning_ , _Phys. Rept._ 841, 1, arXiv:1709.04464 (2020).
  * [13] R. Kogler et al., _Jet Substructure at the Large Hadron Collider: Experimental Review_ , _Rev. Mod. Phys_ 91, 045003, arXiv:1803.06991 (2019).
  * [14] J. Cogan, M. Kagan, E. Strauss, and A. Schwarztman, _Jet-Images: Computer Vision Inspired Techniques for Jet Tagging_ , _JHEP_ 02, 118, arXiv:1407.5675 (2015).
  * [15] L. G. Almeida, M. Backović, M. Cliche, S. J. Lee, and M. Perelstein, _Playing Tag with ANN: Boosted Top Identification with Pattern Recognition_ , _JHEP_ 07, 086, arXiv:1501.05968 (2015).
  * [16] L. de Oliveira, M. Kagan, L. Mackey, B. Nachman, and A. Schwartzman, _Jet-images – deep learning edition_ , _JHEP_ 07, 069, arXiv:1511.05190 (2016).
  * [17] P. Baldi, K. Bauer, C. Eng, P. Sadowski, and D. Whiteson, _Jet Substructure Classification in High-Energy Physics with Deep Neural Networks_ , _Phys. Rev. D_ 93, 094034, arXiv:1603.09349 (2016).
  * [18] D. Guest, J. Collado, P. Baldi, S.-C. Hsu, G. Urban, and D. Whiteson, _Jet Flavor Classification in High-Energy Physics with Deep Neural Networks_ , _Phys. Rev. D_ 94, 112002, arXiv:1607.08633 (2016).
  * [19] J. Pearkes, W. Fedorko, A. Lister, and C. Gay, _Jet Constituents for Deep Neural Network Based Top Quark Tagging_ , arXiv:1704.02124 (2017).
  * [20] S. Egan, W. Fedorko, A. Lister, J. Pearkes, and C. Gay, _Long Short-Term Memory (LSTM) networks with jet constituents for boosted top tagging at the LHC_ , arXiv:1711.09059 (2017).
  * [21] K. Fraser and M. D. Schwartz, _Jet Charge and Machine Learning_ , _JHEP_ 10, 093, arXiv:1803.08066 (2018).
  * [22] G. Louppe, K. Cho, C. Becot, and K. Cranmer, _QCD-Aware Recursive Neural Networks for Jet Physics_ , _JHEP_ 01, 057, arXiv:1702.00748 (2019).
  * [23] T. Cheng, _Recursive Neural Networks in Quark/Gluon Tagging_ , _Comput. Softw. Big Sci._ 2, 3, arXiv:1711.02633 (2018).
  * [24] I. Henrion, J. Brehmer, J. Bruna, K. Cho, K. Cranmer, G. Louppe, and G. Rochette, _Neural Message Passing for Jet Physics_ , _Deep Learning for Physical Sciences Workshop at the 31st Conference on Neural Information Processing Systems (NIPS)_ (2017).
  * [25] Patrick T. Komiske, Eric M. Metodiev, and Jesse Thaler, _Energy flow networks: deep sets for particle jets_ , _Journal of High Energy Physics_ , 1, Jan (2019).
  * [26] H. Qu and L. Gouskos, _ParticleNet: Jet Tagging via Particle Clouds_ , _Phys. Rev. D_ 101, no.5, 056019, doi:10.1103/PhysRevD.101.056019, arXiv:1902.08570 (2020).
  * [27] E. M. Metodiev, B. Nachman, and J. Thaler, _Classification without labels: Learning from mixed samples in high energy physics_ , _JHEP_ 10, 174, arXiv:1708.02949 (2017).
  * [28] P. T. Komiske, E. M. Metodiev, B. Nachman, and M. D. Schwartz, _Learning to classify from impure samples with high-dimensional data_ , _Phys. Rev. D_ 98, 011502(R), arXiv:1801.10158 (2018).
  * [29] A. Andreassen, I. Feige, C. Frye, and M. D. Schwartz, _JUNIPR: a Framework for Unsupervised Machine Learning in Particle Physics_ , _Eur. Phys. J. C_ 79, 102, arXiv:1804.09720 (2019).
  * [30] P. T. Komiske, E. M. Metodiev, and J. Thaler, _An operational definition of quark and gluon jets_ , _JHEP_ 11, 059, arXiv:1809.01140 (2018).
  * [31] ALEPH Collaboration (S. Schael et al.), _Branching ratios and spectral functions of tau decays: Final ALEPH measurements and physics implications_ , _Phys. Rept._ 421, 191-284 doi:10.1016/j.physrep.2005.06.007 arXiv:0506072 (2005).
  * [32] BESIII Collaboration (M. Ablikim et al.), _Precision measurement of the $D^{*0}$ decay branching fractions_, _Phys. Rev. D_ 91 no.3, 031101, doi:10.1103/PhysRevD.91.031101, arXiv:1412.4566 (2015).
  * [33] Particle Data Group (P.A. Zyla et al.), _The Review of Particle Physics (2020)_ , _Prog. Theor. Exp. Phys._ 2020, 083 C01 (2020).
  * [34] Manzil Zaheer, Satwik Kottur, Siamak Ravanbakhsh, Barnabas Poczos, Ruslan Salakhutdinov, Alexander Smola, _Deep Sets_ , arXiv:1703.06114 (2017).
  * [35] Kilian, W., Ohl, T. & Reuter, J., _WHIZARD - simulating multi-particle processes at LHC and ILC_ , _Eur. Phys. J. C_ 71, 1742, doi:10.1140/epjc/s10052-011-1742-y (2011).
  * [36] T. Sjostrand, S. Mrenna and P. Z. Skands, _PYTHIA 6.4 Physics and Manual_ , _JHEP_ 05, 026, doi:10.1088/1126-6708/2006/05/026, arXiv:0603175 (2006).
  * [37] Diederik P. Kingma and Jimmy Ba, _Adam: A Method for Stochastic Optimization_ , _ICLR_ ArXiv:1412.6980 (2015).
  * [38] K. He, X. Zhang, S. Ren and J. Sun, _Delving deep into rectifiers: surpassing human-level performance on ImageNet classification_ , in 2015 IEEE International Conference on Computer Vision (ICCV), IEEE, Santiago, Chile, pg. 1026.c (2015).
  * [39] Laurens van der Maaten and Geoffrey Hinton, _Visualizing Data using $t$-SNE_, _Journal of Machine Learning Research_ 9 2579 - 2605 (2008).
  * [40] H. Qu, Weaver, _a streamlined yet flexible machine learning R &D framework for HEP_, https://github.com/hqucms/weaver.


# The Flare and Warp of the Young Stellar Disk traced with LAMOST DR5 OB-type
stars

Yang Yu(于扬)  Equal first authors  Department of Physics, Hebei Normal
University, Shijiazhuang 050024, China  Hai-Feng Wang(王海峰) Equal first authors
GEPI, Observatoire de Paris, Université PSL, CNRS, Place Jules Janssen 92195,
Meudon, France  South-Western Institute for Astronomy Research, Yunnan
University, Kunming, 650500, China  National Astronomical Observatories,
Chinese Academy of Sciences, Beijing 100101, China  Wen-Yuan Cui(崔文元)
Department of Physics, Hebei Normal University, Shijiazhuang 050024, China
Lin-Lin Li(李林林)  Department of Physics, Hebei Normal University, Shijiazhuang
050024, China  Chao Liu(刘超) National Astronomical Observatories, Chinese
Academy of Sciences, Beijing 100101, China  Bo Zhang(章博) National Astronomical
Observatories, Chinese Academy of Sciences, Beijing 100101, China  Department
of Astronomy, Beijing Normal University, Beijing, 100875, China  Hao Tian(田浩)
National Astronomical Observatories, Chinese Academy of Sciences, Beijing
100101, China  Zhen-Yan Huo(霍振燕)  Department of Physics, Hebei Normal
University, Shijiazhuang 050024, China  Jie Ju(巨洁)  Department of Physics,
Hebei Normal University, Shijiazhuang 050024, China  Zhi-Cun Liu(柳志存)
National Astronomical Observatories, Chinese Academy of Sciences, Beijing
100101, China  Fang Wen(温芳)  Department of Physics, Hebei Normal University,
Shijiazhuang 050024, China  Shuai Feng(冯帅)  Department of Physics, Hebei
Normal University, Shijiazhuang 050024, China Wen-Yuan Cui(崔文元)
<EMAIL_ADDRESS><EMAIL_ADDRESS>
###### Abstract

We present analysis of the spatial density structure for the outer disk from
8$-$14 kpc with the LAMOST DR5 13534 OB-type stars and observe similar flaring
on north and south sides of the disk implying that the flaring structure is
symmetrical about the Galactic plane, for which the scale height at different
Galactocentric distance is from 0.14 to 0.5 kpc. By using the average slope to
characterize the flaring strength we find that the thickness of the OB stellar
disk is similar but flaring is slightly stronger compared to the thin disk as
traced by red giant branch stars, possibly implying that secular evolution is
not the main contributor to the flaring but perturbation scenarios such as
interactions with passing dwarf galaxies should be more possible. When
comparing the scale height of OB stellar disk of the north and south sides
with the gas disk, the former one is slightly thicker than the later one by
$\approx$ 33 and 9 pc, meaning that one could tentatively use young OB-type
stars to trace the gas properties. Meanwhile, we unravel that the radial scale
length of the young OB stellar disk is 1.17 $\pm$ 0.05 kpc, which is shorter
than that of the gas disk, confirming that the gas disk is more extended than
stellar disk. What is more, by considering the mid-plane displacements
($Z_{0}$) in our density model we find that almost all of $Z_{0}$ are within
100 pc with the increasing trend as Galactocentric distance increases.

Milky Way disk (1050); Milky Way evolution (1052); Milky Way formation (1053);
Milky Way Galaxy (1054)

##  1 Introduction

Disk galaxy structure and origins are central problems in galaxy formation and
evolution. The star counts are commonly used method to explore the structure
of the disk galaxy such as our home galaxy, Milky Way (MW), and it also has
been proved as a very effective tool in revealing the structure and
substructure such as flare, warp, truncation, spiral arms, etc, thus then it
will push our understanding of the Galaxy formation and evolution (van der
Kruit & Freeman, 2011; Bland-Hawthorn & Gerhard, 2016). A dual exponential
model is widely used to study the spatial structure of disk: one for thin disk
and the other one for thick disk (Gilmore & Reid, 1983; Jurić et al., 2008;
Bovy et al., 2012; Rix & Bovy et al., 2013; Bovy et al., 2016; Bland-Hawthorn
& Gerhard, 2016). An alternative form, that is, $sech$ function is also often
used in the vertical density profile model (van der Kruit, 1988; van der Kruit
& Freeman, 2011; Wang et al., 2018b).

With the help of Red Giant Branch (RGB) stars, Wang et al. (2018b) revealed
that the Galactic disk radial profile is composed by three sections of
exponential law with the scale length of $2.12\pm 0.26$, $1.18\pm 0.08$, and
$2.72$ kpc respectively at $R<11$, $11\leq R\leq 14$, and $R>14$ kpc.
Meanwhile, they also found the scale length of the thick disk is larger than
that of the thin disk and some density asymmetries might be explained, to some
extent, by shifting either the mid-plane displacements of the thin or the
thick disk. What’s more, they also detected clear flaring features, i.e., the
scale height increases with the distance, from 8 to 19 kpc.

Jurić et al. (2008) found that the scale length of the thin disk is shorter
than that of the thick disk based on Sloan Digital Sky Survey. Similar results
were shown in Chen et al. (2017) by using Photometric data from the Xuyi
Schmidt Telescope Photometric Survey of the Galactic Anti-Centre and the Sloan
Digital Sky Survey. However, Wan, Liu & Deng et al. (2017) suggested that the
younger geometrical thin disk with 4.7 kpc of the exponential scale length is
larger than the older thick disk with 3.4 kpc of the exponential scale length
with two components, which has been argued by Bovy et al. (2012) using mono-
abundance population methods. The geometry definition mentioned in this work
is based on the density profile(Martig et al., 2016).

The scale height of the thin disk is around 220–450 pc and the thick disk is
between 700 and 1200 pc in the solar vicinity which are either from
exponential model or $sech$ function for the vertical component (Bland-
Hawthorn & Gerhard, 2016). As we move forward to the larger Galactic distance,
the value increases due to the flaring effects, in another words, because the
disk flares, the exponential scale height increases with increasing
Galactocentric radius. Recently, with 250,000 OB stars containing all types of
Gaia (Gaia Collaboration: Prusti et al., 2016) and 2MASS (Skrutskie et al.,
2006) photometric catalog, Li et al. (2019) showed that the scale length is
$2.10\pm 0.01$ kpc and scale height is from 132$-$450 pc.

It appears that secular heating of the stellar disk caused by the scattering
of the spiral structures (Sellwood et al., 1984; Carlberg & Sellwood et al.,
1985) or giant molecular clouds (Jenkins & Binney et al., 1990; Jenkins et
al., 1992) might cause the feature that the scale height of the younger
populations is shorter than that of the older ones (Yu & Liu, 2017; Liu et
al., 2012).The well known flare was explained by migration of stars (Solway et
al., 2012), the larger vertical excursions due to the decreased gravitational
pull in the outer disk can lead to the flaring (Bovy et al., 2016). An
alternating explanation was the dynamical heating (Quinn et al., 1993).
However, as first proposed by Minchev et al. (2012), flaring by radial
migration does happen, but the corresponding features are completely wiped out
by mergers shown in Fig. 2 of Minchev et al. (2014). Followed by this work,
Minchev et al. (2015, 2016, 2018) suggested that migration suppresses flaring
in the presence of external perturbations in CDM cosmology (Quinn et al.,
2009). This is in agreement with some other cosmological simulations, e.g.,
Grand et al. (2016); Ma et al. (2017).

As far as we know about the scenarios for the flaring in general, one is the
secular evolution, e.g., mentioned in the works of Minchev et al. (2012,
2014); Narayan & Jog et al. (2006), the other is perturbation, e.g., passing
dwarf galaxies or satellites denoted in Kazantzidiset al. (2008); Villalobos
et al. (2008); Laporte et al. (2018), and so on.

It is almost a consensus that warp exists in the Milky Way, that is describing
the disk bends upwards and downwards in the north and south hemisphere
separately. The kinematical and spatial signals of the stellar warp along with
radius and varies with the azimuth angle are shown in some works (Chen et al.,
2019; Liu et al., 2017c; López-Corredoira et al., 2014). The dynamical signals
of the warp corresponding to the increasing trend of vertical velocity with
vertical angular momentum was displayed in Schönrich et al. (2017); Huang et
al. (2018), and Wang et al. (2020a) also proposed that most likely the main
S-shaped structure of the warp is a long-lived nonsteady feature.

According to our current understanding, it seems that previous studies focused
on the structure of disk tended to use older tracers for the reason that the
number of old stars is quite large, especially for the outer disk, which makes
it easier to describe the density distribution of the disk , e.g., Bovy et al.
(2016); López-Corredoira et al. (2002); López-Corredoira & Molgó (2014); Wang
et al. (2018b) and reference therein. Up to now we are still not sure the
exact flaring mechanisms and the detailed flaring features in the north and
south stellar disk. Before LAMOST (Large Sky Area Multi-Object fiber
Spectroscopic Telescope; also known as Guo Shou Jing telescope) and Gaia, we
could not get large sample for the OB stars due to some reasons such like
sample purity, sampling rates of disk region, precise distance, atmospheric
parameters, line indices, high precision photometry and so on. In this work,
we use young OB-type stars to explore spatial structure of the young thin
disk, which will help us to better understand the formation of the stellar
disk.

The structure of this paper is as follows. In Section 2, we describe how we
select OB-type stars and correct the selection effects. In Section 3, we
introduce the vertical star counts model for this work. In Section 4, we show
the results about the flaring disk and discussions. Finally, we summarize this
work in Section 5.

##  2 sample and selection effects

###  2.1 The OB-type stars catalog

LAMOST sky survey (Cui et al., 2012; Deng et al., 2012; Zhao et al., 2012) has
observed a total of 4,154 sky regions and released 9,026,365 spectra for the
phase-I. The DR5 catalog includes spectral parameters of 5,348,712 stars with
Teff, logg, and [Fe/H]. Now LAMOST is being conducted the second 5$-$year
survey since 2018 with the Medium Resolution of 7,500.

OB-type stars is composed of massive ($M>8M_{\odot}$) stars, hot sub-dwarfs,
white dwarfs, and post-AGB stars, etc. Among them, the massive stars have not
moved out of the star-forming region because of their short life. Liu et al.
(2019) identified 16,032 OB-type stars selected by spectral line indices
space, including 22,901 spectra with signal-to-noise ratio larger than 15 in
the g band, 948 hot subdwarf spectra, and 160 white dwarf spectra. The sample
in this catalog has been tested by independent human eyes one by one and thus
has high accuracy. The completeness for OB stars earlier than B7 contained in
LAMOST DR5 is better than 89 ± 22%, which is used here, it is by now the
largest spectroscopic OB-type star sample. This sample has been used in Cheng
et al. (2019); Wang et al. (2020d) for disk asymmetries.

We cross match the OB-type star catalog with Gaia DR2 (Gaia Collaboration et
al., 2018a) so that the distance of around 16,000 stars can be obtained. We
suggest the small parallax zero$-$point bias will not affect our final
conclusion due to that, our sample are actually mainly within 5 kpc away from
the sun and the parallax is larger than 0.2 mas, the small parallax
zero$-$point bias around 0.05 mas or even small can not change our conclusion,
the zero point bias increases with distance for Gaia, it is very small within
4-5 kpc (Gaia Collaboration et al., 2018b), that is our range. Even though, we
still have already finished a test for the results to prove our conclusions
are robust. The sky coverage and the number counts distributions of the sample
we used in this work are denoted in Fig. 1, respectively.

Figure 1: The top one of the figure shows normalized number counts
distribution of OB type stars without considering selection effects during
this work.The bottom one is the sky coverage of our OB-type stars adopted from
Liu et al. (2019) and the area of the sky mainly comes from the Galactic-Anti
Center direction, the sample is mainly from $l$ = [140 210] deg, $b$ = [-3 5]
deg with the relatively small contribution from other areas

.

###  2.2 The OB-type stars selection correction

We adopt the method to correct the selection effect developed by Liu et al.
(2017a) and it was tested by simulation data, this method has been used in
several works based on LAMOST giant stars, e.g., Wang et al. (2018b); Xu et
al. (2018). Before applying the selection correction, we further clarify the
OB samples by adopting the below criteria:

1) Parallax/Parallaxerror $>$ 2.5 and parallax$>$0;

2) $J$error, $K$error $<$ 0.1 mag;

3) $K$ $<$ 14.3 mag;

After applying these criteria in consideration of the parameter accuracy and
the limited magnitude or the completeness of the photometric survey, there are
13,534 OB-type stars left during this work. Then we assume that, at a given
line-of-sight of Galactic coordinates ($l$, $b$) and color magnitude box, the
probability finding a star at distance $D$ in the LAMOST survey data should be
roughly same as in the photometric 2MASS (Skrutskie et al., 2006) data. Based
on these assumptions and Bayesian theory we could obtain a selection factor,
thus then we could calculate the corrected density. Not only this method could
be used in the spatial density, but also in mono-abundance population
structures. Below we briefly introduce main principle to show how we correct
selection effects:

$p_{ph}(D|J-K,K,l,b)=p_{sp}(D|J-K,K,l,b),$ (1)

where $J-K$ and $K$ stand for the colour index and $K_{s}$-band magnitude in
2MASS (Skrutskie et al., 2006), respectively.

Then the stellar densities for photometric data, $\nu_{ph}$, and for the
spectroscopic data, $\nu_{sp}$, are associated with each other through

$\displaystyle\nu_{sp}(D|J-K,K,l,b)=$
$\displaystyle\nu_{ph}(D|J-K,K,l,b)S(J-K,K,l,b),$ (2)

where $S$ represents for the selection function of the spectroscopic data. The
selection function can be determined by

$S(J-K,K,l,b)={n_{sp}(J-K,K,l,b)\over{n_{ph}(J-K,K,l,b)}},$ (3)

where $n_{sp}$ and $n_{ph}$ are the star counts of the spectroscopic and
photometric data in the $J-K$ vs. $K$ plane from 2MASS (Skrutskie et al.,
2006), respectively. Then we could obtain the corrected stellar density by
combining the selection factor and spectroscopic number density in the
catalog.

Because the spectroscopic stars are usually very limited in a line-of-sight
and the uncertainty of distance to the stars cannot be ignored, a kernel
density estimation (KDE) is applied to derive $\nu_{sp}$ along a line-of-
sight, i.e.

$\nu_{sp}(D|J-K,K,l,b)={1\over{\Omega
D^{2}}}\sum_{i}^{n_{sp}(J-K,K,l,b)}{p_{i}(D)},$ (4)

where $p_{i}(D)$ is the probability density function of $D$ for the $i$th
star.

In general, D is the distance from the sun and we normalize the probability
for D along line of sight from 0 to infinity, the uncertainty of this method
is about 25% including extinction contribution, the size of each bin for star
counts is $\triangle$ (J$-$K) = 0.1 and $\triangle$ K = 0.25, $J-K$ vs. $K$
plane is from 2MASS. Please see more details about this method in our previous
papers (Liu et al., 2017a; Wan, Liu & Deng et al., 2017; Wang et al., 2018b;
Liu, Xu & Wang et al., 2018).

##  3 Star counts model

In this section, we will introduce how we construct the density model of the
disk used in this work. The Galactic disk mainly consists of two components of
thin and thick disks (Gilmore & Reid, 1983). However, the OB-type stars in our
sample are all young stars with lower galactic latitude. Therefore, they
almost could not contribute to the old thick disk so that we choose to only
use a single exponential vertical profile. Wang et al. (2018b) used the non-
parametric method normalizing the radial density part into the mid-plane
density, thus then they could simplify the 2D density plane to 1D density
slice, to unravel the outer disk is more complex than what we thought before.
During this work, we also use similar method but we consider the north and
south sides separately and introduce the mid-plane displacement $Z_{0}$ as an
additional free parameter, then we slice the OB-type star sample into bins
along the direction of the galactic centric distance R, which are centered at
$R$ $=$ 8.5, 9.5, 10.5, 11.5, 12.5, and 13.5 kpc. The width of each bin is 1
kpc by considering the sampling rates and poisson noise.

We use normalized exponential model to fit the stellar density distribution
assuming north-south vertical symmetry as follows:

$\nu(Z|R)=\left\\{\begin{array}[]{rcl}\nu_{0}(R)e^{-\left(\frac{Z-Z_{0}}{h_{n}(R)}\right)}&&{Z>Z_{0}}\\\
\nu_{0}(R)e^{-\left(\frac{|Z-Z_{0}|}{h_{s}(R)}\right)}&&{Z<Z_{0}}\end{array}\right.$
(5)

In which, $\nu_{0}(R)$ is the total volume density when $Z=Z_{0}$, $h_{n}(R)$
is the scale height on the north side of the young disk, $h_{s}(R)$ is the
scale height on the south side of the disk and $Z_{0}(R)$ is the true mid-
plane displacement. We assume $Z=Z_{0}$ is a free parameter in the data and
the sample is separated by $b=0$ or $Z=0$ due to that we don’t know the true
value of the mid-plane displacement and we want to explore it using our
method, then we could explore scale height of the north and south and true
mid$-$plane $Z_{0}$ simultaneously. And up to now we are also not sure
$Z=Z_{0}$ value at different distance and and it is not clear that whether the
scale height of the north and south is similar or not. In order to derive all
the unknown parameters at each $R$ bin, we firstly set up the histogram of the
mean vertical stellar density along $Z$ grid. Then we sample the likelihood
distribution as:

$\displaystyle\mathcal{L}(\\{\nu_{\rm
obs}(Z_{i}|R)\\}|\nu_{0},h_{n},h_{s},Z_{0})=$
$\displaystyle\prod_{i}\exp\left[-{1\over{2}}(\nu_{\rm obs}(Z_{i}|R)-\right.$
$\displaystyle\left.\nu_{\rm
model}(Z_{i}|R,\nu_{0},h_{n},h_{s},Z_{0}))^{2}\right],$ (6)

Where $Z_{i}$ is the $i-$th point of the $Z$ grid, $\nu_{obs}$ is the stellar
density calculated from the observation, $\nu_{model}$ is the stellar density
of the theoretical model, and corresponding errors are calculated by bootstrap
method. Here we have four free parameters in the fitting process:
$ln(\nu_{0}(R))$, $h_{n}(R)$, $h_{s}(R)$, $Z_{0}(R)$.

##  4 Results and discussions

The Markov Chain Monte Carlo (MCMC) method provided by $emcee$ (Foreman-Mackey
et al., 2013) is used for sampling the likelihood distribution. The best
fitting value of the free parameters is selected as the peak of the likelihood
distribution and the estimated uncertainty is determined by the values of 15
and 85 percentiles of the likelihood distribution. The MCMC fitting results
are shown in table 1 listing the four best fitting parameters at different $R$
bins. As an example, we show the likelihood distribution of the parameters for
8$-$9 kpc and the fitting results of model and data in Fig. 2.

Figure 2: The top panel shows the likelihood distribution of the parameters ($ln(\nu_{0})$ ($kpc^{-3}$), $h_{z}$ ( kpc), $Z_{0}$ ( kpc)) for 8$-$9 kpc drawn from the MCMC simulation, the uncertainties in the estimates are determined using the 15th and 85th percentiles of the MCMC samples. The bottom one is the fitting results and relative errors of model and data, it appears that there are some oscillations in the solar neighborhood. Table 1: Best fitting parameters of different $R$ bin with vertical density distribution models. R(kpc) | $ln\nu_{0}(kpc^{-3})$ | $h_{n}(kpc)$ | $h_{s}(kpc)$ | $Z_{0}(kpc)$
---|---|---|---|---
$8<R<9$ | 6.72±0.009 | 0.22±0.002 | 0.14±0.001 | -0.08±0.002
$9<R<10$ | 6.011±0.005 | 0.20±0.001 | 0.15±0.001 | 0.01±0.001
$10<R<11$ | 5.03±0.009 | 0.17±0.002 | 0.17±0.002 | 0.02±0.001
$11<R<12$ | 3.87±0.018 | 0.23±0.003 | 0.20±0.005 | 0.03±0.004
$12<R<13$ | 3.17±0.040 | 0.30±0.011 | 0.30±0.023 | 0.08±0.013
$13<R<14$ | 2.24±0.100 | 0.45±0.062 | 0.50±0.126 | 0.09±0.045

###  4.1 The North-South flaring disk compared with stellar disk

The scale height of the young stellar disk is shown in Fig. 3. As we could
see, the north and south scale height of the young stellar disk is increasing
from 0.14 to 0.5 kpc in the range of 8$-$14 kpc, that is the clear flaring
signatures, so we could say we detect the features of the flaring in the young
stellar disk.Meanwhile, it is shown that the pattern of the north and south
disk is similar or slightly different, the average difference of the scale
height in the north and south is 19 pc and the slope of both sides has no
large difference, suggesting us that the flaring of the disk might be
symmetrical. During this work we use average slope of the scale height vs.
distance to describe the flaring strength, which has been used in Wan, Liu &
Deng et al. (2017).

Previous works are also compared in the figure and we could see, the thickness
of the young stellar disk is similar to or slightly different from the thin
disk traced by RGB stars in Wang et al. (2018b). The clear difference between
the young stellar disk of OB stars and the thin disk or disk distribution of
López-Corredoira et al. (2002); López-Corredoira & Molgó (2014) are also
displayed here. The reason for the difference could be naturally understood
that the different populations and methods are used in these works, e.g., here
we use the young OB-type stars and non-parametric method, López-Corredoira et
al. (2002) used the photometric red clump stars without considering the two
components of the disk, and López-Corredoira & Molgó (2014) adopted the
F8V-G5V stars to study the changes of the scale height of thin and thick disks
with different methods. Importantly, we detect the almost symmetrical flaring
signature in the young stellar disk of OB samples, and it is similar with or
even slightly stronger than the old stellar thin disk traced by RGB samples
for slope value. Please also notice that our OB star sample is distributed in
the Galactic-Anti Center and we just use the projection distance, that is to
say, $R$ and $Z$, and these compared data in Fig. 3 are almost originated from
the same regions, mainly from the north hemisphere. There are some clear
different regions for gas disk but we could still compare the disk scale
length and scale height.

Figure 3: The figure shows the OB-type stars and other tracers scale height
distribution along the radial distance, red pentagram line is northern OB-type
stars scale height, black pentagram line is the southern OB-type stars scale
height. Green triangle line is the thin disk scale height of Wang et al.
(2018b) and green rectangle line is the thick disk scale height of Wang et al.
(2018b). Brown and blue line are the thin disk and thick disk of López-
Corredoira & Molgó (2014) respectively. Yellow line is the results of López-
Corredoira et al. (2002). These data compared here are originated from the
same regions and the average slope describing the strength of flaring is
labeled in the bottom left, different colors are corresponding to different
results.

###  4.2 The North-South flaring disk $h_{z}$ compared with gas disk

The scale height of the young stellar disk compared with HI and H2 gas disk of
Kalberla et al. (2007); Nakanishi & Sofue et al. (2006) is shown in Fig. 4,
for Nakanishi & Sofue et al. (2006) the data is mainly covered by $b$ = [-1.5
1.5∘] and $l$ = [0 to 90 ∘, 270 to 360 ∘] and for Kalberla et al. (2007) the
data range is covered by all sky disk region by discarding $b$ $>$ 30∘ and
Galactic distance $R$ = [-3.5 3.5 kpc]. It is shown that the young stellar
disk, for the overall trend, is slightly thicker than the gas disk but it is
not significantly different for the scale height, which is expected since that
we declare that OB-type stars could be used to traced more gas properties and
dynamics, these will be shown in our series of works. The OB-type stars are
very young and have not moved out of the star formation area so that there is
no doubt the structure of the OB-type stars are inherited some properties of
the gas clouds that are forming stars. We also find that the HI disk has some
differences with the H2 disk due to the temperature difference possibly. The
distance of our OB star sample is from Gaia, the extinction is very small and
we could use its precise distance within 8-14 kpc. If we agree there are some
influence on our results caused by extinction and not perfect selection
effects, then the error bar will be enlarged, thus then the true difference
will possibly be reduced so that strengthen our points to some extent, that
is, we could use OB star to infer some properties of the gas disk, then it
will encourage us to work more.

The most important tracer for molecular gas in galaxies is the CO spin line of
different isotopes, especially the 2.6 mm line of CO and high-sensitivity CO
observations are essential for describing the structure of the disk in
galaxies, sometimes we also choose to convert CO to H2. In Sun et al. (2019),
they also confirmed the gas thin disk (FHWM $\approx$ 90 pc) and revealed the
gas thick disk was around 280 pc consistent with the HI gas disk scale height
(250$-$300 pc), meanwhile they pointed out molecular gas properties might be
related to massive stars evolution, which support our points that the gas disk
properties could be inherited by the young OB-type stars so that the young
stellar disk could be compared with gas disk to investigate more intriguing
details. Recently, Su et al. (2021) find that the molecular gas disk has two
components, gas thin disk and thick disk, with scale height of 85 pc and 280
pc in the region of $l$ from 16∘ to 52∘ and $b$ from -5.1 to 5.1∘,
respectively, these values are belonging to our range of OB stellar disk. We
will have a few more discussions in the latter part.

Figure 4: The figure shows the OB-type stars scale height distribution,
compared with the gas disk, along the radial distance, red line is north side
scale height, black line is the southern scale height. Green line is the
Kalberla et al. (2007) southern scale height of the HI gas disk, the blue
rectangles and pentagram are also the Kalberla et al. (2007) northern HI
thickness distribution, yellow is corresponding to the H2 results in Nakanishi
& Sofue et al. (2006). The average slope of the scale height (hz) vs. distance
(R) describing the strength of flaring is labeled in the bottom left.

###  4.3 The scale length of the OB-type stellar disk compared with gas disk

The stellar density of the disk, general speaking, is exponentially
distributed in the radial direction and gradually decreases with the increase
of the distance and we fit the logarithm of the stellar density in the
direction of the Anticenter belonging to the radial range of R = 8$-$14 kpc.
The scale length of the OB-type star is 1.17 $\pm$ 0.05 kpc and the comparison
with previous works of the H2 gas disk which has been transferred from CO data
(Dame et al., 1987; Bronfman et al., 1988; Digel et al., 1991; Nakanishi &
Sofue et al., 2006; Pohl et al., 2008) is shown in Fig. 5, it should be
pointed out that, for gas density, we directly use the data points from other
literatures. Some works have number density with units cm-3 and this gas
density is often calculated by column density divided by radius size, however,
here we have actually used the surface mass density (M☉pc-2 )which is also
suitable for us to compare scale length. In order to compare with stellar disk
in one figure, the gas density value and stellar density value of Wang et al.
(2018b) is shifted up in y-axis. It shows that the young stellar disk scale
length is shorter than the gas disk when comparing the value labeled in the
figure except the work of Bronfman et al. (1988) with only three points in a
narrower range so it is not reasonable, others are implying that generally the
gas disk is more extended than the stellar disk. We could also attempt to
deduce from here stellar disk might be more compact, similar to the hot dust,
than the gas and cold dust disk. For the gas disk region mentioned in this
part, the range of Dame et al. (1987) is covered by 10 $-$ 20∘ in latitude at
all longitudes and all or nearly all large, nearby clouds at higher latitude.
Latitude is from -2 to 2∘ and longitude is from 300 to 348∘ is corresponding
to the region of the Bronfman et al. (1988), but we only use the north data
beyond 8 kpc. Based on the smoothed$-$particle hydrodynamical simulations and
the complied data (Dame et al., 2001), Pohl et al. (2008) made full use of the
CO survey of which total area covered by the composite survey is 9,353 deg2
$-$ more than one$-$fifth of the entire sky and nearly one half of the area
within 30∘ of the Galactic plane. Digel et al. (1991) used data located in the
outer galaxy within the disk region from $l$ =[65 116∘]. Here H2 gas surface
density is converted from CO surveys, using uniform assumptions regarding the
Galactic rotation curve, solar radius, and the CO$-$to$-$H2 conversion factor,
all these data are from the review (Heyer et al., 2015).

In addition to that, we notice that the scale length of the OB-type stellar
disk here is different from the results given by Li et al. (2019), it is
$2.10\pm 0.01$ kpc and scale height is from 132$-$450 pc, both are
morphological or geometrical definition, we suggest it is due to that Li et
al. (2019) used Gaia DR2 and 2MASS photometric data and it must include the
early type and late-type of OB-type stars, what we want to emphasize here, our
sample is lacking of late type ones and mainly consisted of early type stars
caused by the selection methods, so the scale length has some differences for
these two works.

Figure 5: The figure shows mid-plane density distribution along the distance,
the linear fitting results are shown as dashed lines with different colors,
the scale length of different works for H2 gas disk (Dame et al., 1987;
Bronfman et al., 1988; Digel et al., 1991; Nakanishi & Sofue et al., 2006;
Pohl et al., 2008) is labeled in the bottom left, the H2 value is converted
from CO (Heyer et al., 2015). The left y axis is corresponding to the stellar
disk and the right axis is the gas surface density. The gas density is
adjusted by shifting up 5.5 and the Wang et al. (2018b) is shifted up 10 in
order to compare in one figure.

###  4.4 The mid-plane displacements of the disk

The gas warp was detected as shown in Kerr (1957); Bosma (1981); Briggs et al.
(1990). Recently, with the help of classical Cepheids, Chen et al. (2019) have
revealed an intuitive three-dimensional (3D) map of the stellar warp by star
counts methods, it is covering the all sky until 20 kpc. Meanwhile, Skowron et
al. (2019a, b) have also showed us number density stellar warp signals in
different distance and azimuthal angles located in the north and south sky and
declared that the amplitude of its northern part is very prominent and
stronger than that of the southern part. Meanwhile, the kinematical and height
signal of the stellar warp and gas infall origins are also unraveled in Wang
et al. (2020d) with the help of sample mainly from the LAMOST Galactic-Anti
Center in the north hemisphere. Some other works about warp signals with
different tracers and methods could also be found in Poggio et al. (2018) by
using Gaia sample mainly in the north side, then Poggio et al. (2020) used all
sky Gaia 12 million giant stars to detect the warp precession for the first
time. However, Chrobáková et al. (2021) recently recalculated the warp
precession and found that there is no need for precession using different warp
parameters but same approach and Gaia DR2 kinematic data.

During this work we introduce the mid$-$plane displacement ($Z_{0}$) to
Equation 5 and Equation 3 to do the decomposition and fit the OB-type stars in
radial slice displayed in Fig. 6. Intriguingly, we find that the mid-plane
displacements are existing and different locations have different $Z_{0}$
values, but almost are within 100 pc accompanying by that the value is
approximate to be 0 around solar location in the range of 8-9 kpc.

More importantly , we find the offset is increasing with the distance and
infer that there is a possibility it might reflect suspected signals like the
warp here. If it is true we could say this work detects the similar warp
signal in the young stellar population and warp should happen in the early
time. However, here from only one figure, we could not rule out other
possibilities due to the limited range of the sample and other possible
perturbations. In short, the pattern might be caused by warp, or caused by the
heating of some perturbations which are not the same with warp dynamical
mechanisms, sometimes we actually could not discriminate whether the vertical
signal is from warp or external perturbations (Wang et al., 2020a, c) or other
unknowns, which is also not the main target of this paper. Please notice that
the zero point bias might affect the exact value here but it will not change
the pattern shown in this work. The significance of the signal is still needed
to explore more in the future, but we just don’t want to rule out some
possibilities here, our motivation is flaring.

By using the similar method introduced in Xu et al. (2015), Wang et al.
(2018b) proposed that some disk oscillations are very complicated but could
still be explained, to some extent, by shifting either the thin or the thick
disc, so during this work we do like this in the young stellar disk is
reasonable and the value of the mid$-$plane displacement ($Z_{0}$) is also not
strange at all. In general, our mid-plane shifting pattern is actually similar
with that in Fig. 5 of recent work finished by Xu et al. (2020) based on the
kinematics method.

Figure 6: The figure displays the mid-plane displacement distribution along
the $R$, almost are within the 100 pc and the slope shown in the bottom left
is 0.020.please notice that the slope is just for quantitative analysis, which
is different from flaring strength.

###  4.5 Discussions

####  4.5.1  Comparisons for the disk w/o North and South sides

In order to test the effect of the method with decomposition of the north and
south sides, we have finished a comparison for the disk features between the
north-south scale heights and the values without considering north-south
sides. As shown in Fig. 7, the scale heights, scale lengths and mid-plane
offsets are compared from the top to the bottom. For the scale heights, there
are no large differences and the general pattern is the same with each other
and similarly, the mid-plane offsets are also showing the same trend with some
differences. Interestingly, the mid-plane density and the scale length are
matched very well. All these show that our method with north and south sides
will not change our main conclusion for this work, especially the flaring.

Figure 7: Comparison for the disk in consideration of north and south sides
and the one without considering the scale heights of north and south sides
simultaneously, which are colored by red or black pentagram and green square
respectively. The top one is the scale height, the middle one is the mid-plane
density and, the slope of the scale height and scale length along the distance
is labeled. The bottom one is the mid-plane offsets with its slope label,
clearly, the pattern is stable.

By combining all stars belonging to 8-14 kpc without radial slices, we have
also finished the fitting results for the young stellar disk, the scale height
of the young stellar disk is 0.16 kpc and the mid-plane offset of the disk is
11 pc, the fitting results are quite well and shown by us in Fig. 8.

Figure 8: The panel shows the likelihood distribution of the parameters
($ln(\nu_{0})$ ($kpc^{-3}$), $h_{z}$ ( kpc), $Z_{0}$ ( kpc)) for 8$-$14 kpc
data drawn from the MCMC simulation without considering north and south sides
simultaneously, the uncertainties in the estimates are determined using the
15th and 85th percentiles of the MCMC samples.

####  4.5.2 Possible scenarios for the flaring

With the release of Gaia DR2 data, Antoja et al. (2018) found phase mixing
features of the Milky Way disk possibly due to the disturbance caused by the
vibration of the dwarf galaxy passing through the disk 300-900 Myr ago. The
lifetime of OB stars is mostly $<100$ Myr and the gas disk forming the OB-type
star could possibly be disturbed by the dwarf galaxy so that some imprints,
might be inherited by the OB-type stars (Cheng et al., 2019). Cheng et al.
(2019) also introduced two types of interference sources: First, internal
disturbances, such as spiral structures, huge molecular clouds, or central
rods; The second, is external disturbances, such as dwarf galaxy or satellites
passing through the Milky Way. In addition to the aforementioned disturbances,
Cheng et al. (2019) proposed another case, that is saying, a few decades ago,
the dark matter sub-halo perturbed the disk and left marks on these young
stars. Similar perturbation scenarios to the disk such as warp, stream,
satellites, dwarf galaxies, dark matter sub-halos are also mentioned in Wang
et al. (2018a, 2019, 2020a, 2020d); Khanna et al. (2019); Bland-Hawthorn et
al. (2019); López-Corredoira et al. (2020); Xu et al. (2020) and reference
therein.

Flaring is possibly related to the disk heating, there are many mechanisms
proposed in the past many years. e.g. Spiral arms could be the important
contributor to the in-plane motions but without influence on the vertical
motions (Sellwood et al., 1984, 2013). Giant Molecular Clouds (GMC) could heat
the disk and redirect the velocity of stars out of the plane (Lacey et al.,
1984). Bar might not be important beyond solar neighborhood but can not be
ignorable in the inner bulge region (Moetazedian et al., 2016). Migration
might heat the disk and contribute to the flaring, which was mentioned in
recent work (Bovy et al., 2016), but the different scenarios proposed by
Minchev et al. (2012, 2018) suggest that migration is not important for
heating and flaring. External perturbations such like dwarf galaxy was heating
the disk (House et al., 2011) and is the mechanism for the flaring
(Kazantzidiset al., 2008), in addition, mergers and perturbations from
satellite galaxies and subhalos can induce radial mixing in the outer disk
(Quinn et al., 2009). For our work, we focus on the average slope of scale
height vs. distance which is used to describe the strength of the flaring
features in different populations, then we could attempt to infer which
mechanism of the flaring is more important, we agree the different heating
mechanisms have different size of the thickness but it is not our target of
the paper, hopefully in the future we could combine kinematics and star counts
to pursue more details about this intriguing question.

Here again, as we mentioned in the introduction, we attempt to classify the
mechanisms of the flaring into two types: one is the secular evolution of the
disk, which expects the younger population has weaker flaring strength than
the older population, the other one we think it is perturbation caused by many
agents, which don’t expect the population evolution effect with age for the
origins of the flaring. Against this big background then we could discuss
which one is more important, we also know the physical details might be more
complicated and we are not planing to rule out other physical scenarios only
from current analysis.

Therefore, from current results, we clearly detect the flaring strength of
young OB stellar population is similar with or even slightly stronger than the
old RGB population by focusing on the thin disk by using the slope of the
scale height vs. distance to describe the strength of the flaring, with the
average value of the slop is 0.042 (north-OB) vs. 0.065 (south-OB) vs. 0.020
(RGB), the qualitative conclusion of which is similar with Wang et al.
(2018b), they also found the flaring strength of Red Clump Giant stars (RCG)
is similar with the RGB pattern for the thin disk. By comparing flaring
strength of the thin disk and thick disk, both Wan, Liu & Deng et al. (2017)
and Wang et al. (2018b) detected the similar flaring features for these two
different populations. That is to say, we suggest that the secular evolution
is not the main contributor during our works and for the quantitative analysis
of Wang et al. (2018b) shown in this paper, it is the same value of the slope
for the thin and thick disk, as shown in Fig. 3.

And for the current work, if we look carefully the Fig. 3 by only focusing on
the thin disk stellar population, the purple solid line is the F/G type stars,
the green line with triangles is K-type/RGB stars, the red and black line with
pentagram is OB-type stars, it is not possible for us to see the trend that
the older the populations, the stronger the flaring by focusing on the
variation of the scale height or the slope of the scale height with distance,
the quantitative analysis for the strength of the flaring is that the OB type
is 0.042 or 0.065, F/G type is 0.145, K type is 0.020, which is also not
consistent with the scenario of the secular evolution. Meanwhile, if we
concentrate on the thick disk stellar population, the blue line is F/G type
stars, the green line with rectangles is RGB/K type stars, the former one is
generally younger than the latter one, but it appears that both have similar
flaring features with the slope is 0.023 and 0.020 separately, it is hard to
be explained by secular evolution scenarios. All these evidence again support
that the secular evolution predicting the flaring strength is different in
different populations with that older stellar population should have stronger
flaring, is not the main contributor, although we could not completely rule
out the contributions of this scenario up to now. In other words, the
perturbations such as dwarf galaxy might be more possible according to current
two classifications.

Recently, based on the smooth particle hydrodynamical simulations with clumpy
episodes, Beraldo e Silva et al. (2020) discovered broken exponentials is
existing in the low-[$\alpha$/Fe] population and more importantly, flare
significantly for low-[$\alpha$/Fe] stars but only weakly for
high-[$\alpha$/Fe] populations. Meanwhile they also found that for
low-[$\alpha$/Fe] thin disk stars, the flaring level decreases with age, which
is also contrary to the secular evolution and it is similar with GMC
perturbation scenario.

With the help LAMOST K giant stars and two types of simulations, Xu et al.
(2020) pointed out the external perturbations might influence the outer disk
beyond solar neighborhood and revealed that the flaring could be explained by
perturbations qualitatively, which supports our viewpoints here about the
flaring. Combining young OB type Stars and old K giant stars we could infer
that the sensitive time of the flaring features to the possible perturbations
could be from very early time in OB population to the relatively late time in
K giant population.

Although we are far from the exact flaring mechanisms, we suggest here,
possible perturbations such like the mentioned above could heat the disk and
then, cause the disk flaring, therefore we could detect that the disk displays
clear flaring signatures. More importantly, during this work we discover these
possible heating mechanism might be thickening the disk of both sides in some
similar ways, especially in the outer disk, this will push us to have better
knowledge of the well$-$known flaring disk, especially for theorists.

####  4.5.3  Relations between the young stellar disk and gas disk

We also find the evidence of relationship between the young OB-type stars and
the gas disk by comparing their structure parameters. It is shown in Fig. 4
that the scale height of the OB-type stars is slightly thicker than the gas
thickness by 9 pc in the south and 33 pc in the north for the median value,
the small difference supports the truth that OB stars is very young and only
last a few million years and during this time the stars should not have had
time to significantly move from their birth locations, so the young stars and
gas could share similar properties and we could use young OB stars to trace
some physics about gas. For the reason that the OB stellar disk is a little
thicker than the gas disk, we speculate there are two reasons, one is that the
stars are hotter than gas so it is thicker, the other is the stars are
perturbed by some agents and they are possibly more effective for stars, the
latter one is also consistent with the perturbations scenarios if we assume
that the OB stellar disk is heated by perturbations after it formed from the
gas. As far as we know, we don’t see some other works show the similarities of
the young stellar disk and gas disk so that it might encourage our community
to explore more about these. We are planing to investigate more and farther
for stellar, gas, and dust disk. For this work, due to the sampling rate is
low when the distance is larger than 14 kpc so we only focus a limited range
here.

In short summary, we detect the clear flaring signatures of the young stellar
disk that is possibly mainly contributed by perturbations, this symmetrical
flaring evidence for young stellar disk is clearly of vital importance for us
to understand flaring mechanisms and even gas disk physics.

##  5 Summary

In this work, we explore the spatial structures of young stellar disk with the
LAMOST DR5 OB-type stars catalog, based on the Bayes method for selection
effect correction and non-parametric method for disk decomposition, both are
robustly tested in our previous works. We discover that the north-south
flaring in the outer disk is symmetrical and the scale height is from 0.14 to
0.5 kpc at different distance from 8 to 14 kpc. The flaring strength of the
younger stellar disk is similar to or slightly stronger than the results of
the thin disk traced by older red giant branch stars with the average slope
charactering the strength is 0.042 or 0.065 (OB) and 0.020 (RGB) separately,
which is supporting that the flaring might not be mainly contributed by the
secular evolution scenarios, which are predicting the flaring of the old
populations should be stronger than the younger ones. And meanwhile, we
speculate that the mechanism might have similar influence on the both disk
sides in some similar ways.

We also find that the OB stellar disk is slightly thicker than the gas disk by
9 pc in the south and 33 pc in the north for the median value, implying that
we could attempt to, to some extent, use OB-type stars to infer the gas
properties in the future and what’s more, we unravel that the stellar disk
might be more compact than the gas disk by showing that the scale length of
OB-type stars is shorter than that of gas. We believe that more comparisons of
the young stellar disk and gas disk will be non-trivial to uncover more
secrets of the formation of the disk.

Finally, it is displayed that the mid-plane displacements($Z_{0}$) of the thin
disk have a clear increasing trend with distance from 8-14 kpc, that is
possibly implying the signal like warp, but so far we are still not sure its
origins and significance due to that we can not discriminate the warp
contributions and perturbations contributions, or other unknowns from the
current limited range and analysis. All these observational evidences of the
young stellar thin disk shown here will encourage us to reveal more about the
flaring mechanisms and gas disk physics, more works will be shown in the
future.

We would like to thank the anonymous referee, Lixia Yuan, Yang Su, Miaomiao
Zhang, Ningyu Tang, Yan Sun, Ivan Minchev, Leandro Beraldo e Silva for his/her
very helpful and insightful comments and discussions.This work is supported by
the National Key Basic R&D Program of China via 2019YFA0405500,
2019YFA0405501, the National Natural Science Foundation of China (NSFC) under
grant 11773009, 11873057, 11390371, 11673007, 12003027. W.Y.C. is also
supported by the ”333 talents project” of Hebei Province. H.F.W. is also
supported by the LAMOST Fellow project, funded by China Postdoctoral Science
Foundation via grant 2019M653504 and 2020T130563, Yunnan province postdoctoral
Directed culture Foundation, and the Cultivation Project for LAMOST Scientific
Payoff and Research Achievement of CAMS-CAS. The author gratefully
acknowledges the support of K.C. Wong Education Foundation. The Guo Shou Jing
Telescope (the Large Sky Area Multi-Object Fiber Spectroscopic Telescope,
LAMOST) is a National Major Scientific Project built by the Chinese Academy of
Sciences. Funding for the project has been provided by the National
Development and Reform Commission. LAMOST is operated and managed by National
Astronomical Observatories, Chinese Academy of Sciences. This work has also
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

  * Antoja et al. (2018) Antoja T., Helmi A., Romero-Gómez M., et al. 2018, Nature, 561, 360
  * Bland-Hawthorn & Gerhard (2016) Bland-Hawthorn J., Gerhard O., 2016, ARA&A, 54, 529
  * Bland-Hawthorn et al. (2019) Bland-Hawthorn, J., Sharma, S., Tepper-Garcia, T., et al. 2019, MNRAS, 486, 1167
  * Bovy et al. (2012) Bovy J., Rix H.-W., Liu C., et al. 2012, ApJ, 753, 148
  * Bovy et al. (2016) Bovy, J., Rix, H.-W., Schlafly, E, F., et al. 2016, APJ, 823, 30
  * Bronfman et al. (1988) Bronfman L, Cohen RS, Alvarez H, May J, Thaddeus P. 1988, ApJ, 324, 248
  * Bosma (1981) Bosma, A. 1981, AJ, 86, 1825
  * Briggs et al. (1990) Briggs, F. H., et al.,1990, APJ, 352,15
  * Beraldo e Silva et al. (2020) Silva L. B., Debattista V. P., Khachaturyants, David Nidever D., 2020, MNRAS, 492, 4716
  * Chen et al. (2017) Chen B.-Q., et al. 2017, MNRAS, 464, 2545
  * Cheng et al. (2019) Cheng X. L., Liu C., Mao S. D., et al. 2019, ApJ, 872, L1
  * Cui et al. (2012) Cui, X.-Q., Zhao, Y.-H., Chu, Y.-Q., et al. 2012, RAA, 12, 1197
  * Chen et al. (2019) Chen, X., Wang, S., Deng, L., et al. 2019, Nature Astronomy, 3, 320
  * Carlberg & Sellwood et al. (1985) Carlberg R. G., Sellwood J. A., 1985, ApJ, 292, 79
  * Chrobáková et al. (2021) Chrobáková, Z., López-Corredoira, M., 2021, ApJ in press, arXiv:2105.04348v1
  * Dame et al. (1987) Dame T M, Ungerechts H, Cohen RS, et al. 1987, ApJ, 322, 706
  * Dame et al. (2001) Dame, T.M., Hartmann, D., and Thaddeus, P., 2001, ApJ, 547, 792
  * Digel et al. (1991) Digel S. W. 1991. Molecular clouds in the distant outer galaxy. PhD Thesis, Harvard Univ., Cambridge, MA
  * Deng et al. (2012) Deng, L. C., Newberg, H. J., Liu, C., et al. 2012, Research in Astronomy and Astrophysics, 12, 735
  * Foreman-Mackey et al. (2013) Foreman-Mackey D., Hogg D. W., Lang D., Goodman J., 2013, PASP, 125, 306
  * Gilmore & Reid (1983) Gilmore G., Reid N., 1983, MNRAS, 202, 1025
  * Grand et al. (2016) Grand, Robert J. J., Springel, V., Gómez, Facundo A. et al. 2016, MNRAS, 459, 199
  * Gaia Collaboration: Prusti et al. (2016) Gaia Collaboration, Prusti, T., de Bruijne, J. H. J., et al. 2016, A&A, 595, A1
  * Gaia Collaboration et al. (2018a) Gaia Collaboration, Brown, A. G. A., Vallenari, A., et al. 2018a, A&A, 616, A1
  * Gaia Collaboration et al. (2018b) Gaia CollaborationGaia Collaboration, Katz, D., Antoja, T., et al. 2018b, A&A, 616, A11
  * Huang et al. (2018) Huang, Y., Schönrich, R., Liu, X.-W., et al. 2018, ApJ, 864, 129
  * House et al. (2011) House E. L. et al., 2011, MNRAS, 415, 2652
  * Heyer et al. (2015) Heyer M., Dame T. M., 2015, ARA&A, 53, 583
  * Jurić et al. (2008) Jurić M., et al. 2008, ApJ, 673, 864
  * Jenkins & Binney et al. (1990) Jenkins, A., & Binney, J. 1990, MNRAS, 245, 305
  * Jenkins et al. (1992) Jenkins, A. 1992, MNRAS, 257, 620
  * Kalberla et al. (2007) Kalberla P. M. W. , Dedes L., Kerp J., & Haud U., 2007, Astronomy & Astrophysics, 469, 511
  * Khanna et al. (2019) Khanna, S., et al. 2019, MNRAS, 489, 4962
  * Kazantzidiset al. (2008) Kazantzidis, S., Bullock, J. S., Zentner, A. R., et al. 2008, ApJ, 688, 254
  * Kerr (1957) Kerr F. J., 1957, AJ, 62, 93
  * Liu et al. (2012) Liu, C., $\&$ van de Ven, G. 2012, MNRAS, 425, 2144
  * Liu et al. (2017a) Liu, C., Xu Y., Wan J. C., et al. 2017, RAA, 17, 096
  * Liu et al. (2019) Liu, Z., Cui, W., Liu, C., et al. 2019, ApJS, 241, 32
  * Liu et al. (2017c) Liu C., Tian H. J., Wan J. C., 2017c, arXiv:1702.02233v1
  * Liu, Xu & Wang et al. (2018) Liu, C., Xu, Y., Wang, H.-F., et al. 2018, in IAU Symp. 334, Rediscovering our Galaxy, ed. C. Chiappini et al. (Cambridge: Cambridge Univ. Press), 109
  * Li et al. (2019) Li, C.-D., Zhao G., Jia Y. P.,et al. 2019, APJ, 871, 208
  * López-Corredoira et al. (2002) López-Corredoira, M., Cabrera-Savers, A., Carzón,F., et al. 2002, A$\&$A, 394, 883
  * López-Corredoira & Molgó (2014) López-Corredoira M., Molgó J., 2014, A&A, 567, A106
  * López-Corredoira et al. (2014) López-Corredoira, M., Abedi, H., Garzón, F., & Figueras, F. 2014, A$\&$A, 572, 101
  * López-Corredoira et al. (2020) López-Corredoira, M., Garzón, F., Wang, H.-F., et al. 2020, A$\&$A, 634, A66
  * Laporte et al. (2018) Laporte, C. F. P., Johnston, K. V., Gómez, F. A., et al. 2018, MNRAS, 481, 286
  * Lacey et al. (1984) Lacey C. G., 1984, MNRAS, 208, 687
  * Minchev et al. (2012) Minchev, I., Famaey, B., Quillen, A. C., Dehnen, W., Martig, M., Siebert, A., 2012, A&A, 548, 127
  * Minchev et al. (2014) Minchev, I., Chiappini, C., Martig, M., Steinmetz, M., de Jong, R., et al. 2014, ApJL, 781, L20
  * Minchev et al. (2015) Minchev, I., Martig, M., Streich, D., Scannapieco, C., de Jong, R. S., Steinmetz, M., 2015, ApJL, 804, 9
  * Minchev et al. (2016) Minchev, I., Chiappini, C. and Martig, M. 2016, AN, 337, 944
  * Minchev et al. (2018) Minchev, I., F. Anders, A. Recio-Blanco, C. Chiappini, P. de Laverny, A. Queiroz, M. Steinmetz, V. Adibekyan, et al. 2018, MNRAS, 481, 1645
  * Ma et al. (2017) Ma, X. C., Hopkins, Philip F., Wetzel, Andrew R., et al. 2017, MNRAS, 467, 2430
  * Martig et al. (2016) Martig M., Minchev I., Ness M., Morgan F., Hans-Walter R., 2016, ApJ, 831, 139
  * Moetazedian et al. (2016) Moetazedian R., Just A., 2016, MNRAS, 459, 2905
  * Nakanishi & Sofue et al. (2006) Nakanishi & Sofue et al. 2006,PASJ, 58, 847
  * Narayan & Jog et al. (2006) Narayan, C. A. & Jog, C. J. 2002, A&A, 394, 89
  * Pohl et al. (2008) Pohl M., Englmaier P., Bissantz N., 2008, ApJ, 677, 283
  * Poggio et al. (2018) Poggio, E., et al., Mon. Not. R. Astron. Soc., 2018, 481, L21
  * Poggio et al. (2020) Poggio, E., Drimmel, R., Andrae, R., et al. 2020, NatAs, 4, 590
  * Quinn et al. (1993) Quinn, P. J., Hernquist, L., & Fullagar, D. P., 1993, ApJ, 403, 74
  * Quinn et al. (2009) Quillen, A. C., Minchev, I., Bland-Hawthorn, J., and Haywood, M., 2009, MNRAS, 397, 1599
  * Rix & Bovy et al. (2013) Rix H.-W., Bovy J., 2013, A&ARv, 21, 61
  * Solway et al. (2012) Solway, M., Sellwood, J. A., & Schönrich, R., 2012, MNRAS, 422, 1363
  * Skrutskie et al. (2006) Skrutskie M. F. et al., 2006, AJ, 131, 1163
  * Schönrich et al. (2017) Schönrich, R., & Dehnen, W., 2018, MNRAS, 478, 3809
  * Sun et al. (2019) Sun, Y., Yang, J., Zhang, S. B., et al. 2019, ApJS, 240, 9
  * Skowron et al. (2019a) Skowron, D. M., Skowron, J., Mróz, P., et al. 2019a, Science, 365, 478
  * Skowron et al. (2019b) Skowron, D. M., Skowron, J., Mróz, P., et al. 2019b, Acta Astron., 69, 305
  * Skrutskie et al. (2006) Skrutskie M. F. et al., 2006, AJ, 131, 1163
  * Sellwood et al. (1984) Sellwood, J. A., & Carlberg, R. G. 1984, ApJ, 282, 61
  * Sellwood et al. (2013) Sellwood J. A., 2013, ApJ, 769, L24
  * Su et al. (2021) Su, Y., Yang, J.,Yan, Q.-Z., et al. 2021, in press, arXiv:2102.10727v2
  * van der Kruit (1988) van der Kruit P. C., 1988, A&A, 192, 117
  * van der Kruit & Freeman (2011) van der Kruit P. C., Freeman K. C., 2011, ARA&A, 49, 301
  * Villalobos et al. (2008) Villalobos, A., & Helmi, A., 2008, MNRAS, 391, 1806
  * Wang et al. (2018a) Wang, H. F., López-Corredoira, M., Carlin, J. L., et al. 2018a, MNRAS, 477, 2858
  * Wang et al. (2018b) Wang, H. F., Liu, C., Xu, Y., et al. 2018b, MNRAS, 478, 3367
  * Wang et al. (2019) Wang, H. F., Carlin, J. L., Huang, Y., et al. 2019, ApJ, 884, 135
  * Wang et al. (2020a) Wang, H. F., López-Corredoira, M., Huang, Y., et al. 2020a, MNRAS, 491, 2104
  * Wang et al. (2020c) Wang, H. F., López-Corredoira, M., Huang, Y., et al. 2020c, ApJ, 897, 119
  * Wang et al. (2020d) Wang, H. F., Huang, Y., Zhang, H. W., et al. 2020d, ApJ, 902, 70
  * Wan, Liu & Deng et al. (2017) Wan J.-C., Liu C., Deng L.-C., et al. 2017, RAA, 17, 079
  * Xu et al. (2015) Xu, Y., Newberg, H. J., Carlin, J. L., Liu, C., Deng, L., Li, J., Schönrich, R., & Yanny, B., 2015, ApJ, 801, 105
  * Xu et al. (2018) Xu Y., Liu C., Xue X. X., et al. 2018, MNRAS, 473, 1244
  * Xu et al. (2020) Xu Y., Liu C., Tian H., et al. 2020, ApJ, 905, 6
  * Yu & Liu (2017) Yu J., Liu C., 2017, MNRAS, 475, 1093
  * Zhao et al. (2012) Zhao, G., Zhao, Y.-H., Chu, Y.-Q., Jing Y.-P., Deng, L.-C., 2012, RAA, 12, 723