# ReF - Rotation Equivariant Features for Local Feature Matching

Abhishek Peri∗1, Kinal Mehta∗1, Avneesh Mishra1, Michael Milford2,
Sourav Garg2, K. Madhava Krishna1 *Denotes authors with equal
contribution1Robotics Research Center, IIIT Hyderabad.2QUT Centre for
Robotics, Queensland University of Technology (QUT), Australia.This research
is supported by MathWorks.Code link: https://github.com/abhishek-peri/ReF-
official-code

###### Abstract

Sparse local feature matching is pivotal for many computer vision and robotics
tasks. To improve their invariance to challenging appearance conditions and
viewing angles, and hence their usefulness, existing learning-based methods
have primarily focused on data augmentation-based training. In this work, we
propose an alternative, complementary approach that centers on inducing bias
in the model architecture itself to generate ‘rotation-specific’ features
using Steerable E2-CNNs, that are then group-pooled to achieve rotation-
invariant local features. We demonstrate that this high performance, rotation-
specific coverage from the steerable CNNs can be expanded to all rotation
angles by combining it with augmentation-trained standard CNNs which have
broader coverage but are often inaccurate, thus creating a state-of-the-art
rotation-robust local feature matcher. We benchmark our proposed methods
against existing techniques on HPatches and a newly proposed UrbanScenes3D-Air
dataset for visual place recognition. Furthermore, we present a detailed
analysis of the performance effects of ensembling, robust estimation, network
architecture variations, and the use of rotation priors.

##  I INTRODUCTION

Local feature matching plays a central role in many computer vision and
robotics tasks such as visual place recognition [1, 2, 3], image stitching
[4], and 6-DoF pose estimation [5, 6, 7, 8]. Traditional approaches for
detection and description such as SIFT [9] and ORB [10] first detect the
points of interest and then use the patch around the pixel/keypoint to
generate a descriptor. These approaches fail under severe illumination and
rotation changes. Recent learning-based approaches such as SuperPoint [11],
SuperGlue [6], and RoRD [7] are trained as a ‘general’ system for handling
such challenges using training data that is augmented heavily with
illumination and viewpoint variations.

In this work, we take an alternative but a complementary approach that induces
bias in the model architecture itself to generate ‘rotation-specific’
features. For this purpose, we make use of Steerable CNN [12, 13] and propose
ReF (Rotation Equivariant Features) network which generates feature
descriptors that are group pooled to achieve rotation-invariant local feature
matching. We then show that highly accurate but rotation-specific coverage of
ReF can be expanded to all rotation angles by combining it with data
augmentation-trained RoRD [7] which has better rotation coverage but is often
inaccurate, thus leading to an overall state-of-the-art rotation-robust local
feature matching. We make the following contributions:

  1. 1.

We introduce novel Rotation-equivariant Features, dubbed ReF, based on
steerable E2-CNNs which can be group-pooled to achieve rotation-invariance for
local feature matching, thus leading to $0^{\circ}$-like matching performance
under different rotations including $90^{\circ}$ and $180^{\circ}$.

  2. 2.

We present a ‘correspondence ensemble’ of E2-CNN based ReF with data-augmented
standard CNN based RoRD [7], improving the overall coverage of performance
response to different rotation angles.

  3. 3.

We introduce “UrbanScenes3D-Air”, a new dataset that enables benchmarking
Visual Place Recognition from a drone’s view under varying directions of
travel.

  4. 4.

While evaluating local feature matching both standalone and on the VPR task,
we present a detailed analysis of steerable CNNs for rotation-robustness under
the effect of ensembling, robust estimation (RANSAC), network architecture
variations, and the use of rotation priors.

##  II RELATED WORK

Here we briefly discuss local feature matching and rotation equivariant
networks.

###  II-A Local Feature Matching

Local feature matching plays a central role in many computer vision
applications such as visual place recognition, image stitching, pose
estimation, etc. Traditional approaches like SIFT [9], ORB [10] try to address
specific problems like scale invariance, or rotation invariance in the
descriptors. But their performance heavily depends on the rules used for
keypoint detection and description.

Recent works leverage CNNs for keypoint detection and description. D2-Net [8]
uses detect and describe approach which is jointly trained using triplet loss.
SuperPoint [11] uses a self-supervised approach to learn the keypoints and
descriptors. R2D2 [14] uses prediction of reliability map and repeatability
map to find unique features which are repeatable under some transformations.
These approaches perform well under challenging conditions where standard rule
based systems fail.

Some works for learning rotation invariant descriptors include LBP-HF [15,
16], which leverages local binary pattern histograms to learn invariant
features. RI-LBD [17] uses learning based system to learn invariant local
descriptors. RoRD [7] uses data augmentation to learn rotation robust
descriptors. These approaches try to obtain descriptors that are rotation
invariant. Using augmentation during training does improve the performance
compared to standard methodologies [7]. Our approach, on the other hand, tries
to make the model inherently rotation equivariant. This provides a model which
can generalize better over high rotation changes.

##### Visual Place Recognition (VPR)

VPR is commonly used to describe the ability to recognize a previously visited
location despite significant appearance or viewpoint changes [3]. Local
feature matching is also a significant part of VPR pipelines [1]. Recent works
have leveraged advances in deep networks based patch aggregation [2], saliency
detection [18], and transformers with multi-scale attention [19]. As per our
knowledge, our work is the first to utilize a rotation equivariant system for
the task of VPR and demonstrate state-of-the-art results on real-world scenes
from our UrbanScenes3D-Air dataset.

###  II-B Rotation Equivariance Networks

Equivariance is a property where a transformation in the input space leads to
the corresponding transformation change in the feature space. CNNs are known
to be translation equivariant, which means that a translation change in input
space leads to a corresponding translation change in feature space. Rotation
equivariance is a property desired by a lot of vision applications in medical
images, robotics for drones, etc. Steerable CNN is the tool used to achieve
this property [12, 13]. The idea is to learn group representations that
provide equivariance under the actions of the elements of the group. Jenner
provides a framework unifying such convolutions and differential operators in
[20]. Cesa and Weiler derive a specific implementation of these concepts for
the euclidean group $E(2)$ (group of isometries in $E^{2}$ space - rotation,
translation, mirroring) in [13].

The efficacy of equivariant deep learning has been explored in the field of
medical imaging [21], object detection [22], aircraft detection [23] and
reinforcement learning [24]. To the best of our knowledge, this concept has
not been applied to feature description to obtain descriptors that are
rotation equivariant. We demonstrate the use of this concept in local feature
matching by designing a novel architecture that is trained using concepts
inspired from [14]. Our architecture increases the descriptor matching
performance significantly over a wide range of rotational viewpoint changes.

Figure 1: Model and Steerable CNN Implementation - $C_{4}$ equivariant
($n=4$). Left: The input is presented as a _trivial representation_ (RGB
image), and passed through the rotation equivariant layers $\Phi$. Center
_unpooled_ descriptors: The result of the convolution follows the regular
representation (n-bit in-place rotation along the descriptor channel). Center
_pooled_ descriptors: The group pooling is the max substitution of n-bits in
the descriptor (of each feature). Note that the pooled descriptors of features
in images related by pure rotation are _equal_ (therefore rotation
equivariant). The features themselves are rotation invariant (they do not
change based on rotation). Right (output): The feature _descriptors_ are
obtained by taking the $L_{2}$ norm of the pooled descriptors. The
_repeatability_ and _reliability_ maps (keypoints) are obtained through a
standard CNN and a softmax.

##  III METHODOLOGY

This paper aims to develop rotation invariant descriptors for local feature
matching by inducing bias in the model architecture itself, complementary to
learning invariance through data augmentation based training. We achieve this
by using Steerable CNNs instead of standard/conventional CNNs. In this
section, we first introduce the Steerable CNNs, then describe our specific
model architecture and its training methodology, which is then followed by
local feature matching and ensemble strategy.

###  III-A Steerable CNN

Standard CNNs operating on a 2D image grid are translation equivariant.
Steerable CNNs enable us to extend the property of equivariance to groups
other than $(\mathbb{R},+)$ (action group of 2D translations) using group
theory. Cohen et al. introduced Steerable CNNs in [12]. For an input image
$x$, an action group such as a rotation homography $H$ can be applied such
that the output of the function $f$ is steerable (that is, _H-equivariant_) if
the criterion in the equation 1 is met.

$\left(f\circ\pi_{i}(h)\right)(x)=\left(\pi_{o}(h)\circ f\right)(x)\quad h\in
H$ (1)

where the functions $\pi_{i}$ and $\pi_{o}$ are the input and output group
representations. As described in [25], linear representations of finite groups
allow representing group actions as linear morphism, where rotation can be
implemented as a homography group representation. Specifically, we employ the
action group $H=C_{n}$ which consists of rotations
$\\{\theta_{p}\mid\theta_{p}=\frac{2\pi}{n}p\,,\;\;\;p\in\mathbb{Z}/n\mathbb{Z}\\}$.
The group actions are defined on a particular _field type_.

The raw RGB image input to our network uses a ‘trivial representation’ of the
$C_{n}$ group as its field type, which applies no transformation (identity
map) and presents images as RGB images. The field type of output and
intermediate layers uses the ‘regular representation’, which cyclically
rotates the elements, as shown under Group Transform Law in the model
architecture (in Figure 1). In this form, we refer to the output descriptors
as the _unpooled_ descriptors (the $n$ _bits_ in descriptors rotate - the
representation for $n=4$ is shown in Figure 1). Group pooling, which is in-
place substitution by the maximum of the n-bits, converts these unpooled
descriptors to _pooled_ descriptors. This, as shown in Figure 1 (center region
- from unpooled to pooled), makes the descriptors rotation equivariant where
the rotated versions (of feature descriptors) are only related by a rotation
homography (highlighted in the pooled descriptors in center of Figure 1). Note
that the pooled descriptors of specific pixels have now become invariant to
rotation by the angles in $C_{n}$ (the vector doesn’t change, only its
location in the output space changes).

###  III-B Model Architecture

Our architecture, shown in Figure 1, consists of 5 layers of E2-CNN with
increasing channels: (32, 64, 128, 256, 512). Each E2-CNN layer is followed by
a batch normalization layer and ReLU non-linearity. These layers are followed
by a group-pooling layer (see the center of Figure 1) which pools the feature
field across different rotations to get equivariant features. That is, for
each input field, an output field is built by taking the maximum activation
within that field; as a result, the output field transforms according to a
trivial representation. This is required as the computations of standard CNN
needs trivial representation. Moreover, the feature map obtained after group
pooling is equivariant under rotation as well as translation, this is because
E2-CNNs also strides along the x-axis and y-axis as the standard CNNs. We
implement the network for $C_{8}$ space, therefore group pooling reduces the
number of channels from 512 to 64. These pooled features are then passed
through a standard CNN to get a reliability map. Similar operation is done to
get the repeatability map by passing the pooled features through a standard
CNN. Softmax is applied on the reliability map and softplus is applied on the
repeatability map as non-activation. The pooled features are $L2$-normalized
to get the pixel level descriptors that can be used for feature matching. The
model is trained using training methodology of R2D2 [14]. We use R2D2 as the
base training methodology as it enables us to train for feature detection and
description from scratch. Repeatability map is trained jointly using two
losses, one which tries to maximize the cosine similarity between the
repeatability maps of two images, and the other which tries to maximize the
local peakiness of the repeatability map. Reliability map is basically the
discriminative score of the descriptors. Descriptor matching can be seen as a
ranking optimization problem and hence tries to maximize for differentiable
average-precision defined in [26].

###  III-C Correspondence Ensemble

We follow similar ensemble strategy as mentioned in [7], though we use two
separate networks to generate the keypoints, descriptors, and matches. Our
ensemble strategy is based on extracting the top correspondences from the pool
based on the cosine similarity between their descriptors. Once we have
keypoints and descriptors from our model ReF and the RoRD model [7], mutual
nearest-neighbor based feature correspondences are obtained between an image
pair independently for each model. Using the descriptor similarity of the
combined set of correspondences obtained from the two models, we select the
top 50% of the correspondences which have the highest cosine similarity
values. These correspondences are then filtered by RANSAC-based geometric
verification to obtain a final set of feature correspondences. In Section V-A,
we present results for both the settings, that is, without RANSAC (using only
cosine similarity based ensemble) and with RANSAC.

