11institutetext:  Department of Physics, Hebei Normal University, Shijiazhuang
050024, China 22institutetext: Guo Shoujing Institute for Astronomy, Hebei
Normal University, Shijiazhuang 050024, China 33institutetext: CAS Key
Laboratory of Optical Astronomy, National Astronomical Observatories, Beijing
100101, China 44institutetext: University of Chinese Academy of Sciences,
Beijing 100049, China 55institutetext: Nanjing Institute of Astronomy and
Optical Technology, University of Chinese Academy of Sciences, Nanjing 211135,
China 66institutetext: School of Physics, Astronomy and Mathematics,
University of Hertfordshire, College Lane, Hatfield AL10 9AB, UK

# Deep learning interpretability analysis for carbon star identification in
$Gaia$ DR3

Shuo Ye 112233 Wen-Yuan Cui , Corresponding author,
<EMAIL_ADDRESS>Yin-Bi Li , Corresponding author,
<EMAIL_ADDRESS>A-Li Luo , Corresponding author,
<EMAIL_ADDRESS>Hugh R. A. Jones 66

###### Abstract

Context. A large fraction of Asymptotic Giant Branch (AGB) stars develop
carbon-rich atmospheres during their evolution. Based on their color and
luminosity, these carbon stars can be easily distinguished from many other
kinds of stars. However, a large number of G, K, and M giants are also
distributed in the same region as carbon stars on the HR diagram. Their
spectra have differences, especially in the prominent CN molecular bands.

Aims. We aim to distinguish carbon stars from other kinds of stars using
$Gaia$’s XP spectra, while providing attribution explanations of key features
necessary for identification, and even discovering additional new spectral key
features.

Methods. We proposed a classification model named ‘GaiaNet’, an improved one-
dimensional convolutional neural network specifically designed for handling
$Gaia$’s XP spectra. We utilized the SHAP interpretability model to calculate
the SHAP value for each feature point in a spectrum, enabling us to explain
the output of the ‘GaiaNet’ model and provide further meaningful analysis.

Results. Compared to four traditional machine-learning methods, the ‘GaiaNet’
model exhibits an average classification accuracy improvement of approximately
0.3% on the validation set, with the highest accuracy even reaching 100%.
Utilizing the SHAP algorithm, we present a clear spectroscopic heatmap
highlighting molecular band absorption features primarily distributed mainly
around CN${}_{\text{773.3}}$ and CN${}_{\text{895.0}}$, and summarize five
crucial feature regions for carbon star identification. Upon applying the
trained classification model to the CSTAR sample with $Gaia$ ‘xp_sampled_mean’
spectra, we obtained 451 new candidate carbon stars as a by-product.

Conclusions. Our algorithm is capable of discerning subtle feature differences
from low-resolution spectra of $Gaia$, thereby assisting us in effectively
identifying carbon stars with typically higher temperatures and weaker CN
features, while providing compelling attributive explanations. The
interpretability analysis of deep learning holds significant potential in
spectral identification.

###### Key Words.:

Method: deep learning; Method: interpretability analysis; stars: Carbon; Data:
$Gaia$ DR3 xp_sampled_mean_spectrum.

##  1 Introduction

Carbon stars, first recognized by Secchi (1869), exhibit an inversion of the
C/O ratio (C/O ¿ 1). Compared to other stars, carbon stars possess carbon-
enriched atmospheres . The enrichment can originate from mass transfer in
binary systems or the activation of the third dredge-up process during their
AGB evolution phase (Gaia Collaboration et al. 2023). The majority of carbon
stars in the galaxy are binary systems wherein the host star inherited
atmospheric carbon through mass transfer from its AGB companion, which is now
a white dwarf (Ji et al. 2016). The carbon enrichment of the cool and luminous
carbon stars (C-N type) is due to the pollution by nuclear He fusion products
from the inner to the outer layers in their AGB phase, which are called
intrinsic carbon stars (Ji et al. 2016). This process transports carbon from
the interior region to the stellar surface, leaving the outer layers rich in
carbon-containing molecules and dust particles. Thus, their spectra show
C${}_{\text{2}}$ and CN molecular bands stronger than usual in stars cooler
than 3,800 K (i.e. $G_{\mathrm{BP}}-G_{\mathrm{RP}}\geq 2$) (Gaia
Collaboration et al. 2023).

Carbon stars exhibit strong absorption characteristics due to molecular bands
of carbon-containing compounds such as C${}_{\text{2}}$, CN, CH, and various
hydrocarbons. The presence and intensity of these spectral features provide
valuable insights into the atmospheric conditions and chemical processes
within these stars. Carbon stars are also important contributors to the
interstellar medium and serve as good reference cases for studying the
physical processes affecting the end of the life of low-mass stars (Gaia
Collaboration et al. 2023).

Traditionally, identifying and classifying carbon stars relied on parametric
measurements and manual checking of their spectra. Ji et al. (2016) identified
894 carbon stars from LAMOST DR2 by measuring multiple line indices from the
stellar spectra. Abia et al. (2022) reported the identification of 2,660 new
carbon star candidates through 2MASS photometry, $Gaia$ astrometry, and their
location in the $Gaia$–2MASS diagram. Li et al. (2023) distinguished carbon
stars from M-type giants by selecting the criterion of spectral indices
satisfying [$\rm{CaH3}-0.8\times\rm{CaH2}-0.1$]$<0$. Gaia Collaboration et al.
(2023) screened carbon stars by measuring the carbon molecular band head
strength. Lebzelter et al. (2023) classified 546,468 carbon star candidates
using a method guided by narrow-band photometry (Palmer & Wing 1982). Their
SOS module examines this feature in an automated way by computing the pseudo-
wavelength difference between the two highest peaks in each RP spectrum of a
given source, taking the median value of the results, and storing it in the
parameter of median_delta_wl_rp. They assumed all stars with
median_delta_wl_rp$>7$ to be C rich (i.e., isCstar = TRUE). Li et al. (2024)
identified 3546 carbon stars through line index and near-infrared color–color
diagrams. Through visual inspection of these carbon star spectra, they further
subclassify them into C–H, C–R, C–N, and Ba stars.

To date, many machine-learning methods have been used to identify carbon stars
based on their spectra. Si et al. (2014) applied the label propagation
algorithm to search for 260 new carbon stars from SDSS (Sloan Digital Sky
Survey) DR8. Si et al. (2015) applied the efficient manifold ranking algorithm
to search for 183 carbon stars from the LAMOST pilot survey. Li et al. (2018)
identified 2,651 carbon stars in the spectra of more than seven million stars
using an efficient machine-learning algorithm for LAMOST (Large Sky Area
Multi-Object Fiber Spectroscopic Telescope) DR4. Sanders & Matsunaga (2023)
investigated the use of unsupervised learning algorithms to classify the
chemistry of long-period variables (Lebzelter et al. 2023) from $Gaia$ DR3’s
BP/RP spectra (also called XP spectra; Carrasco et al. 2021) into O-rich and
C-rich groups. They also employed a supervised approach to separate O-rich and
C-rich sources using broadband optical and infrared photometry. In all, they
tagged a total of 23,737 C-rich classifications based on the BP/RP spectra and
identified a small population of C-rich stars in the Galactic bar-bulge
region.

In recent years, the application of interpretable analysis techniques in
astronomy has shown great potential. Qin et al. (2019) used a random forest
(RF) algorithm deriving rank features, then picked out 15,269 Am candidates
from the early-type stars of LAMOST DR5. He et al. (2022) successfully
distinguished between RGB and RC stars using the XGBoost algorithm (Chen &
Guestrin 2016), and used the SHAP interpretable model (Lundberg & Lee 2017) to
obtain the top features that the XGBoost selected. Shang et al. (2022)
presented a catalog of 6,917 CP1 and 1,652 CP2 new candidate sources using
XGBoost followed by the visual investigation, and listed the spectral features
for separating CP1 from CP2 using SHAP.

Before this work, the method of searching for carbon stars based on $Gaia$
spectra was mainly to obtain the matching giant spectra by random forest
classifier (Breiman 2001), and then to screen out carbon stars with strong CN
molecular bands by measuring the molecular band head strength (Gaia
Collaboration et al. 2023), actually most of which belong to the N-type.
Although this method can easily pick out the carbon stars with obvious
spectral features, it also introduces pollution and misses carbon stars that
exhibit relatively weak molecular bands . They are mixed with non-carbon stars
in the head strength map, which makes it difficult to distinguish using this
method. Our proposed method can be used as an improvement of the above method,
which can identify such carbon stars from the spectrum and give a
characteristic attribution explanation.

The organizational logic of this paper is as follows: In Section 2, we cover
the source of the data, the preparation of the training set, and data
processing methods. In Section 3, we introduce the proposed deep learning
classification model and important model parameters, as well as the SHAP
interpretability model and evaluation index respectively. In Section 4, we
show how well the model performs in the validation set, how it compares to
other classification models, demonstrating the effectiveness of model
interpretation and key features, and the newly discovered carbon star
candidates. In Section 5, We compare the results with others, analyze and
interpret the spectral features. Finally, a general summary is provided in
Section 6, and candidates are given in Section B.

##  2 Data

###  2.1 CSTAR sample of $Gaia$ DR3

$Gaia$ Data Release 3 ($Gaia$ DR3) has been released in 2022, and its mean
BP/RP spectral data was released in May of the same year, including four
folders on the archive: rvs_mean_spectrum, xp_continuous_mean_spectrum,
xp_sampled_mean_spectrum and xp_summary. The spectra we adopt are from
xp_sampled_mean_spectrum, which has 34,468,373 BP/RP externally calibrated
sampled mean spectra. All mean spectra were sampled to the same set of
absolute wavelength positions, viz. 343 values from 336 to 1020 nm with a step
of 2 nm, i.e. corresponding to a wavelength range of 3360 to 10200 Å in steps
of 20 Å. To standardize the data scale range, we simply normalized the raw
spectral data as

$x^{*}=\frac{x-x_{\min}}{x_{\max}-x_{\min}}\text{,}$ (1)

where $x$ represents the original spectrum, the $x_{\mathrm{min}}$ and
$x_{\mathrm{max}}$ respectively indicate the minimum and maximum value of the
$x$, and the normalised spectrum is denoted as $x^{*}$.

Min-max normalization is a simple data standardization method that maps the
data to a 0-1 range in equal proportion, thus changing the distribution of the
raw data. It only scales the spectrum without altering its original structural
representation. It also helps to speed up the convergence rate of the loss
function of the deep learning model, improving the efficiency and speed of the
model’s optimal solution, which may otherwise be difficult to converge.

Gaia Collaboration et al. (2023) (Hereafter referred to as C2023) released a
batch of golden samples, which include the golden sample of carbon stars from
‘CSTAR’. Their ESP-ELS module attempted to flag these suspected carbon stars.
The module is based on a random forest classifier trained on the synthetic BP
and RP spectra and a sample of Galactic carbon stars (Abia et al. 2020)
obtained from the $Gaia$ low-resolution spectra
($\mathrm{R}=\lambda/\delta\lambda\approx 25-100$; Carrasco et al. 2021). In
total, 386,936 targets received the ‘CSTAR’ tag. i.e., the potential carbon
star candidates. While most of these stars are M stars rather than carbon
stars, only smaller fraction of the sample exhibit significant
C${}_{\text{2}}$ and CN molecular bands. Most of the candidate carbon stars
have $G_{\mathrm{BP}}$ \- $G_{\mathrm{RP}}$ ¿ 2 mag and have colors consistent
with M stars (Gaia Collaboration et al. 2023).

As shown in Table 1, C2023 considered four molecular band heads in order to
further screening reliable carbon stars.

Table 1: The four molecular band head positions for identifying carbon stars.

Molecular band | Strength |  $\lambda_{1}$ [nm] |  $\lambda_{2}$ [nm] |  $\lambda_{3}$ [nm]
---|---|---|---|---
C${}_{\text{2}}$ | R${}_{\text{482.3}}$ | 462.2345 | 482.3455 | 505.3195
C${}_{\text{2}}$ | R${}_{\text{527.1}}$ | 505.3195 | 527.1080 | 546.5995
CN | R${}_{\text{773.3}}$ | 716.5865 | 773.2905 | 810.7805
CN | R${}_{\text{895.0}}$ | 806.8910 | 894.9855 | 936.6820

Finally, they selected 15,740 golden sample carbon stars based on the
strongest CN molecular band features. Most of these golden sample carbon stars
are AGB stars of which the enriched carbon was produced and then dredged up
from its interior region, i.e., the cold and bright C-N type carbon stars. The
calculation formula for the four molecular band head strengths is as follows:

$R_{\lambda_{2}}=\frac{f\left(\lambda_{2}\right)}{g_{\lambda_{1},\lambda_{3}}\left(\lambda_{2}\right)}\text{,}$
(2)

where $f\left(\lambda_{2}\right)$ is the flux measured at the top of the band
head of the molecular band, and $g_{\lambda_{1},\lambda_{3}}$ is the value
linearly interpolated between wavelengths $\lambda_{1}$ and $\lambda_{3}$
(Gaia Collaboration et al. 2023).

The reason C2023 did not select C${}_{\text{2}}$ molecular band strengths as
filtering criterion may be twofold: firstly, the C${}_{\text{2}}$ features in
the spectra are significantly weaker compared to CN features; secondly, more
than half of the golden sample carbon stars lack prominent C${}_{\text{2}}$
features in their spectra after visual inspection. Therefore, C${}_{\text{2}}$
features lack clear distinction, as can also be seen in the band head
strengths diagram of Figure 4.

To exploit these samples, The identification results of ‘CSTAR’ are saved in
the spectraltype_esphs field of the astrophysical_parameters table. The list
of ’source_id’ for golden sample of carbon stars is available in the
‘gold_sample_carbon_stars’ archive of the performance verification table.

###  2.2 Training data

After cross-matching the source_id of CSTAR with the spectral data of the
xp_sampled_mean_spectrum library in $Gaia$ DR3, we obtained 83,028 spectra.
These spectra will be used to verify the feasibility of our algorithm, which
we refer to as CSTAR_XPM (CSTAR sample with xp_sampled_mean spectra). Of
these, 8288 belong to the 15,740 golden sample of carbon stars and 74,740 are
non-golden sample. We respectively refer to them as CSTAR_XPM_G (for the
golden sample of CSTAR_XPM) and CSTAR_XPM_nonG (for the non-golden sample of
CSTAR_XPM).

Considering that the $T_{\mathrm{eff}}$ provided by $Gaia$ DR3 may not be
reliable, in Figure 1, we also plotted the $G_{\text{BP}}-G_{\text{RP}}$ and
$\log{g}$ diagrams of all CSTAR_XPM in order to get a rough estimate and a
basic idea of the physical properties and evolutionary state of these sources.
Since a large proportion of the sources are missing $T_{\mathrm{eff}}$ and
$\log{g}$, we plotted 3245 CSTAR_XPM_G sources, 47,290 CSTAR_XPM_nonG sources,
and 6763 random sample sources.

(a) $T_{\mathrm{eff}}$ and $\log{g}$ of CSTAR_XPM

(b) Random sample and CSTAR_XPM

(c) Gold and Non-gold from CSTAR_XPM

Figure 1: Panel (a): the $T_{\mathrm{eff}}$ and $\log{g}$ diagram of
CSTAR_XPM; Panel (b): The spatial location and density distribution of the
CSTAR_XPM and random sample in the $G_{\text{BP}}-G_{\text{RP}}$ and $\log{g}$
plane. The grey random sample is randomly selected from xp_sampled_mean. It
can be observed that these stars of CSTAR_XPM are all located in the giant
star branch, with most of them having 2 ¡ $G_{\mathrm{BP}}$ \-
$G_{\mathrm{RP}}$ ¡ 5 mag; The distribution of stars labelled as golden sample
carbon stars (red) and non-golden sample carbon stars (blue) in CSTAR is shown
in the right panel (c).

It is clear to see the CSTAR_XPM_nonG and CSTAR_XPM_G are spatially
intermingled in the distribution and are difficult to distinguish from each
other.

From the panel (a) and (c) of Figure 1, we noted that a few targets (19) have
a $\log{g}$ greater than 3 and $T_{\mathrm{eff}}$ hotter than 6000 K. However,
their $G_{\text{BP}}-G_{\text{RP}}$ is around 3, and their spectra are typical
of AGB carbon stars (showing strong CN bands), which is quite unexpected.
Although the possibility of them being extrinsic carbon stars (such as dwarf
carbon stars) cannot be ruled out, C2023 mentioned that the $T_{\mathrm{eff}}$
of these targets tends to be overestimated. Therefore, it is more likely that
these are due to measurement errors.

In order to further ensure the purity of the carbon star sample, we visually
checked all the CSTAR_XPM_G, finding that some spectra exhibit very
significant contamination, which may be due to the signal-to-noise problem, or
exhibit insignificant CN characteristics. These spectra would have troubled
the purity of our training set. Therefore, we need to eliminate the
contaminated spectra. We show the BP and RP spectra of a standard carbon star
and some of the deleted problematic stars in Figure 2. After screening we
selected 8176 carbon stars as the positive training sample, which showed
significant CN molecular band characteristics.

Figure 2: These images are normalized XP spectra of CSTAR_XPM_G. The upper
panel shows a standard carbon star and the range of four molecular bands
(purple areas) marked, where the red dotted line marks the location of the top
of the band head of the molecular band. The middle panel shows three
problematic spectra and the bottom panel shows three spectra with very faint
molecular band features.

We then randomly selected 9,000 spectra out of 74,740 CSTAR_XPM_nonG as
”negative sample”, the majority of which are M-type stars (Gaia Collaboration
et al. 2023). First, we made sure that none of the 9,000 spectra are included
in several known common carbon star lists (Solar Neighbourhood Carbon stars:
Abia et al. 2020; Galactic Carbon stars: Alksnis et al. 2001; LMC Carbon
stars: Kontizas et al. 2001; Carbon stars in the SMC: Morgan & Hatzidimitriou
1995), and that none of them were labeled C* by SIMBAD or listed as ‘Carbon’
by LAMOST’s pipeline. Second, here we must mention our training method. We
randomly trained several models based on the above dataset. We identified
false positive (FP) spectra resulting from obvious model’s misclassification,
where the true category was 0 (other stars), but the model predicted category
1 (carbon star). We iteratively added these spectra to our ”negative sample”
dataset, finding that this data enhancement method significantly improved the
model’s generalization capability. Finally, during the visual examination of
the ”negative sample”, we found that many spectra exhibited relatively obvious
and prominent molecular bands, particularly two strong CN molecular bands.
Therefore, in order to further ensure the purity of the ”negative sample” and
ensure that the model could learn the correct features, we calculated the band
head strength for all spectra according to the method provided in C2023. We
screened out all the suspicious potential carbon stars under the two
conditions of R${}_{\text{773.3}}\geq 1.03$ and R${}_{\text{895.0}}\geq 1.10$,
which are almost the weakest CN molecular band strength conditions. We
obtained 2255 spectra in total, and then carefully identified the weak CN band
head by eyes to exclude the interference of possible M and K type giant stars.
The reason is that the absorption of M and K giants at the TiO5 (7126-7135 Å;
Li et al. 2023) position is more abundant, and may exhibit a
pseudo-R${}_{\text{773.3}}$ molecular band strength. Through this checking
process, we eliminate 885 objects that exhibit weaker CN molecular bands
relative to the cold and bright golden sample of carbon stars. We showed six
representative spectra in Figure 3.

Figure 3: Six representative spectra were selected from the 885 objects
removed. The two prominent CN bands (purple areas) are located at the
7165-8108 Å and 8160-9367 Å ranges. The red dashed lines mark the top of the
band head. The spectra are normalised and offset from one another by $k\times
1.0+0.2$ for clarity (where $k$ is an integer that varies from 0 to 6 from the
bottom to the top spectrum).

Note that these spectra only show weak CN molecular bands, which does not mean
that they are all potential carbon stars, we just try to eliminate carbon star
contamination as much as possible to obtain pure giant star sample. We also
searched and excluded the problematic spectra with condition
R${}_{\text{773.3}}$ ¿ 1.10 or R${}_{\text{895.0}}$ ¿ 1.30, as some of the
problematic spectra may have large R${}_{\text{773.3}}$ or
R${}_{\text{895.0}}$. Then we screened the obtained spectra for three times to
ensure that all spectra did not show any carbon characteristics and then added
some FP spectra that apparently misclassified by the model mentioned above.
Finally we kept 8556 training negative sample, which can be confidently
identified as non-carbon stars. We plotted the molecular band head strength
distribution of the sample in Figure 4.

Figure 4: Slightly different from C2023, for the 74,740 CSTAR_XPM_nonG (blue)
flagged by ESP-ELS; 8176 selected training positive sample stars (red) from
CSTAR_XPM_G; 8556 selected training negative sample stars (green) from
CSTAR_XPM_nonG and 885 deleted stars (yellow).

All the samples we used were summarised in Table 2. The Venn diagram in Figure
5 depicts the relationship between these samples.

Table 2: The sample sets used to validate our algorithm.

Sample name | Number of spectra | Sample definition
---|---|---
CSTAR_XPM | 83,028 | CSTAR($a$) sample with xp_sampled_mean spectra
CSTAR_XPM_G | 8,288 | Golden sample($b$) of CSTAR_XPM, the vast majority of which are carbon stars
CSTAR_XPM_nonG | 74,740 | Non-golden sample of CSTAR_XPM, most of which are non-carbon stars
Positive sample | 8176 | Carbon stars selected from CSTAR_XPM_G, used for model training
Negative sample | 8556 | Non-carbon stars selected from CSTAR_XPM_nonG, used for model training

  * •

Notes.

  * •

($a$) 386,936 candidate carbon stars received the “CSTAR” tag (Gaia
Collaboration et al. 2023).

  * •

($b$) 15,740 golden sample carbon stars from CSTAR.

Figure 5: Venn diagram of the relationship of all the samples in CSTAR (green
frame). The spectra we used is CSTAR_XPM (purple frame, left in the figure).
CSTAR_XPM_G is a cross between CSTAR_XPM and golden sample of carbon stars
(purplish-red frame, right in the figure), and CSTAR_XPM_nonG is the
complement of CSTAR_XPM and CSTAR_XPM_G. Negative sample (yellow) is selected
from CSTAR_XPM_nonG, and positive sample (blue) is selected from CSTAR_XPM_G.
Grey frame indicate the range of possible carbon stars.

We finally got 16,732 stars as training sample, the total number of which is
approximately 20.15% of the total number of CSTAR_XPM sample. We will use this
training sample to verify the accuracy of our classification model and the
explanatory performance of the interpretable model. We speculate that there
may be some additional carbon stars in CSTAR_XPM_nonG, so we will apply the
trained model to this sample.

##  3 Method

###  3.1 GaiaNet

Convolutional Neural Network (CNN; LeCun et al. 1998) is proposed by the
concept of perceptual field in biology. By setting several convolutional
kernels of different sizes, the model can effectively capture a wider range of
global and local information, thereby enhancing the feature extraction of the
input signal. The neural network model proposed in this work is tentatively
named ‘GaiaNet’.

####  3.1.1 The structure of GaiaNet

The model is a ”light-weight” one-dimensional convolutional neural network, an
improvement based on the TextCNN model proposed in Kim (2014). The main reason
for choosing this algorithm is that the flexible-sized convolutional kernel of
the 1D CNN can effectively capture and extract detail and overall key features
when sliding the convolution in one direction, such a working characteristic
makes the model more convincing and interpretable.

There is experimental evidence that extending the depth of a 1D convolutional
model can effectively enhance the model’s fitting and feature learning
capabilities (Chen 2015). A parallel CNN structure consisting of multiple sets
of convolutional kernels of different sizes can effectively capture key
features of different sizes and achieve similar effects to ensemble learning.
Both the serial deep structure and the wider parallel structure scheme can
effectively capture key features. We have drawn on these advantages, which are
well reflected in the improvement of our model. We first modified the size of
the convolutional kernel in the input layer so that it can receive one-
dimensional data, i.e., data of shape 1 × 343. We added regularization units
consisting of BN (Batch normalization; Ioffe & Szegedy 2015) and Dropout. Both
methods are effective in improving the generalization of the model. For the
data in each batch, the BN layer first normalizes them. The formula for BN is
as following:

$\hat{x}=\frac{x-\mu_{B}}{\sqrt{\sigma_{B}^{2}+\epsilon}}\text{,}$ (3)

where $x$ denotes the input data, $\mu_{B}$ and $\sigma_{B}^{2}$ are the mean
and variance of the datas within that batch, respectively, and $\epsilon$ is a
very small constant to prevent division by zero errors. The normalized data
are then linearly transformed by a set of learnable parameters ($\gamma$ and
$\beta$), which can be optimized by a back-propagation algorithm to obtain the
final output $y$:

$y=\gamma\hat{x}+\beta\text{.}$ (4)

The BN layer avoids drastic changes in the input data distribution caused by
parameter updates in the preceding model layer and ‘pulls’ the data back into
a stable distribution space. Experimental results indicate that it stabilizes
the network, effectively accelerates the training and convergence processes,
and enhances the classification accuracy of the model.

Dropout is a commonly used regularization method to avoid overfitting deep
network models (Hinton et al. 2012). Dropout reduces the risk of model
overfitting by randomly discarding the outputs of neurons to reduce the
complex dependencies within the network during model training. It is
implemented by temporarily erasing the output of a neuron with probability
$p$:

$y_{i}=\begin{cases}\frac{x_{i}}{1-p},&\text{with probability }1-p\\\
0,&\text{with probability }p\end{cases}$ (5)

where $x_{i}$ denotes the original output of neuron $i$ and $y_{i}$ denotes
the output of neuron $i$ after Dropout. During the training phase, each neuron
is randomly selected and turned off with a probability of $p$ and retained
with a probability of $1-p$. The retained neuron needs to be scaled by
dividing its output value by $1-p$ to keep the expectation of the output value
constant. We add Dropout before the final fully connected layer of the model
to prevent overfitting and to improve the generalization of the model.

Our model is made up of 16 parallel network structures, each with different-
sized convolutional input layers to be able to better capture and extract key
features in different fields of view. The same deep hidden layer network
structure is then used, which consists of three convolutional layers, two
pooling layers, a maximum pooling layer, and a average pooling layer. The
deeper hidden layer structure helps to enhance the model’s non-linear fitting
capability, allowing for the extraction of key non-linear features. The
overall architecture of the model is derived from ‘textCNN’, which has been
improved by several iterations. We referenced the NIN network by introducing a
1×1 convolution, and a global average pooling operation (Lin et al. 2013). The
Inception Module network structure from GoogLeNet is borrowed, i.e., a
structure consisting of multiple parallel and smaller dense convolutions with
pooling approximating sparse connections (Szegedy et al. 2015). The flow chart
of the model is as Figure 6.

Figure 6:  The convolution kernels of different sizes are used to extract
features with different field-of-view sizes. The arrows point to the tensor of
the feature map after convolution, pooling, or linear layer operations. The
global maxPooling layer downsamples the feature map to achieve a compressed
feature map size, and captures the most globally sensitive features in the
feature map, stitches them together to obtain a tensor, which then undergoes
global average pooling to share each feature value in the tensor. It then goes
through a fully connected layer to summarise the features to obtain a high-
level abstract representation of the features. Finally, a Sigmoid layer allows
the output of the neural network to be transformed into a probability
distribution and makes the model predictions easier to understand and
interpret. Below the image are labeled the respective convolutions, pooling,
and the shape of the feature map.

The 1×1 convolution (Lin et al. 2013) can further deepen the network,
effectively reduce the number of parameters in the overall convolution layer
of the model, enhance the model’s non-linear fitting ability, allow for
flexible dimensioning up and down of the feature map, and enable cross-channel
information interaction and integration. Global averaging pooling replaces the
original practice of stitching the maximum pooled feature maps directly into
the fully connected layer. These feature maps are passed through the global
average pooling layer and then concatenated to the linear layer to give the
corresponding categories. The purpose of this is to share information from the
multi-channel feature maps so that they can all contribute to the final
result, which we believe makes more sense and is one of the main innovations
of this model. The formula is as below:

$y_{i}=\frac{1}{H}\sum_{j=1}^{H}x_{i,j}\quad i\in[1,n]\text{.}$ (6)

The global average pooling is handled very simply, where $x_{i,j}$ denotes the
feature map value at the $j$-th row position of the $i$-th input tensor and
$H$ denotes the height of the tensor. The global averaging pooling layer
therefore calculates the average of all pixel points within a single tensor
activation feature map. The overall structure of the model is shown in the
Figure 7.

Figure 7: The input data is initially passed through convolutional input
layers of different sizes and a depth hidden layer, and then through a global
average pooling layer before being stitched together, and finally through a
fully connected layer to obtain a probabilistic output of [0,1] using the
Sigmoid layer.

We conducted extensive comparative experiments, and the results indicate that
balancing the depth and width of the network, combined with the above
improvements, can significantly enhance the model’s performance.

####  3.1.2 Important model parameters

The output layer of our model finally passes through a Sigmoid activation
function (Lippmann 1987), which is an ‘S’-type function that converts the real
values of the output into a 0-1 probability distribution to represent the
confidence level of a positive case. We usually set 0.5 as the threshold
value, cases greater than the threshold being classified as positive, and
cases less than or equal to the threshold being classified as negative. It is
in line with the 0,1 label defined when we process secondary classification
tasks. We use the BCELoss (Binary Cross Entropy Loss; Krogh & Hertz 1991)
function to measure the difference between the model output and the true
label, which can be interpreted as a maximum likelihood estimate, i.e.
maximizing the probability of the observed data and finally converging the
loss by an optimizer. It is a natural pairing with the Sigmoid activation
function when dealing with dichotomous problems, which makes the model output
easier to interpret and understand. The formula for calculating the loss of
the output from the Sigmoid activation function using BCELoss is as below:

$\mathcal{L}=-\frac{1}{N}\sum_{i=1}^{N}\left[y_{i}\log(\hat{y}_{i})+(1-y_{i})\log(1-\hat{y}_{i})\right]\text{,}$
(7)

where $N$ represents the number of input cases in a batch, $y_{i}$ denotes the
true label of the $i$th case, and $\hat{y}_{i}$ signifies the predicted
probability (model output value) of the $i$th case. The $p$-parameter of the
Dropout layer is also an important parameter, which refers to the proportion
of randomly ‘lost’ neurons. The output of our several parallel neural network
structures is then spliced together to obtain a vector that is finally
discarded with a given $p$-ratio after the Dropout layer. This approach
promotes the model’s generalization, effectively avoids model overfitting, and
achieves a good integrated learning effect.

In addition, it is important to balance the size of the convolutional kernels.
As larger kernels can effectively expand the model’s receptive field and
capture features in a broader context, while smaller kernels excel at
capturing local details. Therefore, we need to strike a good balance between
the whole and the local depending on the actual task.

Furthermore, the batch_size is the number of input cases in a batch used to
update the model parameters for each iteration. A larger batch_size can speed
up convergence by reducing the number of parameter updates. A smaller
batch_size can improve the generalization ability of the model. By randomly
selecting small batches of case for parameter updates introduces greater
variability in each batch, facilitating the model’s search for a global
optimal solution without overfitting the training set. The selection of
batch_size should be comprehensively considered based on factors such as the
size of the dataset, the nature of the task, memory and computing power.

The number of epochs for training is a critical hyperparameter that requires
careful consideration. Excessive epochs can lead to overfitting, where the
model performs well on the training set with good loss convergence but poorly
on the validation set due to decreasing accuracy rates with increasing
training epochs. On the other hand, insufficient epochs may result in
underfitting, causing the model to inadequately learn key information and
perform poorly on the actual test set. Therefore, it is crucial to select an
reasonable number of epochs based on factors such as model capacity, data
volume, and task complexity.

Finally, there are several parameters that play a crucial role in the SGD
(Stochastic Gradient Descent) optimization algorithm during the training phase
of the model. The $lr$ (learning rate) is a factor that controls the step size
of the model weights update, which represents the size of the updated weights
at each iteration. A large learning rate may lead to fast convergence of the
model during training, but may also lead to the model skipping the optimal
solution. A smaller learning rate may result in slower convergence of the
model, but it is more likely to achieve a better optimal solution. In our
practical tests, we opted for a smaller learning rate to avoid missing the
global optimal solution. During the SGD training shown in Figure 8, the size
of each step is fixed, but with the introduction of the momentum learning
algorithm, the movement of each step depends not only on the magnitude of the
current gradient but also on the accumulation of past velocities. Momentum
uses historical gradient information to adjust the direction and speed of
parameter updates (Polyak 1964), thus accelerating the convergence process of
the SGD gradient descent and reducing oscillations during the gradient update
process. It can be seen as an inertia introduced to leverage the previous
gradient information in parameter updates. The core idea of momentum is to
introduce an accumulated gradient history variable, which is similar to
momentum in physics and records the direction and velocity of the previous
gradient’s motion. During each iteration of the update, the momentum algorithm
considers not only the current gradient but also the trends of previous
gradients. This allows the parameters to be updated with a certain ‘momentum’
in the direction of the gradient, resulting in faster crossing of flat areas
and less invalid oscillation. The parameter update formula in the momentum
algorithm is shown below:

$v_{t+1}=\mu v_{t}-\varepsilon\nabla f\left(\theta_{t}\right)\text{,}$ (8)
$\theta_{t+1}=\theta_{t}+v_{t+1}\text{,}$ (9)

where $\varepsilon$ ¿ 0 represents the learning rate, $\mu\in[0,1]$ is the
momentum parameter, and $f\left(\theta_{t}\right)$ is the gradient at
$\theta_{t}$ (Sutskever et al. 2013). We were surprised to find in our
practical tests that it was excellent for improving model accuracy and further
converging model loss.

Figure 8: As shown in the figure, the contour line depicts a binary loss
function. The left panel depicts the loss update path of a normal SGD, while
the right panel shows the loss update path with the addition of momentum, we
can see that the update path of the SGD with the use of momentum is smoother
and more stable, and it is also easier to reach the global optimal solution
location.

Weight decay is a regularisation technique that reduces the complexity of the
model by reducing the size of the weights. It applies an L2 regularisation
penalty to the weight parameters in the loss function to prevent overfitting,
so that larger weights are penalized more during training, thus encouraging
the model to learn a simpler weight distribution.

###  3.2 SHAP

SHAP (SHapley Additive exPlanations) is an interpretable model based on the
game theory principles approach to explain the output of any machine learning
model. It connects optimal credit allocation with local explanations using the
classic Shapley values from game theory and their related extensions. SHAP can
give the extent to which each input feature affects the output of the model
and offer a comprehensive and interpretable analysis of the anticipated
results of the entire model (Lundberg & Lee 2017). The Shapley value is a
concept in game theory that measures the contribution of each player to the
cooperative game. In machine learning, we can consider different features as
”players” and the combination of their values as a ”cooperative game”. Using
the principle of equivalent assignment of Shapley values, the contribution or
importance of each input feature to the overall model output can be
calculated, and each feature can be assigned a relative contribution value.
And SHAP extends Shapley values to deep learning models.

The prediction value of the model is interpreted as the sum of the importance
of each input feature. The relation between the predicted value and importance
scores is performed below:

$\centering
f(x)=g\left(x^{\prime}\right)=\phi_{0}+\sum_{i=1}^{M}\phi(x_{i}^{\prime})\text{.}\@add@centering$
(10)

Specifically, for a deep learning model $f:X\to Y$, where $X$ denotes the
input features and $Y$ denotes the class or real value of the output, $M$ is
the number of features. The explanation model $g\left(x^{\prime}\right)$
matches the original model $f(x)$ when $x=h_{x}(x^{\prime})$, where
$\phi_{0}=f(h_{x}(0))$ represents the model output with all simplified inputs
toggled off (i.e. missing). In a word, SHAP defines the Shapley value for an
input $x_{i}\in X$ as $\phi(x_{i}^{\prime})$, which represents the marginal
contribution of feature $x_{i}$ to the model prediction.