###  III-D Feature Matching

We use mutual nearest neighbor for finding the correspondences between the
image pair once we have the keypoints and corresponding descriptors. We apply
RANSAC on the obtained matches to filter out noisy matches. We can see in
section V-A that our model is able to perform well on different rotation
variations. RANSAC helps us highlight the fact that our model indeed learns
correct correspondences.

Figure 2: The drone traverses along two trajectories of which the yellow
corresponds to the query and the brown correspond to the reference. At a given
location, 8 different images are recorded by the drone during its time of
flight in the brown trajectory. The sample images corresponding to the 8
cameras can be seen here which have a relative degree of rotation with respect
to the yaw axis of the drone.

##  IV Experimental Setup

In this section, we describe the datasets and the benchmarks used to train and
evaluate our contribution.

###  IV-A Datasets

####  IV-A1 Training Data

Similar to R2D2, our training dataset is a combination of synthetic pairs
generated from random web images and Aachen dataset [14]. Note that we do not
use any image from the HPatches evaluation dataset or our ”UrbanScenes3D-Air”
dataset during training.

####  IV-A2 HPatches

We use this dataset for evaluating our model on the task of feature matching.
The dataset contains a total of 116 scenes, 57 of which relate to illumination
changes and the remaining 59 to scene viewpoint variation. Each scene contains
a set of six images and five ground truth homographies between the first and
the rest of the five images. During the evaluation, we resize the image pairs
to (300,300) in dimensions before matching to overcome the GPU OOM issues. We
also transform the ground-truth homographies according to the scaling done on
the images.

Further, we additionally rotate the rest of the five images in a given scene
incrementally in the range of $0^{\circ}$ to $360^{\circ}$ to verify the
feature matching robustness across the rotation angles. This results in
performing matching on a total of $580$ image pairs at every angle and
$13,920$ image pairs in total across all the angles. This rotation is an
additional homography applied to the existing homographies between the image
pairs.

####  IV-A3 “UrbanScenes3D-Air” Dataset for VPR

While several Visual Place Recognition (VPR) datasets exist in the literature
[27, 28, 3]111[29] is a relevant dataset for this research but not yet
released., the extent of viewpoint variation characteristics they offer is
often either limited or tied to particular applications such as on-road
driving [30]. To complement the existing benchmark, we propose
UrbanScenes3D-Air dataset for UAV/drone based VPR. For this purpose, we
imported the real-world Urbanscenes 3D dataset [31] into the AirSim
environment [32]. We then created a virtual drone equipped with downward-
facing RGB camera and flew it over the buildings in the that environment to
capture image sequences. To generate our UrbanScenes3D-Air dataset, the drone
was flown in two trajectories as shown in the Figure 2. We pick the yellow
trajectory for populating query database with 20 images taken at different
locations. For this query image trajectory, we attach a single camera and fly
the drone heading straight along the yellow path and capture the downward
camera’s views at regular intervals. Further, the reference database is
populated with about 200 images from the brown trajectory for a given query
image. While recording this brown trajectory we place 8 cameras at the bottom
of the drone, each of which has a relative rotation in yaw. In Figure 2, the 8
blue lines represent the setup of the cameras, each having a relative rotation
with respect to drone’s yaw axis. So, at a given location, we are able to
capture 8 images which are relatively rotated in yaw. We do this to simulate
other trajectories that might pass through the original trajectory at an angle
but if done so would not have much trajectory overlap, thus this lets us test
the same thing but with more image pairs in reference database at
intersection. We formulate the problem of VPR here as retrieving the correct
image match for a given query image from a set of images that have a relative
rotation to the query. This is similar to the drone crossing a previous
trajectory at an angle and re-localizing itself at the intersection of the
previous and the current trajectory. The relative angle of rotation in our
dataset varies from 0 to 90 and also includes opposite viewpoint images with
180 degree relative rotation. Figure 2 depicts samples of query images which
have a relative rotation with the the first image which is a sample of the
reference database. We further present our evaluation process of the VPR on
our UrbanScenes-Air dataset in IV-B2

|
---|---
|
|
Figure 3: The first row compares the performance of various methods on
HPatches at different angles with and without RANSAC. The second row shows the
individual performance of only illumination+rotation variations and
viewpoint+rotation variations. ReF (Ours) outperforms all the other methods in
comparison at certain angles in the neighbourhood of discrete C8 space like
$75^{\circ}$, $90^{\circ}$ and $105^{\circ}$ when we apply RANSAC based
filtering on the matches obtained from mutual nearest neighbours. Further, our
proposed ensembling of ReF with RoRD improves the overall performance across
rotation angles. All the scores are calculated for error threshold of 3
pixels.

###  IV-B Tasks and Evaluation

####  IV-B1 Local Feature Matching

We evaluate the local feature matching on HPatches using Mean Matching
Accuracy (MMA). A feature correspondence is considered to be correct if the
back-projection of the keypoint using the ground truth homography matrix falls
within a specific pixel threshold. We use pixel threshold of 3 for plots
showing performance at various angles. Figure 4 shows performance averaged
across all angles for different pixel thresholds. We also compare the
performance after filtering the correspondences with RANSAC since it is
typically used to filter matches for tasks like VPR [2, 7] and 6-DoF pose
estimation [5].

####  IV-B2 Visual Place Recognition

We structured our UrbanScenes-Air dataset so that the ideal match in the
reference database has the same `image_index` as the query image. After
feature matching, our VPR pipeline returns the image in the reference database
with the most inliers for a given query image. For evaluation, the predicted
match is considered to be correct when its `image_index` lies within a
threshold of $\pm 2$ indices from the query image’s index. This is then
repeated by changing the query database images to other angles like 15, 30,
45, 60, 75, 90, and 180. For each generated query database we evaluate and
compare our ReF and the proposed ensemble of ReF+RoRD pipelines with other
state-of-the-art feature matching pipelines. We plot the recall scores as a
metric for VPR evaluation at various angles on our UrbanScenes3D-Air dataset
in Figure 5.

##  V RESULTS

In this section, we first present benchmark results for standalone local
feature matching, then present both quantitative and qualitative results for
VPR, which is finally followed by additional analyses of rotation robustness
of E2-CNNs.

###  V-A Local Feature Matching

####  V-A1 Benchmarking across rotation angles

The quantitative results for local feature matching on HPatches dataset is
presented in Figure 3. This evaluation is performed on a rotated variant of
HPatches as proposed in [7] and described in IV-A2 . Figure 3 presents MMA
scores for both with and without RANSAC (top row) along with split results on
illumination and viewpoint alone (bottom row). It can be observed that ReF
(blue) not just achieves peak performance at certain angles from its discrete
C8 space but also in its immediate neighborhood, for example, $75^{\circ}$ and
$105^{\circ}$ in the neighborhood of $90^{\circ}$. On the other hand, all
methods including RoRD [7] and SuperPoint [11], which are trained with heavy
data augmentation for rotation robustness, exhibit meaningful MMA only under
moderate rotation variations, for example, up to $\pm 60^{\circ}$. RoRD,
however, performs the best among all existing data-augmented CNNs, and is thus
used in our ensemble model (green) which improves the overall performance
coverage across rotation angles with slight reduction of ReF’s peak
performance. This trend is further highlighted in Table I which shows average
MMA across all angles when using RANSAC.

Figure 4: MMA values averaged across a range of angles in [$0^{\circ}$,
$360^{\circ}$], at various pixel thresholds, for multiple methods. Ensemble
outperforms other methods at higher pixel thresholds, where as ReF performs
best at lower pixel thresholds.

Figure 5: VPR on UrbanScenes3D-Air: The ensemble of ReF with RoRD outperforms all the other methods in comparison. Method | Recall(%) | MMA(3px) | MMA(6px)
---|---|---|---
R2D2 | 28.75 | 0.138 | 0.156
D2-Net | 46.88 | 0.0826 | 0.188
SuperPoint | 41.25 | 0.214 | 0.242
RoRD | 77.50 | 0.212 | 0.576
ReF (ours) | 61.25 | 0.469 | 0.499
ReF (ours) + RoRD | 86.88 | 0.393 | 0.669
TABLE I: Recall and MMA scores averaged across rotation angles as evaluated on
UrbanScenes3D-Air dataset for VPR and local feature matching on HPatches
respectively. The average MMA values reported here were computed at pixel
threshold 3 and 6 respectively in the last two columns, corresponding to the
top left plot in Figure 3 using RANSAC based filtering.

##### With and Without RANSAC

In Figure 3 top row, we compare results with and without RANSAC. It can be
seen that while there is general improvement of MMA scores through RANSAC for
all the methods, some methods benefit more than others, e.g. ReF (blue). This
indicates that the feature learnt by ReF can lead to correct correspondences
but with noisy outliers, which can then be successfully filtered by a robust
outlier rejection method like RANSAC, thus validating ReF’s practical use in
tasks like VPR. Other methods do not show a similar behaviour since their
initial correspondences set does not comprise sufficient correct matches. On
the other hand, the ensemble method (green) is able to leverage the correct
correspondences from both RoRD and ReF, where its performance in most cases is
close to the better performing ingredient method (RoRD or ReF).

##### Illumination vs Viewpoint

In Figure 3 bottom row, we compare performance split across illumination and
viewpoint variations (with RANSAC). Here, the viewpoint variations refer to
the perspective viewpoint change originally presented in the HPatches dataset.
It can be observed that performance trends are maintained across the split
results but the combined effect of perspectivity change on top of rotation
variations reduces the absolute MMA.

####  V-A2 Benchmarking at different pixel thresholds

In the previous study, we considered the pixel threshold of 3. Here, we
observe the MMA scores averaged across all the rotation variations plotted
against different pixel thresholds in Figure 4. The ensemble method (solid
green) outperforms all the other methods across the board except the tighter
thresholds where ReF (solid blue) performs the best. It can also be observed
that use of RANSAC (solid) changes the performance trends as compared to
without RANSAC (dashed). This justifies the need for RANSAC to correctly
select the right set of correspondences from an initial set of matches which
is far less noisy than that obtained from the other state-of-the-art methods.

###  V-B VPR on UrbanScenes3D-Air

The performance of various feature matching pipelines on our UrbanScenes3D-Air
dataset for the task of VPR can be seen in Figure 5. As the relative angle of
rotation increases from 0 to 90 degrees, ReF(Ours) performs significantly well
at discrete higher rotations when compared to the other state-of-the-art
methods. Further, the ensemble of our ReF with RoRD improves the recall scores
at specific angles which are not handled by ReF and thus outperforms all the
other methods in comparison. The same trend is also summarized in Table I. The
qualitative matches shown in Figure LABEL:fig:vpr_qualitative highlight the
superior performance of ReF model under large rotation variations with rich
and accurate correspondences where methods like R2D2 [14] and SuperPoint [11]
fail.

Figure 6: We can clearly see that our unpooled version performs better than
the group pooled version when a rotation prior is given at angles where the
group pooled version fails like $45^{\circ}$, $135^{\circ}$, etc.

###  V-C Ablation Studies

####  V-C1 Performance with Rotation Priors

The E2-CNN layers are equivariant under discrete rotations and we are able to
extract feature map specific to a specific angle. So instead of group pooling
the feature fields obtained from E2-CNN in our network, we directly use the
unpooled features. Here, we evaluate the performance in a scenario where we
have a rotation prior available. That is, we know the relative orientation
between the two images for matching, which in some practical use cases could
be known, for example, back and forth last mile deliveries on repeated routes.
Figure 6 shows the MMA results for this setting. It can be seen that the
features learned by C8 space are better highlighted in this scenario where
there is a performance jump at $45^{\circ}$ as well.

Figure 7: Equivariance is lost when standard convolutions are added after the
group pooling layer. Figure 8: Performance comparison of various cyclic groups
$C_{N}$ for $N\in\\{4,8,16\\}$.

####  V-C2 Effect of mixing Standard CNN with E2-CNN

In this section we evaluate the impact of adding standard CNN layers after
group pooling layer for the descriptors. A lot of precious work on
classification [13, 12] have shown various model ablations which combine
different types of layers. For this experiment we add 3 standard CNN layers
with filter size 2 with stride increment of 4,8 and 16 after group pooling
layer. Figure 7 shows that the equivariance property is completely lost, and
the performance is flattened out across the angles, whereas in Figure 3 the
performance of ReF peaks at specific angles. So, directly using the group
pooled features as descriptors is necessary for local feature training to
work. Although our ReF-network has 2 standard CNN layers to the end, those are
only for classification and loss computation, and do not operate on the
predicted descriptors after group pooling, thus preserving the invariance.

####  V-C3 Effect of using different cyclic groups

In Figure 8, we compare the performance of various cyclic groups $C_{N}$ for
$N\in\\{4,8,16\\}$. It can be observed that performance does not change
noticeably across different values. Our choice of $C_{8}$ was based on a
balance between performance and compute given that in the unpooled settings
(that is with a given rotation prior, see Figure 6), one could still achieve
performance response on angles like $45^{\circ}$ which is not possible by
definition for $C_{4}$ space.

##  VI CONCLUSIONS

In this work we present a novel architecture which leverages a steerable CNN
for local feature matching between images under rotation variations. We
further extend the capability of our model by introducing an ensemble approach
combining ReF + RoRD, which combines the benefits of both models and improves
overall performance across all possible rotations. We show that this ensemble
outperforms other methods on the UrbanScenes3D-Air dataset, and on HPatches at
higher pixel thresholds.

We explore the possibility for using steerable CNN for feature matching. This
work can further be extended by studying the effects of different network
architecture nuances like residual connections and combining steerable CNN on
different spaces, e.g. C16, C8, C4, in the same network. Our study shows that
by naively combining steerable CNN with standard CNN the equivariance property
is lost. But more rigorous study on this front could lead to more novel
architectures which help in feature matching under various rotations and
extreme viewpoint changes.

## References

  * [1] M. Cummins and P. Newman, “Appearance-only slam at large scale with fab-map 2.0,” _The International Journal of Robotics Research_ , vol. 30, no. 9, pp. 1100–1123, 2011.
  * [2] S. Hausler, S. Garg, M. Xu, M. Milford, and T. Fischer, “Patch-netvlad: Multi-scale fusion of locally-global descriptors for place recognition,” in _Proceedings of the IEEE/CVF Conference on Computer Vision and Pattern Recognition_ , 2021, pp. 14 141–14 152.
  * [3] S. Garg, T. Fischer, and M. Milford, “Where is your place, visual place recognition?” in _Proceedings of the Thirtieth International Joint Conference on Artificial Intelligence (IJCAI-21)_. International Joint Conferences on Artificial Intelligence, 2021, pp. 4416–4425.
  * [4] T. Liao and N. Li, “Single-perspective warps in natural image stitching,” _IEEE Transactions on Image Processing_ , vol. 29, p. 724–735, 2020. [Online]. Available: http://dx.doi.org/10.1109/TIP.2019.2934344
  * [5] T. Sattler, W. Maddern, C. Toft, A. Torii, L. Hammarstrand, E. Stenborg, D. Safari, M. Okutomi, M. Pollefeys, J. Sivic, _et al._ , “Benchmarking 6dof outdoor visual localization in changing conditions,” in _Proceedings of the IEEE conference on computer vision and pattern recognition_ , 2018, pp. 8601–8610.
  * [6] P.-E. Sarlin, D. DeTone, T. Malisiewicz, and A. Rabinovich, “SuperGlue: Learning feature matching with graph neural networks,” in _CVPR_ , 2020. [Online]. Available: https://arxiv.org/abs/1911.11763
  * [7] U. S. Parihar, A. Gujarathi, K. Mehta, S. Tourani, S. Garg, M. Milford, and K. M. Krishna, “Rord: Rotation-robust descriptors and orthographic views for local feature matching,” 2021.
  * [8] M. Dusmanu, I. Rocco, T. Pajdla, M. Pollefeys, J. Sivic, A. Torii, and T. Sattler, “D2-Net: A Trainable CNN for Joint Detection and Description of Local Features,” in _Proceedings of the 2019 IEEE/CVF Conference on Computer Vision and Pattern Recognition_ , 2019.
  * [9] D. Lowe, “Distinctive image features from scale-invariant keypoints,” _International Journal of Computer Vision_ , vol. 60, pp. 91–, 11 2004.
  * [10] E. Rublee, V. Rabaud, K. Konolige, and G. Bradski, “Orb: an efficient alternative to sift or surf,” 11 2011, pp. 2564–2571.
  * [11] D. DeTone, T. Malisiewicz, and A. Rabinovich, “Superpoint: Self-supervised interest point detection and description,” _CoRR_ , vol. abs/1712.07629, 2017\. [Online]. Available: http://arxiv.org/abs/1712.07629
  * [12] T. S. Cohen and M. Welling, “Steerable cnns,” _CoRR_ , vol. abs/1612.08498, 2016. [Online]. Available: http://arxiv.org/abs/1612.08498
  * [13] M. Weiler and G. Cesa, “General E(2)-Equivariant Steerable CNNs,” in _Conference on Neural Information Processing Systems (NeurIPS)_ , 2019.
  * [14] J. Revaud, P. Weinzaepfel, C. D. Souza, N. Pion, G. Csurka, Y. Cabon, and M. Humenberger, “R2d2: Repeatable and reliable detector and descriptor,” 2019\.
  * [15] T. Ahonen, J. Matas, C. He, and M. Pietikäinen, “Rotation invariant image description with local binary pattern histogram fourier features,” 06 2009, pp. 61–70.
  * [16] T. Ojala, M. Pietikainen, and T. Maenpaa, “Multiresolution gray-scale and rotation invariant texture classification with local binary patterns,” _IEEE Transactions on Pattern Analysis and Machine Intelligence_ , vol. 24, no. 7, pp. 971–987, 2002.
  * [17] Y. Duan, J. Lu, J. Feng, and J. Zhou, “Learning rotation-invariant local binary descriptor,” _IEEE Transactions on Image Processing_ , vol. 26, no. 8, pp. 3636–3651, 2017.
  * [18] N. V. Keetha, M. Milford, and S. Garg, “A hierarchical dual model of environment-and place-specific utility for visual place recognition,” _IEEE Robotics and Automation Letters_ , vol. 6, no. 4, pp. 6969–6976, 2021\.
  * [19] R. Wang, Y. Shen, W. Zuo, S. Zhou, and N. Zhen, “Transvpr: Transformer-based place recognition with multi-level attention aggregation,” 2022.
  * [20] E. Jenner and M. Weiler, “Steerable partial differential operators for equivariant neural networks,” _CoRR_ , vol. abs/2106.10163, 2021. [Online]. Available: https://arxiv.org/abs/2106.10163
  * [21] P. Müller, V. Golkov, V. Tomassini, and D. Cremers, “Rotation-equivariant deep learning for diffusion MRI,” _CoRR_ , vol. abs/2102.06942, 2021. [Online]. Available: https://arxiv.org/abs/2102.06942
  * [22] J. Han, J. Ding, N. Xue, and G.-S. Xia, “Redet: A rotation-equivariant detector for aerial object detection,” 2021.
  * [23] X. Chen, J. Liu, F. Xu, Z. Xie, Y. Zuo, and L. Cao, “A novel method of aircraft detection under complex background based on circular intensity filter and rotation invariant feature,” _Sensors_ , vol. 22, no. 1, 2022\. [Online]. Available: https://www.mdpi.com/1424-8220/22/1/319
  * [24] A. K. Mondal, P. Nair, and K. Siddiqi, “Group equivariant deep reinforcement learning,” _arXiv preprint arXiv:2007.03437_ , 2020.
  * [25] J.-P. Serre, _Linear Representations of Finite Groups_ , 1st ed. Springer, 1977.
  * [26] K. He, Y. Lu, and S. Sclaroff, “Local descriptors optimized for average precision,” 2018.
  * [27] M. Zaffar, S. Garg, M. Milford, J. Kooij, D. Flynn, K. McDonald-Maier, and S. Ehsan, “Vpr-bench: An open-source visual place recognition evaluation framework with quantifiable viewpoint and appearance change,” _International Journal of Computer Vision_ , vol. 129, no. 7, pp. 2136–2174, 2021.
  * [28] F. Warburg, S. Hauberg, M. Lopez-Antequera, P. Gargallo, Y. Kuang, and J. Civera, “Mapillary street-level sequences: A dataset for lifelong place recognition,” in _Proceedings of the IEEE/CVF conference on computer vision and pattern recognition_ , 2020, pp. 2626–2635.
  * [29] A. Vallone, F. Warburg, H. Hansen, S. Hauberg, and J. Civera, “Danish airs and grounds: A dataset for aerial-to-street-level place recognition and localization,” _arXiv preprint arXiv:2202.01821_ , 2022.
  * [30] S. Garg, N. Suenderhauf, and M. Milford, “Lost? appearance-invariant place recognition for opposite viewpoints using visual semantics,” _Proceedings of Robotics: Science and Systems XIV_ , 2018.
  * [31] Y. Liu, F. Xue, and H. Huang, “Urbanscene3d: A large scale urban scene dataset and simulator,” 2021.
  * [32] S. Shah, D. Dey, C. Lovett, and A. Kapoor, “Airsim: High-fidelity visual and physical simulation for autonomous vehicles,” in _Field and Service Robotics_ , 2017. [Online]. Available: https://arxiv.org/abs/1705.05065


# Binary fractions of filed F, G, and K dwarf stars based on the Gaia Early
Data Release 3 and LAMOST Data Release 5

Zexi Niu  National Astronomical Observatories, Chinese Academy of Sciences
20A Datun Road, Chaoyang District, Beijing, China  Haibo Yuan  Department of
Astronomy, Beijing Normal University
19th Xinjiekouwai Street, Haidian District, Beijing, China  Song Wang
National Astronomical Observatories, Chinese Academy of Sciences
20A Datun Road, Chaoyang District, Beijing, China  Jifeng Liu  National
Astronomical Observatories, Chinese Academy of Sciences
20A Datun Road, Chaoyang District, Beijing, China

###### Abstract

Basing on the large volume Gaia and LAMOST data, we are able to estimate the
bias-corrected binary fraction and explore its tendencies with different
stellar populations. A 0.2 million volume-limited sample of the filed late F,
G, and early K dwarfs are well selected. Thanks to the inclusive stellar locus
outlier method, both close and wide binaries are involved in this work. The
general binary fractions decrease toward lower mass and are 44$\%\pm$ 1$\%$
for the solar mass range sample and 42$\%\pm 1\%$ for a lower mass range. A
positive correlation between [$\alpha$/Fe] and the binary fraction is found in
thin disk within a mono-metallicity and mono-mass sample. This tendency
disappears in the thick disk. We conjecture that it is related to the stellar
ages by comparing the age distributions in the disks. Combing with the
previous work, we suggest that the correlation between [$\alpha$/Fe] and the
binary fraction changes with the varying periods, just like the correlation
between [Fe/H] and the binary fraction. We also give the averaged binary
fractions of the thin disk, thick disk, and inner halo. Disk stars are found
to have similar binary fractions. Inner halo stars have a great possibility to
own the highest fraction.

binaries:general - stars:abundances - stars:solar type

††journal: ApJ

##  1 Introduction

Binary systems are ubiquitous in our universe. Understanding star formation
and evolution requires a comprehensive knowledge of binary systems. Being
progenitors and hosts of various interesting systems, they are important in
many fields, such as Type Ia supernovae (e.g. 2000iasn), short-duration Gamma-
Ray Bursts (e.g. 1992grb), kilonova (e.g. 2017kilonova), cataclysmic variables
(e.g. 1995cv), X-ray binaries, stellar population synthesis (e.g. 2002sps),
and gravitational wave astronomy (e.g. 2017gw170817).