The Shapley value of each feature represents the degree of its influence on
the predicted outcome of a single input case. The calculation of the Shapley
value is based on combinatorial game theory, which takes into account the
interactions between each feature and other features and assigns the final
predicted outcome to each feature. For a specific case, the calculation of the
Shapley value $\phi(x_{i}^{\prime})$ of feature $x_{i}$ is performed below:

$\centering\phi(x_{i}^{\prime})=\sum_{S\subseteq X\setminus
x_{i}^{\prime}}\frac{|S|!(M-|S|-1)!}{M!}[f(x_{S}\cup{x_{i}^{\prime}})-f(x_{S})]\text{,}\@add@centering$
(11)

where $S\subseteq X\setminus x_{i}^{\prime}$ is any possible subset of all
input features excluding $x_{i}$. $f(x_{S}\cup{x_{i}^{\prime}})$ represents
the output obtained by adding feature $x_{i}$ to the subset $S$. $f(x_{S})$
represents the output obtained by using the subset $S$ for prediction, and the
result of subtracting the two values represents the effect of the addition of
feature $x_{i}$ on the model output. And the ${|S|!(M-|S|-1)!}\setminus M!$
represents the weight. The overall equation represents the calculated
contribution for each feature subset $S$. This is the contribution from adding
feature $x_{i}$ to the feature subset $S$. Finally, the contribution of each
feature subset $S$ is summed to obtain the Shapley value of the feature
$x_{i}$.

This method allows us to compare the Shapley values of different features to
determine how much they affect the output of the model. It has been flexibly
applied to various deep learning models, including convolutional neural
networks (CNN), recurrent neural networks (RNN; Werbos 1988), and others. We
will use SHAP in our model to interpret the key features that affect the model
output. In addition, SHAP provides several visualization tools and methods,
which will help us understand the model better.

###  3.3 Evaluation index

In addition to the accuracy criterion, we also introduced more refined
evaluation metrics such as recall and precision to effectively measure the
model’s data mining performance during dataset validation. We define TP (true
positive) as the number of cases whose categories are true and predicted
categories are also positive, FN (false negative) as the number of cases whose
categories are true but predicted categories are negative, FP (false positive)
as the number of cases whose categories are false but predicted categories are
positive, and TN (true negative) as the number of cases whose categories are
false and predicted categories are also negative.

Accuracy refers to the proportion of the number of accurate cases classified
by all categories to the total number of cases. We can calculate it as below:

$Accuracy=\frac{TP+TN}{TP+TN+FP+FN}\text{.}$ (12)

Precision refers to the proportion of the number of correctly predicted
positive cases to the total number of predicted positive cases. It is
calculated as below:

$Precision=\frac{TP}{TP+FP}\text{.}$ (13)

Recall refers to the proportion of correctly predicted positive cases relative
to all actual positive cases. It is calculated as below:

$Recall=\frac{TP}{TP+FN}\text{.}$ (14)

We also use confusion matrices to see the performance of the model more
intuitively, and confusion matrices provide a good indication of the
distribution of the various cases. For $Gaia$’s Carbon Star data mining work,
the idea is to have as many TP cases as possible and as few FN cases as
possible, i.e., as large a recall rate as possible, so we use this as our main
validation set measure.

###  3.4 Resampling technique

The resampling technique (Efron 1979) is an important data processing method
used to solve the problem of imbalance in the categories of a data set. For
example, when the number of positive cases is slightly less than the number of
negative cases, resampling techniques can be used to increase the number of
positive cases. This is done by randomly copying or generating new positive
cases, thus making the number of positive and negative cases consistent and
balanced. This can effectively prevent the models from ”collapsing” into
categories with large case sizes during training.

##  4 Results

After validation with several sets of parameters in Section 7, we fixed the
training parameters of the model as follows: $p=0.4$, $lr=0.005$,
momentum$=$0.9, weight_decay$=$0.001.

###  4.1 Analysis and experiments of Dataset

Almost all C-N type carbon stars are in the AGB evolution stage, of which the
carbon abundances in their atmosphere were produced in the stellar interiors,
and were brought to the surface by the third dredge-up process. Almost all of
them are intrinsic, thermally pulsing AGB stars, i.e., TP-AGB stars with large
luminosities (Abia et al. 2002). This kind of carbon stars typically exhibit
stronger C${}_{\text{2}}$ and CN molecular bands than usual stars with
temperatures below 3800 K (Gaia Collaboration et al. 2023). Of course, there
are also some kinds of carbon stars with relatively high surface temperatures,
whose spectral types can be K or even G, and they are often extrinsic carbon
stars. However, there is a large number of giant stars with the spectral type
G, K, and M, which have similar colors, surface temperatures, and luminosities
to the carbon stars, making it impossible to distinguish between the carbon
stars and other giants just from the HR diagrams, but now possible through
$Gaia$’s spectroscopy. Lebzelter et al. (2023) and Sanders & Matsunaga (2023)
have demonstrated its potential for separating C-rich and O-rich stars.

We found that almost all of the CSTAR_XPM_nonG are M giants through the
$G_{\text{BP}}-G_{\text{RP}}$ and $\log{g}$ diagram (see Figure 1), and their
distribution overlaps with the region where the golden carbon stars located.
In order to obtain reliable luminosity and color parameters, we corrected the
photometric data with 3D interstellar extinction (Green et al. 2019), and then
a color-magnitude diagram is shown in Figure 9.

Figure 9: Distribution of CSTAR_XPM and training sample on $(M_{G})_{0}$ vs.
$(G_{\text{BP}}-G_{\text{RP}})_{0}$ diagram. All the samples have been
corrected for 3D extinction. Panel (a): the distribution of the CSTAR sample.
Panel (b): the distribution of the training sample. Panel (c): the
distribution of carbon stars from the training sample. Panel (d): the
distribution of Non-carbon stars from the training sample.

After the extinction correction, we calculate the absolute magnitudes
$(M_{G})_{0}$ and the intrinsic color index of
$(G_{\text{BP}}-G_{\text{RP}})_{0}$ for 4002 positive sample stars, 6,761
negative sample stars, and 64,761 stars of CSTAR_XPM sample, respectively.
Most of the CSTAR_XPM stars have 1 ¡ $(G_{\text{BP}}-G_{\text{RP}})_{0}$ ¡ 3
and $(M_{G})_{0}$ ¡ 2 mag, and most of the carbon stars have 2 ¡
$(G_{\text{BP}}-G_{\text{RP}})_{0}$ ¡ 4 and $(M_{G})_{0}$ ¡ -1 mag. This
suggests that the CSTAR_XPM_G sample may consist mostly of brighter and cooler
C-N type carbon stars. From the corrected color-magnitude diagram, we also
found that about half of the negative sample stars have a
$(G_{\text{BP}}-G_{\text{RP}})_{0}$ $\geq 1.6$, while the
$(G_{\text{BP}}-G_{\text{RP}})_{0}$ of M giants should usually be larger than
1.6. We estimated that there are about half of the G and K giants in the
CSTAR_XPM_nonG, which is corroborated by our later cross-matched results with
the LAMOST DR10 spectroscopic data in Section 5.1. Hence, the CSTAR sample
also includes a large mixture of other types of stars, rather than
overwhelmingly M stars. In addition, a portion of the CSTAR_XPM_G has
$(G_{\text{BP}}-G_{\text{RP}})_{0}$ ¡ 2 and relatively lower luminosities.
From Figure 11 of Li et al. (2024), it can be inferred that they may be
extrinsic warm carbon stars with spectral types of G or K like Ba stars. It
can be seen that there is a slightly clearer distinction between positive
sample and negative sample in the matted color magnitudes, but some of the
stars are still mixed together. However, elemental abundances are also
reflected in the spectra. Specifically, the carbon stars are characterized by
the strong absorption lines of carbon molecules in their spectra. This offers
the possibility of a spectroscopic data-driven approach to distinguish the two
directly from the stellar spectra. C2023 proposed screening carbon stars by
two CN molecular band head strengths (R${}_{\text{773.3}}$ and
R${}_{\text{895.0}}$), which could be partially contaminating through our
inspection, and would result in the loss of some potential carbon stars.

As mentioned in Section 2.2, we carefully selected a negative sample of 8556
spectra from CSTAR_XPM_nonG and a positive sample of 8176 spectra from
CSTAR_XPM_G. We equalized the number of positive and negative samples by
randomly adding negative spectra using the resampling techniques. So we
composed a training dataset after fine visual inspection and machine ranking
to try to make the better distinction in a deep learning way and expect
meaningful interpretable analysis.

####  4.1.1 Model validation

We randomly selected 70 percent of the datasets to train the classification
model, and the remaining 30% was used as the validation set to measure the
model’s performance as the epochs increased.

The best performance of the model occurs when the epochs are 9, 16, 17, 18,
and 19. The final recall rate in the validation dataset is 100 percent, the
accuracy is also 100 percent. The confusion matrix is shown in Figure 10,
showing that our model has extremely high accuracy and reliability for
classifying all stars in the validation dataset. Our model can easily
distinguish the spectra of cold and bright carbon stars with strong carbon
molecular bands from the spectra of giant stars of G, K, M types. We will
provide necessary comparisons and explanations in Section 4.2 .

Figure 10: Label 0 represents a non-carbon star and 1 represents a carbon
star. The abscissa represents the predicted category and the ordinate
represents the actual category.

We also plot the change in the correctness of the training and validation sets
and the number of FP, and FN spectra as the training epochs increase, as shown
in Figure 11. The best performance on the validation set (30%) is 100%
accurate, with excellent convergence and no overfitting of the model.

(a) The accuracy of the training set and the validation set

(b) FN and FP

Figure 11: Panel (a): the performance of the model on the training set and the
validation set changes with the number of epochs; Panel (b): the number of
misclassified spectra of the model varies with the number of epochs.

####  4.1.2 Comparison with other algorithms

In the past decades, many machine learning includes integrated learning
methods (KNN: Cover & Hart 1967; SVM: Cortes & Vapnik 1995; Random Forest:
Breiman 2001; XGBoost: Chen & Guestrin 2016 etc.) have been successfully
applied to stellar spectral classification and physical parameters estimation.
In contrast, convolutional neural networks are well known to have been used
with great success in the field of image recognition, so we see the potential
for the use of one-dimensional convolutional deep learning networks in
spectral recognition.

In Table 3, we compiled a comparison of the accuracy (%), recall, and error
rate (%) of several common machine learning classification models with our
model on the performance of the validation set.

Table 3: For each model, the data was split into 70% for training, and 30% for
validation.

Method | Accuracy (%) | error (%) | Recall
---|---|---|---
Our Model1 | 100 | 0.0 | 1.0
XGBoost | 99.82 | 0.18 | 0.99794
Random Forest | 99.80 | 0.2 | 0.99678
KNN | 99.70 | 0.3 | 0.99518
SVM | 99.62 | 0.38 | 0.99282

  * •

Notes.
1 The best performance in the validation set is selected when the epochs are
9, 16, 17, 18, and 19.

Recall is a very important evaluation index that effectively measures the rate
of retaining positive sample. Our model retains all stars of positive sample
and demonstrates approximately a 0.43% improvement over other algorithms,
resulting in an overall accuracy increase of about 0.3%. The reason for the
error rate of 0 in the validation set is likely to be the excellent overall
architecture of our model, which is mainly reflected in the natural advantages
of the convolutional model in feature extraction, as well as some necessary
innovative improvements. On the other hand, we have conducted ”expert”
screening and purification of the training data to exclude erroneous sprctra
and ensure the high purity and reliability of the training data, combined with
some necessary data enhancement methods to ensure that the model can learn the
real key information. Of course, the spectra of golden sample almost all have
strong molecular bands, making them easy to distinguish from the spectra in
CSTAR_XPM_nonG.

###  4.2 Model interpretation and key features