In the past few decades, the statistics of binary systems has been extensively
studied with various observational methods: visual techniques (1989visual;
1996visual), radial velocity (RV) variation (gao2014; gao2017; tian2018;
2018badnes), spectral fitting (2007wdms; 2018specbi; 2020pyhammer; 2021xiang),
eclipsing binary (2011keplereb1; 2011keplereb2; 2019tesseb), common proper
motion (2004cpm; 2020Hartmancpm; widedr2), displacement in the color-magnitude
diagram (liuchao2019; 2012Milone), astrometric noise excess (2019Kervella;
2020Penoyre), and color locus outlier (2015loci2). Taking advantage of
multiple methods and results above, some review works summarize and present
census of binary systems in the solar neighborhood such as raghavan2010,
moe2017, and moe2019anti. The positive correlation between the binary fraction
and stellar mass is well established as well as the anti-correlation between
the close binary fraction and metallicity.

However, works focusing on the variations of the binary fraction with
[$\alpha$/Fe] and different stellar populations are rare. The reason includes
not only the relatively small data volume but also the bias of the individual
methods in period, mass ratio, and distance limitations, etc (moe2017). To
overcome the biases and limitations of the previous techniques, 2015loci2
proposed a stellar locus outlier (SLOT) method to estimate model-free binary
fraction for large numbers of stars of different populations in large survey
volumes. Applying the method to about 10,000 stars from the Sloan Digital Sky
Survey (SDSS) Stripe 82, 2015loci2 have determined the binary fractions of
different stellar populations and analyzed their dependence on spectral type
and [Fe/H]. They find the highest binary fraction in the Galactic halo and
comparable values in the thin and thick disks. The result is consistent with
moe2019anti, but against those from 1983Carney and 2002Latham. For $\alpha$
abundances, tian2018 suggest a positive correlation between [$\alpha$/Fe] and
binary fraction for solar-type stars in the thin disk. While aopgee2020alpha
perform a detailed study on the close binary fraction as a function of stellar
parameters including $\alpha$ abundances and conclude that $\alpha$ elements
suppress multiplicity at most values of [Fe/H]. More efforts are needed to
study the variations of the binary fraction with [$\alpha$/Fe] and different
stellar populations.

Thanks to the unprecedented data volume and quality of Gaia (gaia2016) and
LAMOST (zhao2012) surveys, we can apply the SLOT method to a more advanced
sample in terms of larger sample size, better data quality, and more stellar
populations. Note that previous works have revealed that Gaia photometry
suffers a magnitude-dependent systematic errors up to 10 mmag (maw2018;
casa2018). Fortunately, niudr2; niudr3 present precise corrections on colors
using the spectroscopy-based stellar color regression method (2015scr),
achieving a precision of about one mmag. Therefore, the LAMOST-Gaia stars
serve as an adequate sample to explore the binary fractions with various
factors, especially [$\alpha$/Fe] and stellar populations, which is the focus
of this paper.

The paper is organized as follows. Section 2 describes the data selection,
mainly focusing on late F to early K dwarfs in the solar vicinity. The SLOT
method is described in Section 3. Section 4 presents the results for the
individual samples and explores the impacts of stellar parameters on binary
factions. The results are discussed in Section 5 and compared with previous
works. Section 6 concludes by summarizing the main findings and discussing
future work.

##  2 Data

###  2.1 Data selection

In this work, we use the Gaia Early Data Release 3 (EDR3; gaia2020) and the
fifth data release of LAMOST (LAMOST DR5; zhao2012; luo2012), accompanying
with the $\alpha$ abundances provided by xiangdd. Gaia EDR3 provides 1.8
billion photometric data of $G$ band and 1.5 billion of $G_{\rm BP}$ and
$G_{\rm RP}$ bands based on its 36 months’ observation. LAMOST database has
accumulated more than 8 million stellar spectra in DR5 with the spectral
resolution R $\sim$ 1800\. Effective temperature $T_{\rm eff}$, surface
gravity log $g$, and metallicity [Fe/H] are delivered from the spectra through
the LAMOST Stellar Parameter pipeline (LASP; wu2011) with the precision of
about 110 K, 0.2 dex, and 0.1 dex, respectively (luo2015). xiangdd uses a
data-driven Payne approach to provide abundances of 6 million stars for 16
elements as well as the [$\alpha$/Fe], which is defined as a weighted mean of
[Mg/Fe], [Si/Fe], [Ca/Fe], and [Ti/Fe]. The results are publicly available as
a value-added catalog of the LAMOST DR5.

As proposed by gaia2020 and phot_cont, the corrected $G$ magnitudes for
sources with 6-parameter astrometric solutions are used in this paper. Further
color and magnitude corrections are applied. The former is from niudr3 as
introduced above. The latter is provided by yl, where a machine learning
technique is applied to training the observed $UBVRI$ magnitudes of about
10,000 Landolt standard stars into the Gaia EDR3 magnitudes. For multiply
observed Gaia objects in LAMOST, we combine their epoch $T_{\rm eff}$, log
$g$, [Fe/H] and [$\alpha$/Fe] values using inverse variance weighting. We
select the main sequence (MS) stars in the $T_{\rm eff}$ and log $g$ diagram
using the same criterion in niudr2. Main sequence turn-off (MSTO) stars are
also eliminated in the Hertzsprung-Russell diagram as shown in Figure 1, as
their intrinsic magnitudes can not be precisely determined by their colors and
metallicities.

Figure 1: Selections of the dwarf stars for different metallicity bins. Blue
points are selected. Metallicities are labelled in each panel.

To grantee the quality of the sample, the following constraints are required:

For Gaia data:

  * 1)

phot$\\_$bp$\\_$mean$\\_$flux$\\_$over$\\_$error $>$ 100

  * 2)

phot$\\_$g$\\_$mean$\\_$flux$\\_$over$\\_$error $>$ 100

  * 3)

phot$\\_$rp$\\_$mean$\\_$flux$\\_$over$\\_$error $>$ 100

  * 4)

phot$\\_$proc$\\_$mode = 0

  * 5)

duplicated$\\_$source = False

For LAMOST data:

  * 6)

the signal-to-noise ratios for the $g$ band (S/Ng) are larger than 20

  * 7)

$T_{\rm eff}>$ 4500 K

For spatial location:

  * 8)

$E(B-V)$ $<$ 0.05 mag according to the sfd (sfd, hereafter SFD) dust reddening
map

  * 9)

galactic latitude $|b|>20$ deg

  * 10)

vertical distance to the galactic disk $|Z|>0.2$ kpc

###  2.2 Sample properties

Stellar mass is estimated by interpolating the color mass relation provided by
the PARSEC (parsecbre; parsectang) model at every 0.1 dex from $-$1.5 to 0.5
dex of [Fe/H], as shown in Figure A17 (see Appendix). Since we have eliminated
MSTO stars, we can have a reliable mass estimation. Our sample mainly contains
late F, G, and early K type stars. To reduce the mass effect on binary
fractions, we divide our sample stars into two groups of 0.55 $<m_{1}<$ 0.8M⊙
(hereafter LM) and 0.8 $<m_{1}<$ 1.05M⊙ (hereafter HM). To avoid Malmquist
bias, i.e., multiple systems can be detected further away than single stars in
a magnitude-limited sample, volume-limited samples are preferred in this work.
Distances are limited to 0.95 kpc for the LM sample and 1.5 kpc for the HM
sample, which lets the cut-off of magnitudes of different mass ranges are
similar of about 17 mag.

Because we have specified the spatial location, we assume that all sources are
beyond the source of reddening. The empirically determined reddening
coefficient $R(G-G_{\rm RP})$ and $R(G_{\rm BP}-G_{\rm RP})$ are dependent on
the temperature and reddening, as provided in niudr3. All colors referred to
hereafter are dereddened using the SFD dust map and the empirical
coefficients.

Figure 2: top: Distributions of distance $d$ and the $G$ magnitude for 0.55
$<m_{1}<$ 0.8M⊙ sample (left) and 0.8 $<m_{1}<$ 1.05M⊙ sample (right). middle:
Distributions of the [Fe/H] and [$\alpha$/Fe] in the same order. bottom:
Distributions of radial Galactocentric distances and vertical distance from
the disk in the same order. Color bars indicate numbers.

The absolute Gaia magnitude in the $G$ band is calculated by $M_{G}=G+5+5{\rm
log_{10}}(\varpi/1000)-A_{G}$, where $\varpi$ is the Gaia parallax and
$A_{G}=1.890$ E$(G_{\rm BP}-G_{\rm RP})$ according to the ws.

Finally, our sample contains 200,827 objects, 86,141 of the LM sample and
114,686 of the HM sample. The distributions of stellar parameters of the HM
and LM samples are demonstrated in Figure 2. From top to the bottom, distance
$d$ against the $G$ magnitude, [$\alpha$/Fe] against [Fe/H], and radial
Galactocentric distances $R_{\rm Gal}$ against vertical distance from the disk
$|Z|$. $G$ magnitudes are almost brighter than 17 mag. The discontinuity of
the magnitude distribution around $G\sim 14$ mag is due to the selection
functions of the LAMOST and has no impact on this work. Our sample covers wide
ranges of [Fe/H] of about 2 dex variations and [$\alpha$/Fe] of 0.5 dex
variations. Thin and thick disks can be clearly distinguished in the abundance
plot for both HM and LM samples. Relatively higher [$\alpha$/Fe] and lower
[Fe/H] group at the top left corner of the plot is the thick disk. The other
group at the bottom right is the thin disk.

##  3 Method

###  3.1 SLOT method overview

Figure 3: Blue histograms are the distributions of $G_{\rm BP}-G$ residuals at
different $G$ magnitudes. Red dashed lines are the results of Gaussian model
fitting.

Proposed by 2015loci2, the SLOT method is developed to provide model-free
estimations of binary factions for large numbers of stars of different
populations in large survey volumes. They have applied the SLOT method to two
samples of Stripe 82 stars by combing the recalibrated SDSS photometric data
with the spectroscopic information from the SDSS and LAMOST surveys. After
fitting the $u-g$, $g-r$, $r-i$, and $i-z$ colors as a function of the $g-i$
color and [Fe/H], they find that the fitting residuals are asymmetric,
pointing to the presence of a significant population of binaries, which is
explained vividly by Figure 1 of their paper. This method works well for not
only close binaries but unresolved wide binaries as well. It is also
insensitive to the assumed mass ratio distribution and does not require
multiple epoch data.

In niudr3, we have fitted the metallicity-dependent stellar loci for the Gaia
colors (see Table 1) and detected the asymmetric $G_{\rm BP}-G$ residuals
(defined as loci predicted colors $-$ observed colors, hereafter the same)
after correcting the magnitude-dependent systematic effects of colors. To be
more clearly, Figure 3 demonstrates the distributions of the $G_{\rm BP}-G$
residuals at different $G$ magnitudes from 11.5 to 16.5 mag. The observed
$G_{\rm BP}-G$ residuals, as shown in the blue histograms, are fitted with
Gaussian models. The results are over-plotted in red dashed lines. There are
significant discrepancies between the symmetric Gaussian distributions and the
observed $G_{\rm BP}-G$ residuals in the first 5 panels. The discrepancies
become weaker in the last panel when $G\sim 16.5$ mag, due to the larger
photometric errors at fainter magnitudes. The sensitivity of the SLOT method
might be challenged at faint magnitudes. Fortunately, almost all the sample
stars in our study are brighter than 16.5 mag in the $G$ band, they does not
suffer from this problem.

The SLOT method estimates binary fraction of a given sample by comparing the
distribution of observed color residuals with the simulated ones. Two sets of
Monte Carlo simulations are performed, assuming that all stars in the sample
are either single or binaries. Then by adjusting the relative fraction of
stars in the two sets of simulations to fit the observed residual
distributions, the binary fraction of the sample is determined. However, the
center of the distribution of the observed $G_{\rm BP}-G$ residuals is not
strictly zero (shifting less than 1 mmag from Figure 3), due to the existence
of binary stars in the sample, while the simulations are centered in zero. To
eliminate this problem on estimating of binary fraction, both offset $\mu$ of
the center of the distribution of the observed $G_{\rm BP}-G$ residuals and
binary fraction $f_{\rm b}$ are set to be free parameters. A minimum
$\chi^{2}$ technique is adopted in the fitting:

$\displaystyle\chi^{2}=\sum_{i=1}^{N_{\rm bin}}\frac{(N_{\rm obs}^{i}-N_{\rm
sim}^{i})^{2}}{{N_{\rm obs}^{i}}\times(N_{\rm bin}-1)}$ (1)