The SHAP value is used to quantify the effect of features on the model output.
We consider each wavelength data point in a spectrum as a feature, and the
corresponding flux to each wavelength as the feature value. The SHAP
interpretable model can explain the output of our deep learning model by
calculating a SHAP value for each feature of a given spectrum, which indicates
the importance impact or contribution of the feature to the model output (the
confidence of a carbon star as mentioned in subsection 3.1.2). A positive SHAP
value indicates that the feature contributes positively to the model’s
prediction of the spectrum being a carbon star, while a negative SHAP value
indicates that the feature contributes negatively to the model’s prediction of
the spectrum being a carbon star. The larger the absolute value of SHAP, the
greater the impact of the feature on the model output. For a spectrum, the sum
of the SHAP values of each feature can be used to approximate the output of
the classification model. We can also calculate the distribution of the global
most important features based on the statistics of SHAP values from a large
sample.

To get an overview of which features are most important and intuitively
understand the influence of them on the model prediction results, we randomly
selected 1024 carbon stars of positive sample and 1024 non-carbon stars of
negative sample as background spectra of the deep explainer interpreter of
SHAP, calculated the SHAP values of these spectra, and plotted the beeswarm in
Figure 12. The top 30 most important features were shown in this figure, which
were ranking by the sum of absolute SHAP value magnitudes of them over all
2048 spectra, and uses SHAP values to show the distribution of the impacts
each feature has on the model output. For a particular feature, SHAP first
normalize its values across all samples, the normalized feature values are
then mapped onto a color gradient bar. The feature point will be assigned
different colors (from blue to red), the color blue represents the feature
point value is relatively small, and red means that feature point value is
relatively large.

Figure 12: The beeswarm visualizing the link between feature values and model
interpretation. Each point in the beeswarm plot represents a feature point
from a spectrum. The horizontal axis describes the SHAP values corresponding
to each feature point, and the vertical axis is the wavelengths of the top 30
most important features for identifying carbon stars (From top to bottom, it
represents the decreasing contribution of features to the model prediction).
The color of each feature point represents the size of the feature value.

Therefore, the beeswarm plot can be very intuitive reflect the influence of
the feature on the output results under different values. As shown in the
figure, we can see that the wavelength range 8900-9000 Å is the location of
the peak arising from the molecular band head. The larger the values of these
feature points are, the more positively the model predicts carbon stars. Near
8100 Å and 9300 Å are the locations of the troughs arising from the molecular
absorption band. The smaller these feature values are, the more negatively the
model predicts carbon stars. It is evident that the most striking
distinguishing feature is the presence of these two troughs.

To better illustrate the differences in key feature contributions within a
single spectrum, the background spectrum for calculating SHAP values can be
chosen as the smooth pseudo-continuum after polynomial fitting using the
approach of Zhang et al. (2021). It serves as a baseline and reference for
individually calculating the SHAP value of each feature in each spectrum. We
believe that contrasting with such a baseline can effectively highlight the
trough and peak features in each spectrum. Subsequently, we can generate an
intuitive and clear feature thermal distribution map for the spectrum shown in
Figure 13, which facilitates observing the distribution locations of the key
features concerned by our model.

(a) Carbon star

(b) Non-carbon star

Figure 13: The baseline is a background pseudo-continuum used for computing
SHAP values in order to better highlight spectral features. The color
represents the size of the SHAP value, and a redder feature indicates its
stronger positive contribution of the SHAP model’s interpretation of the
carbon stars. The bluer the color, the more strongly the feature contributes
negatively to the model interpretation of the carbon star. The sum of the SHAP
values can be considered as the approximation of the model output, and the
larger the sum of the SHAP values, the more inclined to consider the spectrum
as a carbon star.

However, if we want to get a more statistically significant result and reach a
definitive conclusion, we need statistical information on the SHAP values of a
large enough sample, which is more revealing than above. For this purpose, we
calculate the average SHAP values of each feature for all positive and
negative samples. As shown in Figure 14, we can see that the interpretation
ranges of positive features for positive sample are consistent with the CN
features. It is mainly distributed in the ranges 9140-9440 Å, 7900-8360 Å,
8820-9060 Å, and 7600-7820 Å, showing that relatively large SHAP values were
assigned to both peak and trough ranges in the CN molecular bands. However,
compared with the positive carbon star sample, the average SHAP values in
these regions of the negative sample are much smaller. In addition, the
negative sample showed a very strong negative contribution in the range
8420-8900 Å, which is the most significant distinguishing feature of this
sample. We note that the average SHAP values of the positive sample is also
negative in the 8420-8620 Å region. We suspect that this is related to the
absorption feature of Ca II (8484-8662 Å), which is considered by the model to
be a negative feature. We can see that the location of the model’s main
concern is the strong CN molecular band feature at the red end of the
spectrum, which is consistent with the actual interpretation in Section 20.
Next, we will showcase these key features.

Figure 14: The average SHAP values per feature is shown using a wavelength bin
size of 20 Å. Red and blue represent the positive and negative SHAP values
respectively. As a contrast, the average spectra of positive and negative
samples are overlaid as green scatter points in the graph. In the upper panel,
the two prominent CN bands are highlighted with gray shading. The dashed lines
mark the positions of the tops of the two CN molecular band heads.

The sum of SHAP values in Figure 14 represents the sum of the average SHAP
values across the feature dimensions of positive or negative sample. For the
positive sample, the sum of SHAP values is 0.99894. We summarize the main
feature regions of the golden sample of carbon stars, which will serve as the
basis for our subsequent identification work. We calculate the sum of the SHAP
values over 5 main intervals listed in Table 4.

Table 4: The wavelength range and SHAP value score of the 5 most important features from positive sample. Number | Feature area range | Sum of SHAP values
---|---|---
Feature 1 | 9140-9560 Å | 0.46825
Feature 2 | 7900-8360 Å | 0.37919
Feature 3 | 8820-9060 Å | 0.11090
Feature 4 | 7600-7820 Å | 0.05297
Feature 5 | 6780-7460 Å | 0.02585

Among them, the range 6780-7460 Å includes the peak characteristic of
6780-6920 Å and the absorption trough characteristic of 7000-7460 Å.

The reason our model only identifies the prominent CN features in the red end
of the positive sample spectra, while almost failing to detect the
C${}_{\text{2}}$ features in the blue end, can be attributed to two factors.
Firstly, the flux values in the blue end are too small relative to the entire
spectrum, making the C${}_{\text{2}}$ features much more insignificant
compared to the prominent CN features. Secondly, C2023 filters golden sample
carbon stars in the cold (most of them $G_{\mathrm{BP}}$ \- $G_{\mathrm{RP}}$
¿ 2; Gaia Collaboration et al. 2023) CSTAR sample based on the prominent CN
features rather than the much weaker C${}_{\text{2}}$ features. This results
in more than half of the targets in the golden sample of carbon stars having
almost no C${}_{\text{2}}$ features. Therefore, the C${}_{\text{2}}$ features
in Gaia’s XP spectra are not distinguishing for the identification of cold
carbon stars.

###  4.3 451 new carbon star candidates

We combined the model training and validation accuracy, as well as the changes
in FP and FN in Section 4.1.1, and ultimately determined the training epoch
number of the model as 16. We applied the model trained based on the dataset
provided in Section 3 to the 74,740 stars of CSTAR_XPM_nonG, and 451 new
carbon star candidates were obtained. The spectra of these stars exhibit
relatively prominent CN molecular bands compared to typical giant stars but
are weaker than most golden sample stars. Their spectra show reduced flux at
the blue end and hardly any significant C2 features are visible. In Figure
below 15, we show some spectra of these sampl stars.

Figure 15: Eight spectra were selected from the candidates. The spectra are
normalized similarly to those in Figure 3. The yellow areas are the key
features we summarized. The blue dashed line marks the position of the CN
characteristic troughs and the red marks the position of the peaks.

The color-magnitude diagram and distribution position of the CN band head
strengths are shown in Figure 16.

Figure 16: The color-magnitude diagram of new candidates (left), The molecular
band head strength of new candidates (right).

The color-magnitude map shows that our new carbon star candidates occupy the
location of the carbon-star golden sample but mainly concentrate in the hotter
regions. From the strength map, our new carbon star candidates mainly occupy
the position of weak CN band carbon stars mentioned in C2023.

The spatial distribution of the candidates is shown in Figure 17. Since the
positive sample is from the golden sample of carbon stars, they are
distributed throughout the Galactic inner disk, the Large Magellanic Cloud,
and the Small Magellanic Cloud. However, the new candidates are distributed
only in the Galactic inner disk. Therefore, the reason for their low flux
values at the blue end of the spectrum is likely due to severe extinction of
blue light caused by interstellar dust. This also explains why many of them
became warmer ($(G_{\text{BP}}-G_{\text{RP}})_{0}$ ¡ 2) after 3D extinction
correction (Green et al. 2019). Similar to the study of C-rich stars in the
Milky Way’s bar-bulge (Sanders & Matsunaga 2023), these stars are highly
likely to be the result of binary evolution. However, it is also possible they
are the result of single star evolution.

Figure 17: Spatial distributions of our results are drawn in equatorial and
Galactic coordinates. The gray spots denote the 15,400 golden sample of carbon
stars. The blue spots represent the 8176 positive sample from CSTAR_XPM_G. The
red spots are the 451 new carbon star candidates from CSTAR_XPM_nonG.

After cross-matching our 451 new carbon star candidates with LAMOST DR10, four
common sources were obtained, two of which were labeled as carbon stars by
LAMOST’s pipeline. Cross-matching with SIMBAD yielded 73 common sources, 28 of
which were labeled as C* in the main_type or other_types attribute. We will
give the detailed results in Table LABEL:star_table.

##  5 Discussion and Analysis

According to Wallace (1962), we plotted all CN absorption line positions at
the red end of the spectrum, which will be detailed in the dotted line in the
following Figure 18 and Figure 20.

###  5.1 Cross-match with LAMOST DR10

We cross-matched the RA and Dec of all 83,028 CSTAR_XPM sample stars with the
catalog of low-resolution spectra from LAMOST DR10, then we obtain the
subclass types given by the LAMOST’s pipeline (Wei et al. 2014). There are 442
sources labeled as carbon stars, and 1131 are not. Among the sources not
labeled as carbon stars, there are 351 G-type stars (31%), 232 K-type stars
(20.5%), and 509 M-type stars (45%), the vast majority of them come from
CSTAR_XPM_nonG.

To better understand the formation ofcharacteristic troughs in the $Gaia$ XP
spectra of carbon stars, the LAMOST spectrum of the same source was drawn in
Figure 18 to determine the absorption position of CN. Using the fast Fourier
transform and interpolation method (Zhang et al. 2021), the LAMOST spectrum
was convolved to the resolution of the approximate $Gaia$ XP spectra to show
the difference in the position and range of the absorption region caused by
different resolutions. We also plotted the same figures for common sources of
G, K, and M type stars as a comparison to illustrate the areas of feature
differentiation.

Figure 18: The spectra with wavelength range intercept to 4400-9000 Å. The
LAMOST spectra of the common sources are plotted with black lines. The red and
yellow lines are the result of their convolution to 50 resolution. The blue
and green lines are the $Gaia$ spectra. In the bottom panel, the light grey
shaded areas are the peaks and the red dashed lines in the bottom panels mark
several obvious CN absorption locations.

We found that the most significant feature shared by the negative sample of G,
K, and M type stars compared with the positive sample of carbon stars is the
Ca II line absorption feature in the range of 8484-8662 Å (Li et al. 2023),
which can also be identified in the $Gaia$ XP spectra. When the temperature is
low near the M type, the prominent TiO5 (7126-7135 Å; Li et al. 2023) and K I
(7664.9-7699.0 Å; Lépine et al. 2003) absorption line features will be
highlighted, which are also visible in the LAMOST and $Gaia$ XP spectra after
resolution reduction. During the change of spectral type from M to G, the
absorption characteristics of TiO5 and K I gradually weaken in the spectra
until they almost disappear. This is consistent with the Figure 19 interpreted
by our model. However, carbon stars exhibit broader and deeper CN absorption
features at 6900-7300 Å, and 7900-8300 Å, and a set of three peak tops at 6830
Å, 7780 Å, and 9000 Å (Sanders & Matsunaga 2023). While M-type non-carbon
stars also display peak structures, such as the peak top at 7050 Å, and 7600
Å, but their positions are entirely different. Therefore, the characteristic
peaks and troughs in the spectra of carbon stars are crucial distinguishing
features from other types of stars.