$\displaystyle N_{\rm sim}^{i}=N_{\rm binary}^{i}\times f_{\rm b}+N_{\rm
single}^{i}\times(1-f_{\rm b})$ (2)

Considering the maximum effects of the binary influences on the $G_{\rm BP}-G$
residuals and various uncertainties (see the next subsection), when
calculating the $\chi^{2}$, $G_{\rm BP}-G$ residuals are equally binned at the
range from $-$7 to 4 mmag with bin size of 0.275 mmag, because the simulated
residuals barely reach values beyond boundaries and stars beyond boundaries
would contaminate the $\chi^{2}$ calculations. After binning, $N_{\rm
obs}^{i}$, $N_{\rm single}^{i}$, and $N_{\rm binary}^{i}$ are numbers of the
observed, simulated single, and simulated binary residuals in the $i$th bin.
The single and binary sets of simulations are carried out 100 times. For a
given binary fraction $f_{\rm b}$ and an offset $\mu$, the residuals of the
single and binary sets are combined according to Equation 2, and then shifted
toward the bluer side by the offset $\mu$. $N_{\rm sim}^{i}$ is the number of
the final simulated $G_{\rm BP}-G$ residuals divided by 100 in the $i$th bin.
$N_{\rm bin}=40$ is the bin number of the histogram. We vary $f_{\rm b}$ from
0.0 to 1.0 at 0.01 intervals and $\mu$ from 0.0 to 1 mmag at 0.01 mmag
intervals. $\chi^{2}$ is calculated for each set of $f_{\rm b}$ and $\mu$. By
choosing the global minimum $\chi^{2}$, the best solution is determined.

In order to explore the binary fraction of the filed star as a function of
[Fe/H] and [$\alpha$/Fe], the HM and LM samples have been divided into 27 and
35 bins. Bin width for [Fe/H] is fixed to be 0.1 dex, and is flexible for
[$\alpha$/Fe] depending on the number density.

We use Bootstrap Sample method to evaluate formal error of the binary
fraction. One bin with the largest number containing 11,508 objects in the HM
sample is performed random samplings with replacement for 1,500 times. We fit
the distribution of 1,500 $f_{b}$ values by a Gaussian model and get
$\sigma$=0.03 as the formal error of the binary fraction of that bin. Errors
of other bins are then estimated assuming they are inversely proportional to
the square of the ratio of their numbers to 11,508.

###  3.2 Descriptions of the simulation

Here we describe simulations of the binary and single sets step by step. A
schematic description is shown in Figure 4. For each object, the observed
$G_{\rm BP}-G_{\rm RP}$, [Fe/H], and absolute $G$ magnitude are the inputs.
Error in this work is assigned and added through a normal distribution. We use
superscript ′ to indicate values after adding errors. For example, [Fe/H]′ is
sampled by adding the random error $\sigma_{\rm ran}$([Fe/H]) and systematic
error $\sigma_{\rm sys}$([Fe/H]) to the [Fe/H].

When simulating the single set, steps 1 $\sim$ 6 are as below:

Step1: Given the intrinsic $G_{\rm BP}-G_{\rm RP}$ and [Fe/H]′, the intrinsic
$G_{\rm BP}-G$ is computed by the metallicity-dependent loci in Table 1.

Step2: Combining the absolute $G$ magnitude and the intrinsic colors, the
absolute $G_{\rm BP}$ and $G_{\rm RP}$ magnitudes are obtained.

Step3: Random errors in the three magnitudes are added.

Step4: Then $G_{\rm BP}^{\prime}-G^{\prime}$ and $G_{\rm BP}^{\prime}-G_{\rm
RP}^{\prime}$ are derived.

Step5: When yielding the simulated $(G_{\rm BP}-G)^{\prime}$, the effect of
calibration error on the $G_{\rm BP}-G$ residuals is acquired according to the
$G$ magnitude and added to the $G_{\rm BP}^{\prime}-G^{\prime}$. $(G_{\rm
BP}-G)_{\rm loci}^{\prime}$ is computed using $G_{\rm BP}^{\prime}-G_{\rm
RP}^{\prime}$ and [Fe/H]′ by the loci at the same time.

Step6: Simulated $G_{\rm BP}-G$ residual defined as the difference between
$(G_{\rm BP}-G)_{\rm loci}^{\prime}$ and $(G_{\rm BP}-G)^{\prime}$ is the
output.

When simulating the binary set, we assume that the observed $G_{\rm BP}-G_{\rm
RP}$ is belong to the primary as it is the brighter component. Then the
primary mass is determined based on the PARSEC model. Mass ratio distribution
is assumed to follow a power law of index $\gamma$= 0.3 for mass ratios
between 0.05 and 1 (2013araa). Mass of the secondary is obtained after
generating a mass ratio $q$. When $M_{2}<0.08M_{\odot}$, the binary system is
supposed to be identical to the single star, where steps 1 $\sim$ 6 is
performed. When $M_{2}\geq 0.08M_{\odot}$, the intrinsic $G_{\rm BP}-G_{\rm
RP}$ of the secondary is interpolated based on the same PARSEC model assuming
the metallicity of the binary system is uniform. Next, the absolute magnitudes
of the primary and the secondary in three bands are separately calculated as
steps 1 and 2 do, except the absolute $G$ magnitudes are interpolated from the
PARSEC model (see Appendix). Then the combined absolute magnitudes of the
binary system are derived. Finally steps 3 $\sim$ 6 are applied to the
combined magnitudes and colors, $G_{\rm BP}-G$ residual are yielded.

Figure 4: Schematic illustration of the SLOT method in Gaia photometry. The
observed [Fe/H], absolute $G$ magnitude, and $G_{\rm BP}-G_{\rm RP}$ are
adopted as inputs. $G_{\rm BP}-G$ residuals of the single set and binary set
are the output. Each step is described in the main text in details.

Following the above flow chart, we could estimate the $G_{\rm BP}-G$ residual
for binary systems of different combinations of $G_{\rm BP}-G_{\rm RP}$ and
mass. In Figure 5, the residuals are estimated assuming all errors equal zero
and [Fe/H]=0. We can see that the maximum binary influence on the $G_{\rm
BP}-G$ residual is about 3 mmag, happening at the intermediate mass ratios.
That helps to set the boundaries when calculating the $\chi^{2}$.

Table 1: Fitting Coefficients of the Metallicity-dependent Stellar Color Loci

Coeff.1 |  dwarf2 |  RGB3 | Coeff. | dwarf | RGB
---|---|---|---|---|---
$a_{0}$ | 6.01E-03 |  $-$1.41E-02 | $a_{8}$ |  $-$3.08E-03 | 6.63E-04
$a_{1}$ |  $-$3.22E-03 | 1.34E-02 | $a_{9}$ | 1.62E-01 | 1.51E-02
$a_{2}$ | 5.98E-03 |  $-$4.10E-03 | $a_{10}$ |  $-$1.63E-02 | 2.45E-02
$a_{3}$ | 3.92E-03 |  $-$5.19E-04 | $a_{11}$ | 8.05E-04 |  $-$1.81E-03
$a_{4}$ | 5.70E-04 |  $-$6.29E-05 | $a_{12}$ |  $-$3.68E-02 | 5.72E-02
$a_{5}$ | 2.82E-01 | 3.75E-01 | $a_{13}$ | 3.15E-03 |  $-$8.66E-03
$a_{6}$ | 2.21E-02 |  $-$2.33E-02 | $a_{14}$ | 5.28E-03 |  $-$1.54E-02
$a_{7}$ |  $-$5.86E-03 | 7.35E-03 |  |  |

  * 1

1 $G_{\rm BP}-G$=$f(x,y)=a_{0}+a_{1}\times y+a_{2}\times y^{2}+a_{3}\times
y^{3}+a_{4}\times y^{4}+a_{5}\times x+a_{6}\times x\times y+a_{7}\times
x\times y^{2}+a_{8}\times x\times y^{3}+a_{9}\times x^{2}+a_{10}\times
x^{2}\times y+a_{11}\times x^{2}\times y^{2}+a_{12}\times x^{3}+a_{13}\times
x^{3}\times y+a_{14}\times x^{4}$, where $x$ is $G_{\rm BP}-G_{\rm RP}$ and
$y$ is [Fe/H].

  * 2

2 Coefficients are fitted under the $G_{\rm BP}-G_{\rm RP}$ range from 0.6 to
1.4 mag and [Fe/H] range from $-$1.5 to 0.5 dex.

  * 2

3 Coefficients are fitted under the $G_{\rm BP}-G_{\rm RP}$ range from 0.75 to
1.7 mag and [Fe/H] range from $-$1.5 to 0.5 dex.

Figure 5: Distributions of differences between the combined $G_{\rm BP}-G$
colors of binary systems and that predicted by the metallicity-dependent
stellar loci for the combined $G_{\rm BP}-G_{\rm RP}$ color of the system.
Binary systems are assumed to be composed of two single dwarf stars with the
solar abundance.

###  3.3 Error treatments

As described in the previous subsection, different kinds of errors are
considered in the simulation, including random and systematic errors of LAMOST
[Fe/H] and Gaia photometry.

The random errors $\sigma_{\rm ran}$([Fe/H]) are fitted using multi-observed
objects as a function of S/Ng and [Fe/H] with a third-order polynomial. Note
that we only take S/N${}_{g}<100$ in the fitting. When assigning $\sigma_{\rm
ran}$([Fe/H]) for stars with S/N${}_{g}>100$, S/N${}_{g}=100$ are used. For
S/Ng = 20, the $\sigma_{\rm ran}$([Fe/H]) are about 0.045, 0.034, and 0.027
dex for [Fe/H] = $-$1, $-$0.5, and 0 dex, respectively. For S/Ng = 50, the
$\sigma_{\rm ran}$([Fe/H]) are about 0.029, 0.022, and 0.018 dex for [Fe/H] =
$-$1, $-$0.5, and 0 dex, respectively. Systematic errors $\sigma_{\rm
sys}$([Fe/H])$=0.05-0.05\times$ [Fe/H] (2015loci2) are also considered.

Gaia EDR3 catalog provides flux errors of three Gaia bands. They are converted
into magnitude errors by the follow equations:

$\displaystyle m=-2.5\times log_{10}(F)+ZP$ $\displaystyle
m+\sigma_{m}=-2.5\times log_{10}(F-\sigma_{F})+ZP$ (3)
$\displaystyle\sigma_{m}=m+\sigma_{m}-m$

where $F$ and $\sigma_{F}$ are respectively the published flux and flux error,
$ZP$ is the photometric zero point, $\sigma_{m}$ is the random error.

The calibration errors of the photometric zero points (ZP) in the VEGAMAG,
ZPVEG = 2.8 mmag for all the three bands, are provided in phot_cont. However,
the calibration errors in the three bands are probably strongly correlated.
Therefore, the calibration errors in terms of colors should be much smaller
than involving ZPVEG of two bands independently. Besides, the calibration
errors of the ZPVEG are very likely to be magnitude-dependent, considering the
unique observation mode of Gaia and the magnitude-dependent systematic effects
discovered previously (niudr2; maw2018; casa2018).

We are mostly concerned with the effect of calibration error on the $G_{\rm
BP}-G$ residuals in this work. To get an empirical estimation of this effect,
which is represented by $\sigma_{\rm ZP,sys}$ in this work, red giant branch
(RGB) stars are used. Considering that the possibility of RGB-RGB binaries is
very small, RGB stars could be regarded as ideal single stars in terms of
photometry, since their luminosities are hundreds of times brighter than those
of dwarfs. Therefore, their $G_{\rm BP}-G$ residuals should be fully accounted
by errors of [Fe/H], random errors of the three bands, and the $\sigma_{\rm
ZP,sys}$.

Figure 6: Substitutions of the calibration error of the ZPVEG derived from the
single star model against the $G$ magnitude.