Figure 19: The wavelength range is intercepted from 6500 to 10000 Å. The
average distribution of SHAP values per features for these three types was
plotted respectively.

From Figure 19, we can also see that the absorption characteristics near 8300
Å and 9400 Å are gradually enhanced and the sum of SHAP values increases with
the temperature increase. Therefore, we hypothesize that in the
CSTAR_XPM_nonG, G, and K type stars are more likely to be carbon star
candidates. After a visual inspection of the 1131 stars, we estimate that G
stars may contain 10% candidates and K type stars may contain 2.6% candidates.
The reason is that if a low-temperature star has a high abundance of C/O, its
red-end molecular band characteristics will be more obvious. It will be easier
to select by the method of C2023. On the contrary, if at high-temperature
conditions, the carbon molecular band characteristics would be relatively
weak. In other words, if a source can still show the carbon characteristics at
high temperatures, its C abundance can be considerable.

###  5.2 Comparison with the results of Sanders & Matsunaga (2023)

Lebzelter et al. (2023) released a candidate catalog of 1,720,588 long-period
variable stars (LPVs). After cross-matching with $Gaia$ DR3 data, we found
that there are 13,513 stars in the golden sample of carbon stars and only 1289
stars in CSTAR_XPM_nonG are classified as LPVs. Based on the UMAP unsupervised
algorithm, Sanders & Matsunaga (2023) have presented 23,737 C-rich candidates
from the LPV candidate catalog above. There are 40 of the above 1289 LPVs in
the list of our 451 new carbon star candidates, and 26 are also identified as
C-rich stars by Sanders & Matsunaga (2023). There are 20 common objects
between our above 40 carbon candidates and 26 Carbon stars of Sanders &
Matsunaga (2023), where six objects are not included in our results. We show
their heat maps (see Figure 23).

We can see that Gaia 963568667249950592 (the first subplot of Figure 23) is an
obvious problematic spectrum, our model also assigns it a small confidence
level. The spectra in subplot 2-6 of Figure 23 have relatively weak absorption
features in the region around 8000 Å and 9300 Å, i.e., Feature 1 and Feature 2
in Table 4. Among these objects, Gaia 3586939158411536768 and
4056303080211083392 lack distinct CN molecular absorption band, then our model
gives them lower confidence. Gaia 4314316918298388992 is weak in all Features
1-5, so the confidence level is also low. Gaia 4339500235641264256 has a
stronger Feature 2, but a weaker Feature 1, resulting in a confidence level
not exceeding 0.5. The molecular absorption bands of Gaia 5546572717982752640
are not prominent enough, so the confidence level is low.

The majority of the golden sample of carbon stars are C-rich TP-AGB stars,
typical of long-period variable stars, so it is not surprising that they are
overwhelmingly classified as LPVs. Whereas the proportion of CSTAR_XPM_nonG
labeled as LPV candidates is small. This may be attributed, on the one hand,
they have a low signal-to-noise ratio or a lack of data points in the G time
series, small variability amplitudes, or short periods. Therefore they are not
detected as period signals by the SOS module of (Lebzelter et al. 2023). On
the other hand, they have not yet evolved to the intense third dredge-up phase
of TP-AGB like LPVs, which is consistent with their location site on the HR
diagram (Figure 9), or the dredge-up intensity of TP-AGB phase is
insufficient.Regarding why only a few dozen of the 1289 LPVs are classified as
C-rich stars, this could be due to selection effects. In CSTAR, the majority
of C-rich LPVs were selected by the algorithm of C2023 and included in the
golden sample, resulting in the remaining 1289 stars being predominantly
O-rich LPVs. Therefore, mining possible carbon stars from CSTAR_XPM_nonG is
not a simple task. Our candidates contain the vast majority of C-rich stars of
Sanders & Matsunaga (2023), which is a side effect of the high completeness of
our method. The other six stars are significantly different from the golden
sample, but there is no denying that they still are possible carbon stars. Our
model scores them with a not-low confidence level, just not above the 0.5
threshold. Although the positive sample our model learned is from the golden
sample of carbon stars with strong molecular bands, it is still able to
identify quite a few carbon stars exhibiting relatively weak CN molecular band
features. This indicates the good generalization ability of our model and
provides highly encouraging and informative results for our future systematic
exploration of carbon stars in all of $Gaia$’s XP spectra.

###  5.3 Explanation of key features

The positiona of the CN absorption lines with the wavelength greater than 6900
Å in the literature (Wallace 1962) were marked in the Figure 20.

Figure 20: In the common source, the model explains the distribution of the
most critical features (light gray areas) in the Gaia spectrum. The red
gradient against a light grey background and the deep pink energy bar at the
upper panel indicate the intensity of the features. The corresponding CN
molecular absorption positions (light blue dotted lines) are plotted in the
LAMOST spectrum.

From the figure, we can see that there are some differences in the position
and shape of the molecular absorption band regions in the spectrum at
different resolutions. From a positional perspective, the absorption regions
of the Gaia carbon star spectrum are relatively more backward than that of
LAMOST. From the perspective of shape, the continuous and dense CN absorption
line region in the LAMOST spectrum shows the form of wide and deep CN
characteristic trough in the Gaia spectrum, such as Feature 5 and Feature 2.
We posit that this is the combined absorption effect of superimposed multiple
molecular absorption lines in the spectrum at $Gaia$ spectral resolution, we
call it the ”mixed-molecular-absorption-band” in the Gaia spectrum. It can be
seen that Feature 1 and Feature 2 are consistent with the dense CN absorption
line regions except for the small SHAP value assigned by the absorption region
of Feature 5. which may be caused by the weaker signal flux values in the
spectra compared to Feature 1 and Feature 2.

As for the peak observed in a spectrum, itself is a reflection of the relative
decrease of opacity. When the molecular absorption on both sides of the peak
becomes stronger, it will increase opacity and thus form the characteristic
trough. Consequently, the region with less opacity stands out as a peak in the
spectrum. Hence, the characteristic peaks of Feature 3 and Feature 4 in the
shadow region are are also meaningful for distinguishing the positive and
negative samples.

###  5.4 Reliability analysis of new candidates

To verify the reliability of the new candidates, we mainly analyzed the
relationship between the molecular band head strengths (R773.3 and R895.0) and
the confidence levels, as well as the relationship between the strengths and
the SHAP values. We also calculated the average sum of SHAP values for the new
candidates.

Figure 21: The upper panels show the relationship between the confidence given
by the classification model and R773.3, R895.0; The bottom panels show the
relationship between the sum of SHAP values interpreted by the SHAP model for
each spectrum and R773.3, R895.0.

From Figure 21, it is clear that there is an approximate linear relationship
between the strengths and the confidence, SHAP value given by the models.
Moreover, it is evident that smaller values of R773.3 and R895.0 correspond to
lower confidence and SHAP values from the color-coded scatter plots.
Conversely, as the values of R773.3 and R895.0 increase, so do their
confidence and SHAP scores (dark color). Such strong positive correlations
proved that our candidates are considered to be reliable.

The average SHAP values distribution of features is shown in Figure 22. As can
be seen from this figure, the sum of average SHAP values (0.76164) of the 451
carbon star candidates is slightly lower than that of the positive sample
(0.99894; see Figure 14), but much larger than that of the negative sample
(0.00056), which also confirm the reliability of our candidates. We also see
that the contribution of Feature 5 of the candidates almost disappears and
Feature 2 becomes weaker compared to the positive sample, which is consistent
with the relatively weaker carbon signatures of their spectra. While the
negative feature near 8500 Å (Ca II absorption lines) becomes stronger, this
seems to imply that the candidates have greater metallicity (see Fig.8; Mallik
1997). It makes sense that our model treats this as a negative contribution to
carbon star identification, this is because metallicity tends to be positively
correlated with the oxygen levels of a star’s atmosphere (Sanders & Matsunaga
2023). Another possible reason is that the neighboring molecular bands become
weaker due to temperature or C/O ratio effect, thereby highlighting the Ca II
absorption line features. Although there is a slight difference in the
statistical distribution of SHAP values between the new candidates and
positive sample, the overall characteristics suggest that the model is more
inclined to classify them as carbon stars, indicating that these candidates
are more similar to the golden sample of carbon stars.

Figure 22: The same wavelength ranges from 6500 to 10000 Å. The average SHAP
values across feature distributions are calculated for 451 new carbon
candidates.

###  5.5 Origins of the new carbon star candidates

Investigating the origins of the 451 new carbon star candidates will be very
intriguing, we try to discuss and explain. First, according to the dividing
line between carbon dwarfs and carbon giants proposed by Li et al. (2024),
(i.e., $(M_{G})_{0}$ = 5.0 mag), their luminosities are all less than 3, thus
not meeting the criteria for dwarf carbon stars. Additionally, only 40 of
these candidates have been reported as LPV candidates by Lebzelter et al.
(2023). This implies that many of them might not have evolved to the intense
third dredge-up phase of the TP-AGB, and their extinction-corrected (Green et
al. 2019) color indicates that many of them are warmer than CN carbon stars.
These reasons suggesting they are likely extrinsic carbon stars, products of
binary evolution. However, it cannot be excluded that some of them are in the
TP-AGB phase, but their third dredge-up event and the loss of mass are
relatively mild (so warmer). Consequently, their variability amplitude is
small, and they were not detected by the algorithm of Lebzelter et al. (2023).
The mild dredge-up intensity indirectly suggests their lower masses (Sanders &
Matsunaga 2023) or that they are in the early stages of the TP-AGB phase,
where the dredge-up process has not fully developed. These could explain why
their spectral carbon features (such as CN and C2) are weak, although the C/O
ratio has already begun to exceed 1. Of course, the higher initial oxygen
abundance of the star is also a contributing factor, as can be glimpsed from
the speculation on metallicity in the previous section. Finally, an accretion
origin is also possible (Sanders & Matsunaga 2023), where the stars accumulate
carbon by accreting surrounding carbon-rich interstellar material, thereby
becoming carbon stars. Equatorial and Galactic coordinate maps indicate that
they are mostly distributed in a dense, band-like region within the inner
Galactic disk.

In summary, the 451 new candidates are likely external-origin carbon stars,
having evolved to the giant star phase. However, the possibility of single
star evolution cannot be excluded, where these carbon stars might have
undergone mild third dredge-up processes or inherited carbon through accretion
from surrounding interstellar material.

##  6 Conclusions

In this work, we have demonstrated the capability of the proposed ‘GaiaNet’
model in classifying carbon stars within the CSTAR_XPM sample. Additionally,
we have provided compelling feature attribution explanations for the model’s
outputs using the SHAP method. This ”magic” empowers the originally
incomprehensible black-box model with the ability to proactively ”speak up”
and make its discriminations traceable. Our main results and conclusions are
as follows:

  1. (i)

The method presented in C2023 for calculating the strength of molecular bands
has proven effective in selecting carbon stars with strong CN features and
temperatures below 3,800K. However, it lacks effectiveness in identifying
warmer carbon stars with weaker CN features and may introduce contamination.
In contrast, the proposed ‘GaiaNet’ model demonstrates superior capability in
distinguishing between carbon and non-carbon stars with good robustness from
the training datasets, particularly in identifying carbon stars with weak CN
features from CSTAR_XPM_nonG. Compared to four conventional machine learning
methods, the ‘GaiaNet’ model exhibits an average accuracy improvement of
approximately 0.3% on the validation set, with the highest accuracy reaching
up to 100%. The flexible parameter architecture provides the model with
enhanced fitting capability and stability achieved through the incorporation
of operations such as 1×1 Convolution, Global Average Pooling, Batch
Normalization (BN), and Dropout layers. Furthermore, Optimised training
methods including momentum and weight decay were employed. By utilizing deep
learning models with varying convolutional kernels during the 1D convolution
process, the ‘GaiaNet’ model effectively captures both local and global
features from low-resolution smoothed spectra. Coupled with the nonlinear
capabilities provided by an appropriate model’s depth, the model possesses
natural advantages in handling $Gaia$’s XP spectra.

  2. (ii)

We utilized SHAP for feature attribution interpretation on the ‘GaiaNet’
classification model based on positive sample. We provided clear feature
heatmaps for each spectrum, which serve as a basis for judgment and
identification of each spectrum. Based on the statistical analysis of SHAP
values from the spectra, we found that the C${}_{\text{2}}$ features are
almost non-distinguishing for the identification of cold carbon stars due to
their weak signal at the blue end of low flux values and the minimal presence
in many cold carbon stars. We summarised five key spectral feature regions:
6780-7460 Å, 7600-7820 Å, 7900-8360 Å, 8820-9060 Å, and 9140-9560 Å. Among
them, the 7900-8360 Å and 9140-9560 Å regions are the most obvious trough
areas. These two regions contribute over 80% due to dense CN molecular
absorption. The 7600-7820 Å and 8820-9060 Å regions contribute about 16%,
which are are also important distinguishing features due to the peak areas
from the strong CN molecular absorption on both sides. The 6780-7460 Å region
contains relatively weak peaks and troughs, making the smallest contribution
less than 4$\%$. We hypothesize that this is related to the spectral values,
which seems to be addressable through improved spectral normalization methods.
These results suggest that our model can effectively learn and capture the
critical distinguishing features that align well with CN molecular absorption
lines and visual identifications by astronomers when identifying carbon star
spectra. These five key feature regions will be essential for future
systematic mining of carbon stars from $Gaia$’s XP spectra.

  3. (iii)

There are 451 new carbon star candidates finally found from CSTAR_XPM_nonG
using the trained classification model. Compared with the method of simply
calculating band head strength, our approach effectively identifies more
carbon stars with weaker CN features. Compared to the golden sample carbon
stars, these new candidates exhibit relatively weaker CN molecular band
features and occupy the positions of weak CN feature carbon stars illustrated
in Fig. 21 presented in C2023. This means that the carbon abundances in their
atmospheres may be relatively low, with C/O ratios not significantly greater
than 1. These candidates are mainly concentrated in the Galactic inner disk
regions near the center. Additionally, we also noted that most of them have
higher temperatures (smaller $(G_{\text{BP}}-G_{\text{RP}})_{0}$) and few of
them are labelled as LPVs by Lebzelter et al. (2023), indicating the
possibility that they could be mostly extrinsic and less evolved carbon stars.
This is similar to the findings in Sanders & Matsunaga (2023). Moreover, our
candidates exhibit stronger negative contributions caused by Ca II absorption
features near 8500 Å compared to the positive sample from the statistical
analysis of SHAP values. Therefore, some of these stars are also likely to be
AGB stars with relatively mild third dredge-up processes and higher
metallicity, while also not excluding the possibility that they were
contaminated by a carbon rich environment near the Galactic center (Abia et
al. 2002).

  4. (iv)

We cross-matched the 451 new candidates with LAMOST DR10, resulting in a total
of 4 common sources. Among them, 2 were identified as carbon stars by the
LAMOST 1D pipeline. By cross-matching with the SIMBAD database, we obtained 73
common sources. Among these, 28 had already been reported as carbon stars or
carbon star candidates in previous work. A reference sample of 1,289 LPVs was
constructed to compare the results in Sanders & Matsunaga (2023). In this
work, we have discovered 40 potentially carbon stars from it, and they
identified 26 C-rich stars including only 6 stars that are not in our
candidates. This is due to the significantly different spectral features of
the six stars from the positive sample, which results in their confidence
lower than 0.5. These results indicate that the reliability and completeness
of our candidates are high. Furthermore, the assigned confidence and SHAP
values for each spectrum were used to quantify the model’s output, revealing a
strong correlation with the band head strengths of R773.3 and R895.0. This
suggests that the classification model primarily focuses on the crucial CN
features, and the learned features make sense. Hence, our model can
efficiently learn the key features between carbon and non-carbon stars in the
CSTAR_XPM sample, and the model’s reliability is further validated.

Compared to the positive sample, Most CSTAR_XPM_nonG exhibit distinctly
different distributions (much weaker) on the band head strength diagrams.
However, our trained model still identifies 451 credible potential carbon
stars within it, which also demonstrates the generalization ability of the
model. These carbon stars exhibit higher temperatures and weaker CN features,
making them undetectable by the C2023 method. Due to the limited diversity of
the training positive sample, there are still some carbon stars with weaker CN
features that are not recognized by our algorithm. Their probabilities given
by the model’s output do not exceed the threshold of 0.5 but significantly
larger than 0, thus the proportion of candidates in CSTAR_XPM_nonG could be
higher, especially the G and K type stars, predictably if we set a lower
threshold. we roughly estimate that the loss percentage is at least about 5%.
In the next work, we will make reasonable improvements to the classification
algorithm and incorporate more carbon stars with weaker CN bands into positive
sample to search for and identify more potential warm candidates from CSTAR.
We will then expand the algorithm to the entire $Gaia$ XP spectral dataset,
thereby significantly enriching the number of known carbon stars. Indeed it
would be interesting to perform further classification and meaningful
analysis. Given the resolution of $Gaia$’s XP spectra, our method can be
directly transferred to the upcoming China Space Station Telescope (CSST),
especially for the interpretation of different types of stars whose spectral
properties are still unclear.

## Acknowledgments