A number of 56,871 RGB stars are selected in the same way as niudr2.
Corrections of the Gaia magnitudes and colors mentioned in Section 2 are also
performed. The matallicity-dependent stellar loci of the RGB stars are listed
in Table 1 as well. In order to explore whether the $\sigma_{\rm ZP,sys}$ is
magnitude-dependent, the RGB stars are divided into bins of 0.2 mag width in
the range from 11 to 16.4 mag. For each bin, we use the single star model
shown in Figure 4 to simulate their $G_{\rm BP}-G$ residuals. The $\sigma_{\rm
ZP,sys}$ is varied from 0.0 to 2 mmag at 0.01 mmag interval and the
corresponding $\chi^{2}$ values are calculated. By choosing the global minimum
$\chi^{2}$, the best solutions of the $\sigma_{\rm ZP,sys}$ against the $G$
magnitudes are determined. Figure 6 plots the smoothed curve of the
$\sigma_{\rm ZP,sys}$. We can see that the $\sigma_{\rm ZP,sys}$ is magnitude-
dependent, showing the minimum at around $G\sim 13.5$ mag where owns the best
data quality (phot_cont). It varies around 0.75 mmag when $G<15.5$ mag, and
then increases quickly towards fainter magnitudes. Only a few RGB stars
(dwarfs as well) are of $G>16.4$ mag, therefore the value of $\sigma_{\rm
ZP,sys}$ at $G=16.4$ mag is adopted for stars of $G>16.4$ mag.

###  3.4 Correction for resolved wide binaries

In the above simulation, we have assumed that binary stars are unresolved.
However, due to the angular resolution of telescopes, 2$"$ for Gaia
(2018Arenou) and 2.5$"$ for LAMOST (liuchao2019), a binary star would be
resolved and mistakenly identified as two single stars when its two companions
have a recognizable spatial separation $\rho>2"$. The binary fractions would
be underestimated. We correct for the effect of resolved binaries in this
subsection.

According to the target selection process of LAMOST (yuan2015lamost), the
fainter secondary has a small chance of being targeted in the LAMOST survey.
Therefore, we assume that LAMOST only targets the primary stars of resolved
binaries. We use $N_{s}$, $N_{b}^{\rm unresolved}$, and $N_{b}^{\rm resolved}$
to represent the number of single stars, unresolved binary systems, and
primaries of resolved binary systems in our sample, respectively. The binary
fractions in the simulation can be described as:

$f_{b}=\frac{N_{b}^{\rm unresolved}}{N_{b}^{\rm unresolved}+N_{b}^{\rm
resolved}+N_{s}}$ (4)

Given the orbit distribution of the binary systems, the fractions of the
resolved ones could be estimated via Monte Carlo simulations. As suggested by
raghavan2010, we choose the period distribution to follow a log-normal
Gaussian profile with a mean of log$P$ = 5.03 and $\sigma_{{\rm log}P}$ =
2.28, where $P$ is in unit of day. A total mass of 1 M⊙ is also assumed. We
have verified that the result changes slightly with the assumed total mass. We
further adopt circular orbits and uniform distributions of orbital
inclinations and phases. The cut-off period against distance is firstly
calculated and plotted in black in Figure 7. Based on the above assumptions
and the angular resolution of Gaia, we calculate the fraction of the resolved
binary systems as a function of distance. It is plotted in red and can be
explained as $k=\frac{N_{b}^{\rm resolved}}{N_{b}^{\rm resolved}+N_{b}^{\rm
unresolved}}$. Most stars in our sample have distance around 0.6 kpc (see
Figure 2), so that $k$ is typically about 10$\%$. Consequently, the corrected
binary fraction for a given sample is:

$\displaystyle f_{b}^{\rm corrected}=\frac{N_{b}^{\rm resolved}+N_{b}^{\rm
unresolved}}{N_{b}^{\rm unresolved}+N_{b}^{\rm resolved}+N_{s}}$ (5)
$\displaystyle=f_{b}\times\frac{1}{1-k}$

where $k$ is the averaged value weighted by the distance distribution.

Figure 7: The black line aligned with primary $y$-axis shows the logarithmic
cut-off period as a function of distance, assuming a total mass of 1 M⊙. The
red lines aligned with the secondary $y$-axis shows the fraction of missing
resolved binaries.

##  4 Result

###  4.1 Binary fraction of the individual bin

Figure 8: Distributions of the observed $G_{\rm BP}-G$ residuals (blue
histograms) and the best-fit models (orange histograms) for the 0.8 $<m_{1}<$
1.05M⊙ sample of the individual bins of [Fe/H] and [$\alpha$/Fe] as indicated
at the top left. The derived binary fractions ${f_{b}}$, $\chi^{2}$, offsets
$\mu$, and star numbers of each bin are texted in sequence. The best-fit
single models (green solid lines) and binary models (red dashed lines) are
overplotted as well.  Figure 9: Same plots as Figure 8 but for the 0.55
$<m_{1}<$ 0.8M⊙ sample.

Figures 8 and 9 show the results of the best-simulated models for 27 bins in
the HM and 35 bins in the LM samples. In each panel, the blue histogram shows
the distribution of the observed $G_{\rm BP}-G$ residual. Green and red dashed
lines plot the modeled residual of the single and binary sets, with $f_{b}$
and $\mu$ (in unit of mmag) texted at the top left. Both sets have been
systematically shifted toward the bluer side by the offset $\mu$. Orange
histogram is the sum of the green and red dashed lines. The minimum
$\chi^{2}$, ranges of [Fe/H] and [$\alpha$/Fe], star numbers of the bin are
marked as well. We have noticed a few imperfect simulated histograms near the
end of the redder side in Figures 8 and 9. Their reasons and effects are
discussed in the next Section.

The fits are reasonably good with typical minimum $\chi^{2}$ values of 1.0
$\sim$ 5.0. The corresponding $\chi^{2}$ distributions on the $f_{b}$ and
$\mu$ panel of each the individual bin are demonstrated in Figures 10 and 11,
where red stars donate the resultant set of $f_{b}$ and $\mu$, error bars
indicate the error of $f_{b}$ from the Bootstrap Sample method. To highlight
the result, we only contour sets with $\chi^{2}$ less than 1.4 times the
minimum $\chi^{2}$. The general $\mu$ are quite small, most are less than 1
mmag. It tells that the $f_{b}$ and $\mu$ are tightly correlated in each the
individual bin, the higher the $f_{b}$, the larger the $\mu$, as we expected.
As a result, the global minimum could be simply figured out through a
ellipsoid fitting.

Figure 10: Same [Fe/H] and [$\alpha$/Fe] bins as Figure 8 for the HM sample,
distributions of $\chi^{2}$ between the observed $G_{\rm BP}-G$ residual and
model predicted $G_{\rm BP}-G$ residual for combinations of different
unresolved binary fractions ${f_{b}}$ and offsets $\mu$. Red stars represent
the global $\chi^{2}$ minimum. Contour labels indicate the $\chi^{2}$. For
good looking, only combinations with $\chi^{2}$ less than 1.4 times the
minimum $\chi^{2}$ are plotted in the individual bins. Figure 11: Same plots
as Figure 10 but for the 0.55 $<m_{1}<$ 0.8M⊙ sample. Figure 12: Summary of
results for all bins of the 0.8 $<m_{1}<$ 1.05M⊙ sample. For each box,
simulated unresolved binary fraction $f_{b}$ and bias-corrected binary
fraction $f_{b}^{\rm corrected}$ are texted in black at the left. Offset $\mu$
in the unit of mmag is in red at the right. Orange and turquoise frames
separately donate the thin and thick disks. The $f_{b}^{\rm corrected}$ as
well as formal error of every [Fe/H] interval is also texted in the same
color.  Figure 13: Same plots but for the 0.55 $<m_{1}<$ 0.8M⊙ sample.

Given $f_{b}$, $f_{b}^{\rm corrected}$ can be subsequently obtained using
Equation 5 and the curve shown in Figure 7. We summarize $f_{b}$ and the
corresponding formal error, $f_{b}^{\rm corrected}$, and $\mu$ in Figures 12,
13 and Table 2. For each the individual bin, its ranges of [Fe/H] and
[$\alpha$/Fe] are illustrated by a box in Figures 12, 13 and listed in Table
2. The $f_{b}$, $f_{b}^{\rm corrected}$, and $\mu$ of each bin are texted
within the box as exampled by the inset plot. The typical complement of the
resolved binaries are 2 $\sim$ 6$\%$ for the HM sample and 3 $\sim$ 7$\%$ for
the LM sample depending on their distance distributions. The 2D histograms
imply the number distribution of [Fe/H] and [$\alpha$/Fe], where the thin and
thick disks could be easily identified. We use colors of orange and teal to
mark the thin and thick disk, respectively. The [Fe/H]$<-1$ bin in the LM
sample is counted as the inner halo. Notice that the LM and HM sample cover
different ranges of [Fe/H] and [$\alpha$/Fe]. So do the scales of Figures 12
and 13.

Table 2: Result of the individual bin

[Fe/H]l |  [Fe/H]u |  $[\alpha/Fe]$u |  $[\alpha/Fe]$l | $\mu^{1}$ | $f_{b}$ |  error2 | $f_{b}^{\rm corrected}$ | Mass
---|---|---|---|---|---|---|---|---
$-$0.4 |  $-$0.3 | 0.125 | 0.033 | 0.520 | 0.37 | 0.07 | 0.41 | HM
$-$0.4 |  $-$0.3 | 0.200 | 0.125 | 0.530 | 0.44 | 0.07 | 0.49 | HM
$-$0.4 |  $-$0.3 | 0.300 | 0.200 | 0.640 | 0.23 | 0.11 | 0.25 | HM
$-$0.3 |  $-$0.2 | 0.100 | 0.025 | 0.550 | 0.44 | 0.05 | 0.49 | HM
$-$0.3 |  $-$0.2 | 0.150 | 0.100 | 0.510 | 0.46 | 0.05 | 0.51 | HM
$-$0.3 |  $-$0.2 | 0.200 | 0.150 | 0.550 | 0.42 | 0.08 | 0.47 | HM
$-$0.3 |  $-$0.2 | 0.275 | 0.200 | 0.590 | 0.43 | 0.09 | 0.48 | HM
$-$0.2 |  $-$0.1 | 0.075 | 0.000 | 0.510 | 0.38 | 0.04 | 0.43 | HM
$-$0.2 |  $-$0.1 | 0.125 | 0.075 | 0.500 | 0.45 | 0.04 | 0.51 | HM
$-$0.2 |  $-$0.1 | 0.175 | 0.125 | 0.530 | 0.45 | 0.06 | 0.50 | HM
$-$0.2 |  $-$0.1 | 0.250 | 0.175 | 0.360 | 0.36 | 0.10 | 0.40 | HM
$-$0.1 | 0.0 | 0.050 |  $-$0.025 | 0.460 | 0.33 | 0.04 | 0.37 | HM
$-$0.1 | 0.0 | 0.100 | 0.050 | 0.450 | 0.40 | 0.03 | 0.45 | HM
$-$0.1 | 0.0 | 0.150 | 0.100 | 0.480 | 0.41 | 0.05 | 0.46 | HM
$-$0.1 | 0.0 | 0.200 | 0.150 | 0.480 | 0.22 | 0.10 | 0.25 | HM
0.0 | 0.1 | 0.025 |  $-$0.040 | 0.400 | 0.37 | 0.05 | 0.42 | HM
0.0 | 0.1 | 0.075 | 0.025 | 0.450 | 0.40 | 0.03 | 0.45 | HM
0.0 | 0.1 | 0.125 | 0.075 | 0.490 | 0.40 | 0.05 | 0.45 | HM
0.0 | 0.1 | 0.200 | 0.125 | 0.450 | 0.46 | 0.11 | 0.51 | HM
0.1 | 0.2 | 0.025 |  $-$0.040 | 0.330 | 0.24 | 0.04 | 0.27 | HM
0.1 | 0.2 | 0.075 | 0.025 | 0.450 | 0.39 | 0.03 | 0.44 | HM
0.1 | 0.2 | 0.150 | 0.075 | 0.440 | 0.51 | 0.07 | 0.57 | HM
0.2 | 0.3 | 0.025 |  $-$0.040 | 0.420 | 0.29 | 0.05 | 0.33 | HM
0.2 | 0.3 | 0.100 | 0.025 | 0.530 | 0.42 | 0.04 | 0.47 | HM
0.3 | 0.4 | 0.025 |  $-$0.040 | 0.520 | 0.33 | 0.06 | 0.37 | HM
0.3 | 0.4 | 0.100 | 0.025 | 0.590 | 0.50 | 0.05 | 0.56 | HM
0.4 | 0.5 | 0.100 |  $-$0.025 | 0.550 | 0.45 | 0.07 | 0.51 | HM
$-$1.5 |  $-$1.0 | 0.375 | 0.200 | 0.60 | 0.81 | 0.09 | 0.91 | LM
$-$1.0 |  $-$0.9 | 0.375 | 0.200 | 0.57 | 0.48 | 0.13 | 0.54 | LM
$-$0.9 |  $-$0.8 | 0.375 | 0.200 | 0.58 | 0.55 | 0.09 | 0.62 | LM
$-$0.8 |  $-$0.7 | 0.275 | 0.200 | 0.58 | 0.43 | 0.09 | 0.48 | LM
$-$0.8 |  $-$0.7 | 0.375 | 0.275 | 0.59 | 0.40 | 0.08 | 0.45 | LM
$-$0.7 |  $-$0.6 | 0.200 | 0.100 | 0.71 | 0.46 | 0.13 | 0.52 | LM
$-$0.7 |  $-$0.6 | 0.275 | 0.200 | 0.60 | 0.49 | 0.07 | 0.55 | LM
$-$0.7 |  $-$0.6 | 0.375 | 0.275 | 0.65 | 0.51 | 0.07 | 0.58 | LM
$-$0.6 |  $-$0.5 | 0.200 | 0.050 | 0.61 | 0.52 | 0.06 | 0.59 | LM
$-$0.6 |  $-$0.5 | 0.275 | 0.200 | 0.59 | 0.43 | 0.05 | 0.49 | LM
$-$0.6 |  $-$0.5 | 0.350 | 0.275 | 0.57 | 0.45 | 0.08 | 0.51 | LM
$-$0.5 |  $-$0.4 | 0.150 | 0.050 | 0.68 | 0.39 | 0.05 | 0.44 | LM
$-$0.5 |  $-$0.4 | 0.200 | 0.150 | 0.62 | 0.46 | 0.06 | 0.52 | LM
$-$0.5 |  $-$0.4 | 0.275 | 0.200 | 0.58 | 0.42 | 0.05 | 0.47 | LM
$-$0.5 |  $-$0.4 | 0.350 | 0.275 | 0.45 | 0.33 | 0.10 | 0.37 | LM
$-$0.4 |  $-$0.3 | 0.125 | 0.050 | 0.72 | 0.35 | 0.05 | 0.40 | LM
$-$0.4 |  $-$0.3 | 0.175 | 0.125 | 0.65 | 0.40 | 0.05 | 0.45 | LM
$-$0.4 |  $-$0.3 | 0.225 | 0.175 | 0.58 | 0.39 | 0.06 | 0.44 | LM
$-$0.4 |  $-$0.3 | 0.325 | 0.225 | 0.53 | 0.46 | 0.06 | 0.52 | LM
$-$0.3 |  $-$0.2 | 0.100 | 0.025 | 0.73 | 0.30 | 0.05 | 0.34 | LM
$-$0.3 |  $-$0.2 | 0.150 | 0.100 | 0.67 | 0.37 | 0.05 | 0.42 | LM
$-$0.3 |  $-$0.2 | 0.200 | 0.150 | 0.59 | 0.34 | 0.07 | 0.38 | LM
$-$0.3 |  $-$0.2 | 0.300 | 0.200 | 0.38 | 0.31 | 0.07 | 0.35 | LM
$-$0.2 |  $-$0.1 | 0.075 | 0.000 | 0.71 | 0.27 | 0.05 | 0.31 | LM
$-$0.2 |  $-$0.1 | 0.125 | 0.075 | 0.69 | 0.30 | 0.05 | 0.34 | LM
$-$0.2 |  $-$0.1 | 0.175 | 0.125 | 0.74 | 0.44 | 0.08 | 0.50 | LM
$-$0.2 |  $-$0.1 | 0.250 | 0.175 | 0.51 | 0.40 | 0.12 | 0.45 | LM
$-$0.1 | 0.0 | 0.050 |  $-$0.025 | 0.65 | 0.25 | 0.06 | 0.29 | LM
$-$0.1 | 0.0 | 0.100 | 0.050 | 0.81 | 0.35 | 0.05 | 0.40 | LM
$-$0.1 | 0.0 | 0.175 | 0.100 | 0.61 | 0.32 | 0.10 | 0.36 | LM
0.0 | 0.1 | 0.050 |  $-$0.025 | 0.73 | 0.20 | 0.07 | 0.23 | LM
0.0 | 0.1 | 0.125 | 0.050 | 0.84 | 0.25 | 0.08 | 0.28 | LM
0.1 | 0.2 | 0.050 |  $-$0.025 | 0.93 | 0.26 | 0.09 | 0.30 | LM
0.1 | 0.2 | 0.125 | 0.050 | 1.10 | 0.32 | 0.13 | 0.36 | LM
0.2 | 0.3 | 0.125 |  $-$0.025 | 1.19 | 0.19 | 0.16 | 0.22 | LM

  * 1

1 $\mu$ is in the unit of mmag.

  * 2

2 Formal error deduced from the Bootstrap Sample method.

###  4.2 Explore the dependency of the stellar parameters

We are now ready to explore the impacts of stellar parameters on the binary
fractions. Current studies based on both observations (e.g. 2013araa; gao2014;
2015loci2; moe2019anti; aopgee2020alpha) and theoretical simulations (e.g.
bate2019; hurley2002) have shown that the mass of the primary, metallicity
[Fe/H], and [$\alpha$/Fe] contribute to the binary fractions. Examining the
impact from any of them on the binary fractions needs decomposing the
degenerate parameters. For instance, investigating [$\alpha$/Fe] within a
mono-mass and mono-metallicity sample. Figure 14 plots the bias-corrected
binary fraction $f_{b}^{\rm corrected}$ against [$\alpha$/Fe] of both LM (red
solid lines) and HM (black dashed lines) samples that are equally binned into
12 [Fe/H] intervals from $-$0.8 to 0.4 dex. When studying one parameter, other
parameters are within the tolerance that their contributions on the binary
fraction are tiny. The $f_{b}^{\rm corrected}$ values of points in Figure 14
are identical to those in Figures 12 and 13. The vertical error bars are same
as the horizontal error bars in Figures 8 and 9. The [$\alpha$/Fe] values of
points are the median values of all stars contained in the individual bins.
Filled dots are for the thin disk and open circles for the thick disk.

Among the 12 panels, 6 panels contain points from both LM and HM samples with
comparable [$\alpha$/Fe] and [Fe/H], showing that the binary fraction
increases toward higher mass, which has been widely confirmed by many
investigations (e.g. 2015loci2; 2013araa; moe2017).

We tend to explore the relationship between [$\alpha$/Fe] and binary fractions
under situations of the thin and thick disks. Nine panels ($-0.5<$ [Fe/H]
$<0.4$ dex) contain two thin disk points (filled dots), 8 pairs from the HM
sample and 7 pairs from the LM sample. We generate a Gaussian likelihood
distribution of the binary fraction for each point based on its error. For
each pair, we estimate the possibility that the point with a larger
[$\alpha$/Fe] owns a larger binary fraction than that owned by the point with
a smaller [$\alpha$/Fe], by comparing their likelihood distributions. From
metal poor to metal rich, the possibilities are 0.77, 0.60, 0.89, 0.92, 0.69,
1.0, 0.98, and 0.98 for the HM sample, 0.81, 0.77, 0.85, 0.68, 0.91, 0.70, and
0.66 for the LM sample. Therefore, we think there is a strong positive
correlation between [$\alpha$/Fe] and the binary fraction of the thin disk
stars within within a mono-mass and mono-metallicity sample.

While in the thick disk, this tendency is not clear, some show positive
correlation while others show negative. Also, the relationship of the binary
fractions between the thin and thick disk is hard to simply recognize.

Tendency on the $\alpha$ elements is consistent with the result from tian2018,
in which LAMOST DR4 are used to select binaries with RV variations. They find
that stars with high [$\alpha$/Fe] have higher binary fraction than that of
stars with low [$\alpha$/Fe] within a mono-temperature (by proxy, mass) and
mono-metallicity sample. They speculate that it is probably related to the age
of stars, suggesting that earlier forming stars with high [$\alpha$/Fe]
possess a higher binary fraction than newly forming ones with low
[$\alpha$/Fe].

Figure 14: Binary fractions of mass range 0.55 $<m_{1}<$ 0.8M⊙ (red solid
lines) and 0.8 $<m_{1}<$ 1.05M⊙ (black dashed lines) against [$\alpha$/Fe] at
each [Fe/H] range. Grids from the thin and thick disks are respectively
scattered in filled and hollow dots.

However, estimating precise ages of the dwarf stars is kind of a world-class
problem. To double check whether the relation between $\alpha$ elements and
binary fractions is due to age, evidences of relationships between age,
metallicity, and $\alpha$ elements from other populations of stars are
demanded. It has been well studied with data from different surveys.
2019Feuillet use RGB stars from APOGEE data and explore the [$\alpha$/M]-age
relations at different spatial locations in the Milky Way disk. From their
Figure 6, we can see that [$\alpha$/M]-age relation is nearly flat in the
range of [$\alpha$/M] $>$ 0.15, which could be counted as the older thick disk
stars. But for the younger thin disk, [$\alpha$/M] systematically increase
with increasing age. Meanwhile, xiang2017; wu2019; 2020ApJS..249...29H present
the same result using MSTO, RGB stars, and Red clump stars from the LAMOST
data, respectively.

Likely, our sample shows a positive correlation between [$\alpha$/Fe] and
binary fraction for the thin disk stars but not for the thick disk stars.
Combing these two phenomena, it is reasonable to suggest that the binary
fractions has a positive correlation with the stellar age with hints of the
$\alpha$ abundances, at least for the field F-K dwarf stars in the solar
vicinity. The different dependence of binary fractions in the thin and disks
also imply that they have very different chemical enrichment histories, as
validated by many published works (e.g. 2001Matteucci; agemetal; xiang2015).

Our sample also covers wide [Fe/H] ranges in the thin and thick disks. Binary
fractions $f_{b}^{\rm corrected}$([Fe/H]) at 1 dex intervals of [Fe/H] in both
thin and thick disks are then calculated as well as errors. They are texted at
the top and bottom of Figure 12 and 13 in orange and teal. Figure 15 scatters
$f_{b}^{\rm corrected}$([Fe/H]). Grey and brown colors donate different mass
ranges. $f_{b}^{\rm corrected}$([Fe/H]) roughly decreases with increasing
[Fe/H] as we have known before (e.g. 2015loci2; moe2019anti).

Figure 15: Binary fractions $f_{b}^{\rm corrected}$([Fe/H]) for the thin disk
(top panel) and thick disks (bottom panel). Pink and gray dashed lines
separately represent their tendencies of the LM and HM samples in each panel.

We apply linear fittings to points having [$\alpha$/Fe] decrease with
increasing [Fe/H] and performing 2-sigma clipping. Results are drawn in red
and black dashed lines in Figure 15. Slopes of lines are marked beside.
Obviously, $f_{b}^{\rm corrected}$([Fe/H]) in the thin disk suffers a steeper
slope than that in the thick disk, especially for the LM sample. The steeper
slope of the thin disk is contributed by both increasing [Fe/H] and decreasing
[$\alpha$/Fe] since we have suggested a positive correlation between
[$\alpha$/Fe] and the binary fraction. While in the thick disk, where the
decreasing binary fraction only comes from the increasing [Fe/H] given that
[$\alpha$/Fe] do not affect the binary fraction. For the HM sample, the trend
is not outstanding because of the fewer points in the thick disk.

For [Fe/H] $\geq$ 0.2, as [Fe/H] increasing, the binary fractions are constant
instead of decreasing. We suppose this might be due to three possible reasons:
(1) [$\alpha$/Fe] has no contribution for a roughly constant distribution; (2)
only HM sample left where mass may plays an role; (3) the super-metal-rich
stars ([Fe/H] $>0.2$ by convention smr2001) are more closer to the galactic
center (Andrievsky2002; Anders2014) and have significantly larger probability
of radial migration (2018Minchev; 2009Chiappini; 2009Asplund) than others.
Therefore the formation and evolution process could be different consequently
the binary fraction.