We thank Bo Zhang for his effective spectral smoothing method and helpful
discussions. We thank Xiao Kong for his help with the data. This study is
supported by the National Natural Science Foundation of China (NSFC) under
grant Nos. 12173013, 12273078, and U1931209; the National Key Basic R$\&$D
Program of China (No. 2019YFA0405500); the project of Hebei provincial
department of science and technology under the grant number 226Z7604G, and the
Hebei NSF (No. A2021205006). W.Y. thanks the science research grants from the
China Manned Space Project, and we also acknowledge the grants from the China
Manned Space Project under Nos. CMS-CSST-2021-A10. 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). This research has used
the SIMBAD database, operated at CDS, Strasbourg, France. This research also
used the low-resolution spectral data of LAMOST DR10 (the Large Sky Area
Multi-Object Fiber Spectroscopic Telescope DATA RELEASE 10;
https://www.lamost.org/dr10/v1.0). LAMOST is operated and managed by the
National Astronomical Observatories, Chinese Academy of Sciences.

## References

  * Abia et al. (2020) Abia, C., de Laverny, P., Cristallo, S., Kordopatis, G., & Straniero, O. 2020, A&A, 633, A135
  * Abia et al. (2022) Abia, C., de Laverny, P., Romero-Gómez, M., & Figueras, F. 2022, A&A, 664, A45
  * Abia et al. (2002) Abia, C., Domínguez, I., Gallino, R., et al. 2002, ApJ, 579, 817
  * Alksnis et al. (2001) Alksnis, A., Balklavs, A., Dzervitis, U., et al. 2001, Baltic Astronomy, 10, 1
  * Breiman (2001) Breiman, L. 2001, Machine learning, 45, 5
  * Carrasco et al. (2021) Carrasco, J. M., Weiler, M., Jordi, C., et al. 2021, A&A, 652, A86
  * Chen & Guestrin (2016) Chen, T. & Guestrin, C. 2016, 785
  * Chen (2015) Chen, Y. 2015, Master’s thesis, University of Waterloo
  * Cortes & Vapnik (1995) Cortes, C. & Vapnik, V. 1995, Machine learning, 20, 273
  * Cover & Hart (1967) Cover, T. & Hart, P. 1967, IEEE transactions on information theory, 13, 21
  * Efron (1979) Efron, B. 1979, The Annals of Statistics, 7, 1
  * Gaia Collaboration et al. (2023) Gaia Collaboration, Creevey, O. L., Sarro, L. M., et al. 2023, A&A, 674, A39
  * Green et al. (2019) Green, G. M., Schlafly, E., Zucker, C., Speagle, J. S., & Finkbeiner, D. 2019, ApJ, 887, 93
  * He et al. (2022) He, X.-J., Luo, A. L., & Chen, Y.-Q. 2022, MNRAS, 512, 1710
  * Hinton et al. (2012) Hinton, G. E., Srivastava, N., Krizhevsky, A., Sutskever, I., & Salakhutdinov, R. R. 2012, arXiv preprint arXiv:1207.0580
  * Ioffe & Szegedy (2015) Ioffe, S. & Szegedy, C. 2015, in International conference on machine learning, pmlr, 448–456
  * Ji et al. (2016) Ji, W., Cui, W., Liu, C., et al. 2016, ApJS, 226, 1
  * Kim (2014) Kim, Y. 2014, in Proceedings of the 2014 Conference on Empirical Methods in Natural Language Processing (EMNLP) (Doha, Qatar: Association for Computational Linguistics), 1746–1751
  * Kontizas et al. (2001) Kontizas, E., Dapergolas, A., Morgan, D. H., & Kontizas, M. 2001, A&A, 369, 932
  * Krogh & Hertz (1991) Krogh, A. & Hertz, J. 1991, Advances in neural information processing systems, 4
  * Lebzelter et al. (2023) Lebzelter, T., Mowlavi, N., Lecoeur-Taibi, I., et al. 2023, A&A, 674, A15
  * LeCun et al. (1998) LeCun, Y., Bottou, L., Bengio, Y., & Haffner, P. 1998, Proceedings of the IEEE, 86, 2278
  * Lépine et al. (2003) Lépine, S., Rich, R. M., & Shara, M. M. 2003, The Astronomical Journal, 125, 1598
  * Li et al. (2023) Li, J., Long, L., Zhong, J., et al. 2023, ApJS, 266, 4
  * Li et al. (2024) Li, L., Zhang, K., Cui, W., et al. 2024, ApJS, 271, 12
  * Li et al. (2018) Li, Y.-B., Luo, A. L., Du, C.-D., et al. 2018, ApJS, 234, 31
  * Lin et al. (2013) Lin, M., Chen, Q., & Yan, S. 2013, arXiv preprint arXiv:1312.4400
  * Lippmann (1987) Lippmann, R. 1987, IEEE Assp magazine, 4, 4
  * Lundberg & Lee (2017) Lundberg, S. M. & Lee, S.-I. 2017, in Advances in Neural Information Processing Systems 30, ed. I. Guyon, U. V. Luxburg, S. Bengio, H. Wallach, R. Fergus, S. Vishwanathan, & R. Garnett (Curran Associates, Inc.), 4765–4774
  * Mallik (1997) Mallik, S. V. 1997, A&AS, 124, 359
  * Morgan & Hatzidimitriou (1995) Morgan, D. H. & Hatzidimitriou, D. 1995, A&AS, 113, 539
  * Palmer & Wing (1982) Palmer, L. G. & Wing, R. F. 1982, AJ, 87, 1739
  * Polyak (1964) Polyak, B. T. 1964, Ussr computational mathematics and mathematical physics, 4, 1
  * Qin et al. (2019) Qin, L., Luo, A. L., Hou, W., et al. 2019, ApJS, 242, 13
  * Sanders & Matsunaga (2023) Sanders, J. L. & Matsunaga, N. 2023, MNRAS, 521, 2745
  * Secchi (1869) Secchi, A. 1869, Astronomische Nachrichten, 73, 129
  * Shang et al. (2022) Shang, L.-H., Luo, A. L., Wang, L., et al. 2022, ApJS, 259, 63
  * Si et al. (2014) Si, J., Luo, A., Li, Y., et al. 2014, Science China Physics, Mechanics, and Astronomy, 57, 176
  * Si et al. (2015) Si, J.-M., Li, Y.-B., Luo, A. L., et al. 2015, Research in Astronomy and Astrophysics, 15, 1671
  * Sutskever et al. (2013) Sutskever, I., Martens, J., Dahl, G., & Hinton, G. 2013, in International conference on machine learning, PMLR, 1139–1147
  * Szegedy et al. (2015) Szegedy, C., Liu, W., Jia, Y., et al. 2015, in Proceedings of the IEEE conference on computer vision and pattern recognition, 1–9
  * Wallace (1962) Wallace, L. 1962, ApJS, 7, 165
  * Wei et al. (2014) Wei, P., Luo, A., Li, Y., et al. 2014, AJ, 147, 101
  * Werbos (1988) Werbos, P. J. 1988, Neural networks, 1, 339
  * Zhang et al. (2021) Zhang, B., Li, J., Yang, F., et al. 2021, ApJS, 256, 14

##  Appendix A The Heat Map of Six diffent stars

Figure 23: From top to bottom, left to right are subfigures 1-6. The green
dashed lines mark the positions used to calculate the CN band head strengths.
7732.9 Å and 8950.0 Å are the top of the band head, others are used to
identify the location of band troughs.

##  Appendix B New Carbon star candidates from CSTAR_XPM_nonG

The catalog has 451 entries, which are all the most likely carbon star
candidates selected by our model. The model obtains this result by training
positive sample selected from CSTAR_XPM_G and negative sample selected from
CSTAR_XPM_non_G. All candidates are sorted in descending order of confidence,
and we show the top 100 candidates.

Table 5: The Catalog of carbon star candidates in CSTAR_XPM_nonG Source ID (GAIA) |  R.A (deg) |  Dec. (deg) |  Main type (SIMBAD) |  Other types (SIMBAD) | Confidence | SHAP value
---|---|---|---|---|---|---
5938465054747237248($S$) | 255.244434 | -48.950652 | LPV* | *—C*?—LP*—LP?—MIR—NIR | 1.0 | 0.99955
5835730639669683968($S$) | 241.084489 | -57.364918 | LPV* | *—C*?—LP*—LP?—MIR—NIR | 0.99999 | 0.99957
5313143689277626880($S$) | 140.040299 | -52.485428 | C* | *—C*—C*?—LP?—MIR—NIR | 0.99994 | 0.99955
5489622585506231552 | 116.27471 | -53.090801 | Star | *—IR—NIR | 0.99991 | 0.9703
656271591242470528 | 125.428401 | 17.285119 | S* | C*—IR—LP*—MIR—Mi*—NIR | 0.99972 | 0.99313
5941474349332334976 | 245.800041 | -48.712454 |  |  | 0.99755 | 0.9972
5854048743845734528 | 212.221097 | -62.928023 | Star | *—IR | 0.99754 | 0.99715
5854243902825896064 | 215.206128 | -62.274869 |  |  | 0.99737 | 0.99694
5980488835879061888 | 259.325419 | -30.927344 |  |  | 0.99676 | 0.99641
5604270350862378752($S$) | 106.535192 | -31.87378 | LPV* | *—C*?—IR—LP*—LP?—NIR | 0.99597 | 0.99524
4061196422397873920 | 262.590675 | -28.135145 |  |  | 0.99586 | 0.99566
4062631113240876672 | 268.808798 | -28.434779 |  |  | 0.99498 | 0.99456
3047772534862552320($S$) | 110.696578 | -8.815217 | Mira | *—C*—C*?—IR—LP*—LP?—Mi* | 0.99465 | 0.99262
4094930504172297984 | 271.316536 | -20.239457 |  |  | 0.99451 | 0.99423
5865171467859742848 | 200.505129 | -63.584571 |  |  | 0.99428 | 0.99374
4059056016839338752 | 261.375858 | -30.012123 | LP*_Candidate | *—LP?—NIR | 0.99413 | 0.99307
514263548487890176 | 36.376851 | 63.054981 |  |  | 0.9938 | 0.99351
5862198598937747456 | 199.179033 | -63.386326 |  |  | 0.99335 | 0.993
354762791118252288($S$) | 33.533472 | 47.662447 | C* | *—C*—C*?—IR—LP*—MIR—Mi* | 0.99306 | 0.99121
5716413183908004224($S$) | 114.258845 | -19.548336 | C* | *—C*—C*?—IR—LP*—LP?—Mi* | 0.99174 | 0.99082
4070093150664753792 | 271.441386 | -21.348017 |  |  | 0.9917 | 0.99131
2014593550230928896 | 342.495667 | 60.299062 | RedSG* | *—IR—LP*—MIR—Mas—NIR—pA? | 0.99072 | 0.99034
4045821465839438976 | 274.444504 | -32.26074 | Mi*_Candidate | *—IR—MIR—Mi?—NIR—V* | 0.99022 | 0.98829
5990083449153819392 | 243.082181 | -46.963369 |  |  | 0.98988 | 0.98948
5937319015683955968 | 250.688189 | -50.815719 |  |  | 0.98822 | 0.98793
5255092670818771968 | 153.907807 | -59.55129 | deltaCep | *—MIR—NIR—V*—cC* | 0.98572 | 0.98533
2178053614615673856 | 324.821403 | 55.932396 |  |  | 0.98554 | 0.9851
4058721456082780544 | 262.301994 | -29.531504 |  |  | 0.98474 | 0.98036
2029261378964662272 | 300.317391 | 28.786557 | Star | *—IR—MIR—NIR | 0.98369 | 0.98328
6057834256157778432 | 182.445267 | -61.527859 |  |  | 0.98285 | 0.98253
5258718997589107712 | 154.820297 | -57.32384 | post-AGB* | *—AB?—IR—LP*—LP?—MIR—NIR | 0.98075 | 0.98049
4056545208914430720 | 268.940437 | -29.089513 |  |  | 0.98008 | 0.97941
4065758154277969408 | 271.411365 | -24.856418 |  |  | 0.97982 | 0.97961
4054558494476007296 | 264.153127 | -32.681742 |  |  | 0.97892 | 0.97856
2026948006494690048 | 298.042068 | 25.93909 |  |  | 0.97884 | 0.97853
6054537473621712128 | 185.729092 | -63.137329 |  |  | 0.97847 | 0.9781
4056219684743385600 | 269.089299 | -30.304209 |  |  | 0.97829 | 0.97679
4104914825724879360 | 278.259201 | -12.901917 |  |  | 0.97615 | 0.9758
5932639326833555328 | 242.63607 | -53.942337 |  |  | 0.97602 | 0.97562
2026983981142199296 | 297.299261 | 25.949449 |  |  | 0.97599 | 0.97557
4259206577970407936 | 282.889751 | -1.597829 |  |  | 0.97575 | 0.97536
5960473527976734208 | 261.260894 | -39.265847 |  |  | 0.97508 | 0.97463
4044123652438371456 | 269.878542 | -30.884273 |  |  | 0.97454 | 0.97414
5862947950437062656 | 191.303003 | -63.410645 |  |  | 0.97201 | 0.9716
4062586784782656896 | 269.371786 | -28.626903 |  |  | 0.97174 | 0.97085
4057240787453141376 | 267.31741 | -29.262486 |  |  | 0.96676 | 0.9664
4067075609721316992 | 266.855157 | -25.982856 |  |  | 0.96669 | 0.9663
5972266168138911488 | 258.397097 | -39.766671 | Star | * | 0.96648 | 0.96618
4066701402179072640 | 272.329522 | -22.653578 |  |  | 0.96605 | 0.96562
4062660525189492608 | 269.324272 | -28.23621 |  |  | 0.96576 | 0.96525
3451032236956436096 | 89.311802 | 32.377569 | S* | *—C*—C*?—IR—LP*—LP?—MIR | 0.96564 | 0.96266
448608341933318656 | 52.575367 | 56.563998 |  |  | 0.96447 | 0.96409
4068826998652116480 | 266.535475 | -23.246602 |  |  | 0.96299 | 0.96256
486835887330942464 | 54.455116 | 62.415504 |  |  | 0.96165 | 0.95998
5942937868055138304 | 249.816206 | -46.113176 |  |  | 0.96146 | 0.96111
4060674841494267904 | 266.021275 | -27.667819 |  |  | 0.96144 | 0.961
5886865420714514944 | 230.322616 | -54.987715 |  |  | 0.96101 | 0.96063
4062342693202139136 | 269.461274 | -29.236824 |  |  | 0.96024 | 0.95963
5311959480914027392 | 138.009931 | -51.782806 |  |  | 0.95941 | 0.95839
4064870814038165504 | 271.674298 | -25.950435 |  |  | 0.9592 | 0.95895
5972008611885974016 | 258.9437 | -40.182199 |  |  | 0.95847 | 0.95822
5338251170445346560 | 162.625494 | -60.525464 |  |  | 0.95468 | 0.95429
5972260022101577216 | 259.209972 | -39.244572 |  |  | 0.95202 | 0.95162
5618457383755161216($S$) | 112.738987 | -23.464088 | C* | *—C*—C*?—IR—LP*—LP?—MIR | 0.95178 | 0.95018
5942980542846977408 | 249.414567 | -45.784273 |  |  | 0.9516 | 0.95094
2030353670740758656 | 298.576264 | 29.883676 |  |  | 0.95156 | 0.95101
5968203992396813056 | 250.512108 | -41.840297 |  |  | 0.95154 | 0.95131
5974841263083171456 | 263.673464 | -35.429062 |  |  | 0.95146 | 0.95108
526430675440161536 | 15.541814 | 66.977695 |  |  | 0.95074 | 0.95003
4063963721314562304 | 269.949861 | -26.829915 |  |  | 0.9505 | 0.94625
4322164334588338816 | 293.372593 | 17.214312 | LPV* | *—LP*—MIR—NIR | 0.95031 | 0.94997
5878424087663723904 | 217.675214 | -60.971208 | Star | * | 0.95018 | 0.94967
4069866513830121856 | 270.629017 | -22.40119 |  |  | 0.94908 | 0.94872
4256920826432562944 | 279.872356 | -4.237485 |  |  | 0.94807 | 0.94733
4118583194001521920 | 267.444473 | -21.508462 | SB* | *—NIR—SB* | 0.94748 | 0.947
2201585087589197696 | 335.620434 | 60.982812 |  |  | 0.94573 | 0.94533
5334981845684641280 | 177.503771 | -61.224682 |  |  | 0.94472 | 0.94413
2031872508602623360 | 297.895438 | 29.923054 |  |  | 0.94363 | 0.94314
5884051972285306880 | 236.587553 | -56.50199 |  |  | 0.94315 | 0.94265
4319411299229958784 | 290.048676 | 13.672668 |  |  | 0.94273 | 0.94232
4056588678371817600 | 268.51291 | -29.083316 |  |  | 0.94228 | 0.94124
4255740779922180608 | 283.449916 | -3.247059 |  |  | 0.94114 | 0.94066
4117181282312558208 | 265.719645 | -22.047087 |  |  | 0.93976 | 0.93926
5940367789632743040 | 248.633299 | -50.116494 |  |  | 0.93939 | 0.9386
5972190026958636416 | 258.096814 | -39.644554 |  |  | 0.9393 | 0.93895
2137619448150914432($S$) | 297.93139 | 53.692079 | S* | *—C*—C*?—Em*—IR—LP*—MIR | 0.93918 | 0.93707
5878466526231182208 | 219.454802 | -61.054658 |  |  | 0.93825 | 0.93786
4104919842247578240 | 278.542311 | -12.778434 |  |  | 0.93789 | 0.93758
887462018563120128($L$)($S$) | 104.758731 | 28.814589 | C* | *—C*—C*?—IR—LP*—NIR | 0.93704 | 0.93209
4255321483686501504 | 282.459084 | -3.991967 |  |  | 0.93671 | 0.93612
5972269432314624256 | 258.456518 | -39.618216 |  |  | 0.93528 | 0.93473
4110745157943191936 | 259.92529 | -25.261413 |  |  | 0.93506 | 0.93471
5964551444805061632 | 253.839017 | -44.03463 |  |  | 0.93337 | 0.93295
5969733756664102912 | 253.454426 | -40.520653 |  |  | 0.93309 | 0.93164
4266241042905220352 | 283.816115 | 0.179313 |  |  | 0.9328 | 0.93237
4068727046189164416 | 266.400872 | -23.76739 |  |  | 0.93238 | 0.93201
5961670105863745792 | 263.017942 | -39.635199 |  |  | 0.93092 | 0.93027
4040149880018961920 | 266.545071 | -36.492316 |  |  | 0.9303 | 0.92944
5934333021361329792 | 247.918148 | -50.699238 |  |  | 0.92446 | 0.92401
5932728936967577216 | 240.743333 | -54.404695 |  |  | 0.92275 | 0.92243

  * •

Notes.
($L$) represents a carbon star labeled by LAMOST’s pipline.
($S$) represents a carbon star identified by Sanders & Matsunaga (2023).
Other types are given a maximum of 6 types due to table width limitations.