Table 3: Binary Fractions of the Field in Galactic Disks and Halos  | 0.8 $<m_{1}<$ 1.05M⊙ | 0.55 $<m_{1}<$ 0.8M⊙
---|---|---
Thin Disk | $0.43\pm 0.01$ | $0.39\pm 0.02$
Thick Disk | $0.46\pm 0.02$ | $0.49\pm 0.02$
Halo | - | $0.91\pm 0.09$
Total | $0.44\pm 0.01$ | $0.42\pm 0.01$

To sum up, the averaged binary fractions $f_{b}^{\rm corrected}$ and their
errors of the HM samples are $0.43\pm 0.01$ for the thin disk and $0.46\pm
0.02$ for the thick disk. As for the LM sample, $f_{b}^{\rm corrected}$ is
$0.39\pm 0.02$ for the thin disk and $0.49\pm 0.02$ for the thick disk. For
0.8 $<m_{1}<$ 1.05M⊙, thick disk owns a little bit higher binary fraction than
that of thin disk in total. For the lower mass range, $f_{b}^{\rm corrected}$
in the thick disk is clearly higher, which could be understood by the larger
difference of distributions of abundances in the thin and thick disks. Binary
fraction of halo has great probability to be larger than others. Notice that
there are less objects in the halo in our sample. The averaged binary
fractions $f_{b}^{\rm corrected}$ of the disk stars are $0.44\pm 0.01$ for the
HM sample and $0.42\pm 0.01$ for the LM sample.

##  5 Discussion

###  5.1 Redder excess of the color residual

In this paper, based on the estimations of the binary effect on the $G_{\rm
BP}-G$ residuals, the lower and upper limits of the $G_{\rm BP}-G$ residuals
in the fitting produce are respectively $-7$ and 4 mmag. Objects beyond bounds
are ignored in our method yet worthy to be discussed. The median ignoring
fraction of objects at $G_{\rm BP}-G<-7$ mmag is 0.3 per cent for both HM and
LM samples. While at $G_{\rm BP}-G>4$ mmag, the median ignoring fraction is 5
per cent for both HM and LM samples. That leads several imperfect fitting
histograms near the end of the redder side in Figures 8 and 9. Reasons can be
attributed to two: (1) contamination from the background and (2) flux excess
due to the unresolved multiple systems. The first effect should be symmetric
and barely changes the results since only 0.3 per cent at the bluer side. The
second influence is expected to be same as effects from the binary and even
greater.

However, it is hard to give a convinced and specific modification on the
binary fractions we obtained above when considering effects from the multiple
systems. Even though the fractions of them have been studied (e.g.
raghavan2010; dm91; 2012bate; 2014Tokovinin1) and conducted to be $0.11\sim
0.17$ for the solar mass range in general, it changes with the period
(2012allen; 2008Makarov) in a complex way. Works also reveal that the
structures of multiplicity are variegated and require further research. The
mass ratio distribution of the short and long period subsystems are
distinguishing (raghavan2010; 2008Tokovinin; 2012bate). The period
distribution of subsystems varies with different kinds of hierarchies
(2014Tokovinin).

Meanwhile, a substantial portion of the multiple systems contain the
third/fourth components of little mass, whose performances on color residuals
are close to the binary systems. These kind of systems have been accounted in
our results. We suppose that the cut-off of the fitting range barely affect
our results.

###  5.2 Comparison with other works

Comparing with works based on other approaches, the biggest advantage of our
method is the less limitations on the binary period. Summarized illustrations
of the ranges of binary periods and mass ratios that can be detected with
different methods are shown in moe2017 and 2018specbi. For instance, as one of
the most widely used techniques, the spectroscopic binaries and eclipsing
binaries could just be capable of close binary. The former is suitable for
Log10($P$) $\leq 3$ day. The later gathers at even a shorter period of
Log10($P$) $\leq 2$ day. The common proper motion, X-ray, and visual binaries
are only sensitive to wide binaries. The spectral fitting method could achieve
Log10($P$) $\leq 8$ day but only works for the intermediate mass ratio. In
addition to periods, other methods also suffer from other specific
limitations. Taking several examples, spectroscopic technique requires close
binaries with having time-domain data. Photometric technique is confined to
high inclination systems. Astrometric techniques like from Hipparcos and Gaia
are biased to the bright binaries having intermediate separations based on the
angle resolution and time baseline of the serveries.

As for the SLOT method, its realization is solely taking advantage of the
difference between the color deduced by the metallicity-dependence stellar
loci (of the single star) and combined color (of binary system). It covers
both close and wide binaries, supporting Log10($P$) less than $6\sim 8$ day
depending on the distance distribution and instrument spatial resolution.
Besides, it is independent of the mass ratio distribution and efficient for
single epoch data.

We stress that most samples of other binary fraction studies are restricted to
the thin disk stars due to the selection bias. It would be reasonable to
compare their results with ours of the thin disk. The thin disk stars of the
HM sample have the similar spectra type range with the G2–K3 sample in
raghavan2010. Their result $0.41\pm 0.03$ agrees well with ours $0.43\pm
0.01$. Our work basically uses the the same method as 2015loci2 except for
taking the resolved binary into consideration. The HM sample roughly
correspond to the $0.6<g-i<0.9$ mag and $-0.5<$ [Fe/H] $<0.5$ sample in
2015loci2, where the binary fractions is 39$\%$. The LM sample roughly
correspond to the $0.9<g-i<1.2$ mag and $-1<$ [Fe/H] $<0$ sample, where the
binary fractions is 35$\%$. The average fractions of the resolved binary in
our work is 39$\%$ for the HM sample and 37$\%$ for the LM sample. Two results
agree well supposing the spatial resolutions of Gaia and SDSS are comparable.

The spectroscopic APOGEE survey has been widely used in detecting binary
system and exploring the close binary fraction ($a<$ 10 AU) (2018badnes;
2018specbi; aopgee2020alpha; moe2019anti). They have suggested an
observational result that the negative correlation between the stellar
multiplicity and [Fe/H] weakens and even disappears with the increasing
period. Our result also supports this conclusion. In aopgee2020alpha;
moe2019anti, the close binary fraction is about 0.22 at [Fe/H]=$-$0.25 and
halve at [Fe/H]=0.25. While for the binary fraction in this work without
limitations on the period, it is about 0.4 at [Fe/H]=$-$0.25 and 0.3 at
[Fe/H]=0.25. That supports the weaker negative correlation for the wide binary
compared to the close binary.

We also explore the impact of [$\alpha$/Fe] on the binary fraction, which is
insufficient in the previous researches and worthy to be explored. In
aopgee2020alpha, an anti-correlation between $\alpha$ abundances and
multiplicity is found and substantially steeper than that observed for Fe,
which is contrary to our results. The opposite result might be due to the
differences of the period distributions between samples in this work and
aopgee2020alpha. As mentioned, they focus on the close binary exclusively,
roughly Log10($P$) $\leq$ 4 day. While our work covers a larger volume of both
wide and close binary. Close binaries are believed to form via disk
fragmentation, while wide binaries are via fragmentation of molecular
cores/filaments (2001White; 2007Goodwin; moe2017). As we have known that the
relation between the binary fraction and [Fe/H] could change with the binary
period, we conjecture that it happens to [$\alpha$/Fe] as well. We think it is
tightly related to the binary formation and evolution and should be carefully
investigated in the future. Up to now, the relationship between the binary
fraction and [$\alpha$/Fe] varies from a negative correlation to a positive
correlation with the increasing binary period.

##  6 Summary

In this work, the binary fraction of individual stellar population is derived
using the SLOT method and corrected for the effect of resolved binaries. We
use a 0.2 million volume-limited sample of the filed late F, G, and early K
dwarf stars selected from Gaia and LAMOST data, in which Gaia photometery has
been delicately re-calibrated by niudr3 and yl. We divide the sample into two
mass ranges and further divide each subsample into 62 bins based on their
[Fe/H] and [$\alpha$/Fe]. All bins are classified into the thin disk, thick
disk, and halo when exploring tendencies of the binary fractions. The bias-
corrected binary fraction is $0.44\pm 0.01$ for the 0.8 $<m_{1}<$ 1.05M⊙
sample and $0.42\pm 0.01$ for the 0.55 $<m_{1}<$ 0.8M⊙ sample. The results
suggest that the Galactic thin and thick disks have comparable binary
fractions. Halo stars have a great probability of the highest binary fraction.

By investigating the mono-metallicity and mono-mass sample, [$\alpha$/Fe] is
found to have a positive correlation with the binary fraction in the thin disk
and become independent in the thick disk. We suppose it is caused by the
distributions of the stellar age. Plenty of works (2014Bergemann; xiang2017;
Mackereth2018; wu2021) that conducted from diverse data sources have revealed
that stars of intermediate-to-old ages from the thick disk show a contiguous
distribution of age within a mono-metallicity sample. While for the young
stars in the thin disk, stellar ages clearly exhibit a gradient with
[$\alpha$/Fe]. Given the binary fractions as a function of [$\alpha$/Fe], we
further suggest that investigating the changing of the binary fractions with
[Fe/H] should subtracting the contributions from [$\alpha$/Fe]. We find the
steeper observed negative correlation between the binary fraction and [Fe/H]
in the thin disk compared with that in the thick disk. By comparing with the
slope of the close binary (moe2019anti; aopgee2020alpha), our result agrees
with the previously announced phenomenon that the effect from the [Fe/H] on
the binary fractions changes with the increasing binary period.

The SLOT method is described in detail in the Section 3. We compare the
approximate parameter space of detection abilities of the SLOT method with
that of other methods, such as spectroscopic binaries and eclipsing binaries.
We stress that the biggest advantage of our method is less limitations on the
period. It suffers smaller restrictions on epoch data, mass ratio, and
inclination as well.

The mass ratio distribution is assumed to follow a single-power-law of index
$\gamma$= 0.3, although some recent works suggest that the index varies with
stellar mass and metallicity (liuchao2019) and a more complex broken-power-law
is preferred (moe2017; 2019gaiatwin). We are going to investigate mass ratio
distribution in the future by combining both the color offsets and magnitude
offsets together.

We acknowledge helpful discussions with Prof. Liu Chao, Zhang Haopeng, and
Wang Yilun. This work is supported by National Science Foundation of China
(NSFC) under grant numbers 11603002, 11988101, and 113300034, National Key
Research and Development Program of China (NKRDPC) under grant numbers
2016YFA0400804, 2019YFA0405503, and 2019YFA0405504. This work has made use of
data products from the Guoshoujing Telescope (the Large Sky Area Multi-Object
Fiber Spectroscopic Telescope, LAMOST). 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 the National Astronomical Observatories, Chinese
Academy of Sciences. 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.

##  Appendix A PARSEC model

The PARSEC model is applied in this work when calculating the stellar mass and
absolute $G$ magnitude. We extract isochrones of different sets of [Fe/H] with
step of 0.1 dex from $-$1.5 to 0.5 dex and age with step of 1 Gyr from 1 to 10
Gyr. In Figure A16, grey points are the sample we used, colored tracks are the
color-magnitude tracks from the PARSEC model. Red lines that laying on the
PARSEC model are empirically fitted using the grey points, in which 2-sigma
clipping is performed. Some grey points scatter above the red line as well as
the tracks are binaries. Our sample is in good agreement with the PARSEC
model. We can see that the width of tracks caused by ages among the color
ranges of our sample is narrow. Hence we just interpolate the tracks to have
the $M_{\rm G}$.

Figure A17 is the color-mass model of the PARSEC model. Stellar mass is
derived in the same way as the $M_{\rm G}$.

Figure A16: Grey: dwarf samples used in this work. Red lines: fitted
$M_{G}=f(G_{\rm BP}-G_{\rm RP},{\rm[Fe/H]})$ relation as texted in the figure.
Colored tracks: Color-magnitude model of the PARSEC isochrones used in this
work. Ages are indicated in the color bar. Metallicities are labelled in the
top center of each panel. Variations of the relation of different
metallicities are displayed in the last panel.  Figure A17: Color-mass model
of the PARSEC isochrones used in this work. Ages are indicated in the color
bar. Metallicities are labelled in the top center of each panel.