Unsupervised Machine Learning for the Classification of Astrophysical X-ray Sources

The automatic classification of X-ray detections is a necessary step in extracting astrophysical information from compiled catalogs of astrophysical sources. Classification is useful for the study of individual objects, statistics for population studies, as well as for anomaly detection, i.e., the identification of new unexplored phenomena, including transients and spectrally extreme sources. Despite the importance of this task, classification remains challenging in X-ray astronomy due to the lack of optical counterparts and representative training sets. We develop an alternative methodology that employs an unsupervised machine learning approach to provide probabilistic classes to Chandra Source Catalog sources with a limited number of labeled sources, and without ancillary information from optical and infrared catalogs. We provide a catalog of probabilistic classes for 8,756 sources, comprising a total of 14,507 detections, and demonstrate the success of the method at identifying emission from young stellar objects, as well as distinguishing between small-scale and large-scale compact accretors with a significant level of confidence. We investigate the consistency between the distribution of features among classified objects and well-established astrophysical hypotheses such as the unified AGN model. This provides interpretability to the probabilistic classifier. Code and tables are available publicly through GitHub. We provide a web playground for readers to explore our final classification at https://umlcaxs-playground


INTRODUCTION
X-ray astronomy is embarking on a promising era for discovery, facilitated by cutting-edge missions like eROSITA.This mission, which scans the entire X-ray sky, has already made its initial data release accessible to the scientific community (Predehl et al. 2021;Merloni et al. 2012).Prospective missions such as Athena and Lynx offer enhanced capabilities.Specifically, Athena combines the sensitivity of Chandra with the expansive coverage area of eROSITA (Nandra et al. 2013).Meanwhile, Lynx plans to surpass Chandra in terms of sensitivity and resolution (Gaskin et al. 2019).These missions will provide a wealth of new data and detections, catalyzing the identification of previously unrecognized sources.
The Chandra X-ray Observatory is one of the NASA's great observatories and its flagship mission for X-ray astronomy.Since it was launched to space in 1999, the telescope has observed several sources with two instruments, the Advanced CCD Imaging Spectrometer (ACIS) and the High Resolution Camera (HRC) (Wilkes & Tucker 2019).Accreting black holes in the center of galaxies, supernova remnants, X-ray binaries, and young, rapidly rotating magnetic stars are some of the most common targets that Chandra has identi-★ E-mail: samuelperez.di@gmail.comfied.As one of NASA's Great Observatories, Chandra has contributed to numerous groundbreaking discoveries in high-energy astrophysics over its 23-year lifespan.The Chandra Source Catalog (CSC) collects the X-ray sources detected by the Chandra X-ray Observatory through its history (Evans et al. 2010).It represents a fertile background for discovery, since many of this sources have not been studied in detail.A significant proportion of these sources, which range from young stellar objects (YSO) and binary systems (XB) to distant active galaxies housing supermassive black holes (AGN), remain largely unexplored.Furthermore, Chandra's data contains traces of exceptional phenomena like extrasolar planet transits, tidal disruption events, and compact object mergers.Despite the potential scientific wealth within Chandra's data, only a fraction of CSC sources have been classified.In order to conduct a thorough investigation of CSC sources, and to gear up for forthcoming large-scale X-ray surveys, we need to classify as many catalog sources as possible.
As larger astronomical surveys become available, researchers are increasingly adopting more sophisticated statistical models, including machine learning and data science methods, for classification tasks.A comparison of various supervised and unsupervised learning methods for the statistical identification of XMM-Newton sources was presented in Pineau et al. (2010).This study employed a probabilistic cross-correlation of the XMM-Newton Serendipitous Source Catalog (2XMMi) with catalogs such as Sloan Digital Sky Survey (SDSS) DR7 and Two Micron All-Sky Survey (2MASS).They compared algorithms like k-Nearest Neighbors, Mean Shifts, Kernel Density Classification, Learning Vector Quantization, and Support Vector Machines.
In Lo et al. (2014), an automatic classification method for variable X-ray sources in 2XMMi-DR2 using Random Forest was presented.The authors achieved approximately 97% accuracy for a 7class dataset.In Farrell et al. (2015), a catalog of variable sources in 3XMM was classified using Random Forest, training the classifier on manually classified variable stars from 2XMMi-DR2 and obtaining approximately 92% accuracy.
As interest in this research field continues to grow, unsupervised learning techniques have gained prominence in recent years.In Rostami Osanloo et al. (2019), an automated machine learning tool for classifying extra-galactic X-ray sources using multi-wavelength data, particularly from the Hubble Space Telescope, was proposed.In Ansari, Zoe et al. (2021), a probabilistic assignment was performed using mixture density networks (MDN) and infinite Gaussian mixture models to classify objects in the dataset as stars, galaxies, or quasars, achieving a 94% accurate split.The training data consisted of magnitudes from SDSSDR15 and the Wide-field Infrared Survey Explorer (WISE).
In the same line, in Logan, C. H. A. & Fotopoulou, S. (2020), an alternative unsupervised machine learning method for separating stars, galaxies, and QSOs using photometric data was presented.This approach employed Hierarchical Density-Based Spatial Clustering of Applications with Noise (HDBSCAN) to identify different classes in a multidimensional color space.Using a constructed dataset of approximately 50,000 spectroscopically labeled objects, the authors achieved F1 scores of >0.92 for star, galaxy, and QSO classes.
In recent years, significant efforts have been made to improve the classification of X-ray sources using supervised machine learning methods.Notably, Yang et al. (2022) developed a novel approach for classifying CSC version 2 (CSC2) sources.They successfully classified 66369 CSC2 sources (21% of the entire catalog) and performed focused case studies, demonstrating the versatility of their machine learning approach.Additionally, they provided insights into the biases, limitations, and bottlenecks encountered in these types of studies.Similarly, Kumaran et al. (2023) aimed to classify point Xray sources in the Chandra Source Catalog 2.0 into categories such as active galactic nuclei (AGN), X-ray emitting stars, young stellar objects (YSOs), high-mass X-ray binaries (HMXBs), low-mass X-ray binaries (LMXBs), ultra-luminous X-ray sources (ULXs), cataclysmic variables (CVs), and pulsars.Using multi-wavelength features from various catalogs, including CSC2, they used the Light Gradient Boosted Machine algorithm and achieved scores of ≥91% in precision, recall, and Matthew's Correlation coefficient.
Despite the recent efforts, supervised classification of X-ray sources remains a complex endeavor.Building a reliably labeled training set often involves manual review of the existing literature, which can sometimes contain ambiguous results.Moreover, a large number of these sources are missing optical or infrared data that could provide additional insights into their physical nature.Thus, unsupervised learning offers a distinct advantage as it eliminates the need for a pre-labeled training set.It can also reveal hidden patterns in the data that may not be immediately apparent.In this study, we introduce an unsupervised learning approach designed to classify as many CSC sources as possible, leveraging solely available X-ray data.By clustering the source observations by their similarities, and then associating these clusters with objects previously classified spectroscopically, we propose a new methodology to provide probabilistic classification for a numerous amount of sources.Furthermore, we offer not only a classification, but also a comprehensive discussion on the astrophysical implications arising from the observations made at each stage of the pipeline.We evaluate the strenghts and limitations of the method.
In Section 2, we describe the primary dataset used in this work, the Chandra Source Catalog, detailing the query and preprocessing steps required to make the data suitable for analysis with unsupervised learning methods.We also discuss the feature selection process and final properties used.In Section 3, we explain the methods employed, including Gaussian Mixtures and the development of a probabilistic classification algorithm based on the Mahalanobis Distance.We also introduce Soft/Hard voting classifiers that were used for the final output of this study.This section includes an overview of the techniques, a detailed explanation of the algorithms, relevant parameter choices, and a validation approach for the method.In Section 4, we provide a comprehensive exploration of results for each step of the algorithm.We report per-detection classifications, final master classes for 8,756 sources, and significant trends or patterns observed during the analysis.In Section 5, we discuss the implications of our results, including comparisons with previous work and catalogs.We analyze the astrophysical nature revealed by the data and evaluate the possible limitations of our study, addressing potential biases and sources of uncertainty.We also suggest directions for future research.Finally, in Section 6, we summarize the conclusions of the work, highlighting the main outcomes and their relevance.

DATA
In this section, we describe the data set used in the work at hand, and the properties chosen for the subsequent analysis.Our main data consists of Chandra Source Catalog version 2 per-detection registers.

The Chandra Source Catalog
The Chandra Source Catalog (CSC) collects and presents summarized properties for the X-ray sources detected by the Chandra X-ray Observatory through its history.Version 2.0 (CSC2), which we use here, is the second major release of the catalog, including properties for 317167 X-ray sources in the sky.Properties include source photometry (brightness), spectroscopy (energy), and variability (change over time).
The Chandra Source Catalog includes properties for 928280 source detections, which were detected in 10382 Chandra observations until 2014.The catalog presents information through 1700 columns of tabular data, distributed across five energy bands (broad, hard, medium, soft, and ultra-soft) for the Advanced CCD Imaging Spectrometer (ACIS), and a single band for the High Resolution Camera (HRC) in the wide spectrum.Chandra has been observing the universe in the 0.5-8 KeV band (Wilkes & Tucker 2019).The master source properties present summary measurements of stacked detections.Properties for detected sources are measured for both the stacked detection and the individual observations in all of the bands (Wilkes & Tucker 2019).
This study utilizes the individual properties measured for detections in the Chandra Source Catalog.Our initial classification is performed on a per-observation basis, before we evaluate the most probable class for each master source.This approach takes into account that CSC sources may have multiple detections.A subset of the CSC2 per-observation detection table was extracted using CSCView1 , with selection criteria restricting source detections to those with a broad band flux significance greater than 5 (flux_significance_b > 5), and non-null values for power-law gamma (powlaw_gamma) and blackbody temperature (bb_kt).These criteria aimed to ensure observations possessed sufficient significance to allow for fitting their spectra to a model, thereby enabling the extraction of statistical relationships from the data.This process yielded a table of 37878 unique observational entries.

Properties and preprocessing
In this subsection, we discuss the properties used in our analysis.We also outline the preprocessing steps applied to the dataset.These steps include data cleaning, transformation, and normalization.We conducted an iterative exploration process for feature selection.We selected 12 properties, primarily involving variability and spectral measurements.Most properties are extracted directly from the per-observation level of the CSC, except for var_ratio_* and var_newq_b.These were computed by us to summarize source information from variability mean, standard deviation, min and max.Table 1 lists the chosen properties and their descriptions.Our initial step was to cleanse the original data, selecting rows with valid entries (not NaN) for the chosen properties.This reduced the dataset to 29655 rows.Following this, we carried out either normalization exclusively or a combination of logarithmic transformation and normalization, depending on the distribution and range of the selected properties.Feature preprocessing was tailored to achieve scale uniformity, preventing large-scale features, like outliers, from dominating the learning process and facilitating faster convergence.Details of the chosen processing approach for each property are presented in Table 2.When required, the log transformation is performed prior to normalization, feeding its result into the normalization process.

Normalization
We employ the MinMaxScaler method from the scikit-learn Python library (Pedregosa et al. 2011) for data normalization.Using this, we scale each selected feature to the range [0, 1] based on its minimum and maximum values.The transformation can be represented as: (1) In equation 1, X denotes the original feature values, while X min and X max represent the minimum and maximum values of the feature, respectively.The resulting X  maintains the relative relationships between the original values.

Log transformation
We apply a natural logarithm to each data point using the numpy library (Harris et al. 2020).To handle zero values, we add one-tenth of the non-zero minimum value of each property to the data before taking the logarithm: (2) In the equations above,  * min represents the minimum non-zero value in the feature values X.The log-transformed property is represented by X log .

Feature selection
The efficacy of our proposed algorithm relies significantly on the selection of CSC properties used as features.Feature selection in unsupervised machine learning is an active field of research (Solorio-Fernández et al. 2020).For this study, we prioritized features that enhance information gain, i.e, those that maximize the separation between clusters.While dimensionality reduction techniques are often used in feature selection, we chose to follow a different criteria to prioritize interpretability.Additionally, feature scales vary significantly, adding complexity to the task.Selected CSC properties are described in Table 1, and they were chosen using the following criteria: (i) We start from astrophysical domain knowledge.Features that are more likely to inform the classification are those related with the spectral and time-domain properties, because they are associated to specific astrophysical processes (e.g.accretion).
(ii) Besides using astrophysically significant features, we also employ an empirical strategy.Here, we evaluate feature importance by the degree to which a given set of features yields distinct clusters.We do this by observing how the properties distribute in each cluster for various feature selections and combinations.
(iii) Then, we choose properties based on maximizing the amount of information in clusters while simultaneously minimizing redundancies.By minimizing redundancy, we can reduce the number of features required to accurately classify sources, which may lead to simpler and more interpretable models.At the same time, maximizing the amount of information in clusters ensures that the clusters formed are distinct and informative.This approach addresses the challenge of determining the optimal number of astrophysical x-ray properties required for accurate source classification, and provides insights into the most important properties for distinguishing between different source types.

METHODS
In this section we present the unsupervised learning and cluster-wise classification techniques employed in our analysis.We use Gaussian Mixtures Models (GMMs), which allow dataset clustering based on the multi-dimensional distribution of data points.In this study, selected CSC column properties serve as the features for our analysis.We first run a GMM on those features to find initial clusters.We then cross-match our CSC sources with the SIMBAD database in order to extract any existing labels assigned to the observations.We assign probabilistic classes to each unlabelled observation inside each group.The assignment is accomplished by examining the distance between the unlabelled observation and the centroids of labeled observations (extracted from SIMBAD) that share the same class within the cluster.Based on the detection classifications, we assign a final class to the master source.The next subsections offer a more detailed account of each step in this algorithm.

Gaussian Mixture Models and the EM algorithm
Clusters are first identified using a Gaussian Mixture Model (GMM) in the multi-dimensional space of the chosen CSC properties.In order to understand the assumptions and limitations of this method, we first have to describe the basics of its functionality as a clustering technique.A Gaussian Mixture represents a linear combination of  different Gaussian distributions (Bishop 2006) where   are the mixture coefficients or mixture weights, and The Gaussian distributions N (x | µ  ,   ) that compose the mixture are called components, each having its own covariance   and mean µ  .We use maximum likelihood estimation in order to adjust the parameters of the mixture of Gaussians to optimally represent clusters in the data.From 3, the log likelihood function is given by In the case of a single Gaussian distribution, an analytical solution for the maximum likelihood can be obtained easily.Indeed, in the simplest case (one dimensional) it corresponds to the mean and variance of the data.However, maximizing the likelihood of a Gaussian Mixture is not an easy task, lacking a closed-form analytical solution.A broadly applicable method solution, which we used in this paper, is the expectation-maximization algorithm, also abbreviated as EM algorithm.
The EM algorithm is a method for estimating the parameters that maximize the likelihood in unobserved latent variable dependent models.In this case, the latent variable is the cluster assignment of each data point.The EM algorithm was initially proposed in Dempster et al. (1977).A complete review of the method could be found in McLachlan & Krishnan (2007).
The EM algorithm consist of two main states, as the name suggests: the expectation (E) step and the maximization (M) step.The initial parameters of the Gaussian Mixture could be chosen in different approaches.The most common one would use random initial parameters.However, a more convenient routine is to assume identity covariances in the Gaussians and perform K-means clustering to find the means.A summary of the EM algorithm steps is described in online Appendix A.
We have provided a general description of the GMM and the EM algorithm.For the mathematical details, we refer the reader to Bishop (2006); Deisenroth et al. (2020).An analysis of how the log-likelihood increases in each of the steps of the EM routine is presented in Neal & Hinton (1998).
In this work, we use the GMM architecture provided in the scikitlearn library (Pedregosa et al. 2011) for Python.There is no direct evidence supporting the assumption that our data come from a mixture of gaussians.However, GMM is a flexible method capable of producing reliable results even when this assumption is not strictly met.We expect this minimizes the impact of the possible more complex distributions present in the data.The Gaussian Mixture serves as a foundational step in our classification process.Throughout our analysis, components and clusters will mean the same.The selection of hyperparameters for our study is explained in Section 3.3.

Per-detection classification
The initial clustering serves as a preliminary step towards classification.Given the complexity of astrophysical classes, they likely exceed the optimal number of clusters we will derive.This can be partly attributed to the potential overlap of various astrophysical classes within the feature space, particularly when relying solely on X-ray data.Consequently, a single cluster may encompass objects of multiple classes.However, it is expected that object detections within the same cluster will exhibit similar properties.
The SIMBAD Astronomical Database, operated at CDS, Strasbourg, France, (Wenger et al. 2000) provides a relatively solid knowledge base source for a systematic class extraction (Oberto et al. 2020a) 2 .SIMBAD employs a hierarchical system for classifying objects, which relies on catalogue identifiers.SIMBAD provides a base for associating probabilistic classses to unclassified object detections.Initially, we cross-match our dataset with classified objects in the SIMBAD database, leading to a minority of classified objects per cluster.We then compute the distance between each unclassified object and the classes centroid within the cluster.These distances facilitate probabilistic assignment, providing an indication of the detection's likelihood of belonging to a given class.
The crossmatch was performed using the CDS Upload X-Match functionality of TOPCAT (Taylor 2005), which is an interface for the crossmatch service provided by CDS, Strasbourg.The match radius was set to ≤1 ′′ and the Find mode was set to Each to obtain a new table with the best match from the SIMBAD database, if any, within the selected radius.In the absence of a match, the SIMBAD columns in the table would be blank.We specifically looked for the otype property, which refers to the main object type assigned to the sources (Oberto et al. 2020b).It is expected that a significant number of sources would not have a matching entry in the SIMBAD database or may have an ambiguous otype classification, such as peculiar emitters (Radio, IR, Red, Blue, UV, X, or gamma), assigned when no further information is available about the nature of the source.
We performed the classification of sources using the extracted labels and cluster information.To achieve this, we measured the similarity between the unlabeled and the labeled observations within a cluster using the Mahalanobis distance metric.
The Mahalanobis Distance was introduced in Mahalanobis ( 1936), as a distance metric of a point and a distribution.It is defined as: where µ and Σ are the mean and covariance matrix of a probability distribution  respectively, and x is a data point.
We chose to use the Mahalanobis distance as our metric for analysis because it takes into consideration the distributions of each class group within each cluster.This feature allows us to calculate the distance of a source detection without a label with respect to each class group in the cluster based on the shapes of their distributions.The computation is carried out using Equation 6.
After measuring the distance of a point  to all classes in a cluster space, we apply the softmin function to convert the distances to a probability distribution.The softmin is a variation of the well known softmax function (Goodfellow et al. 2016), and it is defined as: where   is a real number in a vector x = { 1 , . . .,   }. softmin(  ) is equivalent to softmax(−  ).After applying this function to a vector, every element will be in the range [0, 1] and   softmin(  ) = 1.We apply the softmin function for all unlabeled points in all clusters, resulting in smaller distances yielding larger probabilities.The final outcome is a table of source detections with a distribution of probabilities over classes.
It is worth highlighting that a single source detection may receive multiple probabilistic classifications.This process can be represented mathematically as follows.Suppose we have  detections unistra.fr/guide/otypes.htxfor a source, denoted by  1 ,  2 , . . .,   .For each detection   , let   be the set of possible classifications that it can receive.Then, we can define a collection of sets  1 ,  2 , . . .,   , where each set   corresponds to the possible classifications for detection   .Thus, the relation is summarized as where (  ) represents the probability that detection   belongs to each possible classification in the set   .This formalism represents classifications at the per-detection level, treating each detection as an individual target to classify.We refine this process to produce a single classification for each source.

Source master class
A fraction of the sources in the CSC have been observed more than once during the lifetime of Chandra.A source with multiple detections will appear more than once in the catalog.The classification scheme that we have described so far operates on individual detections, which implies that sources with multiple observations get multiple classifications.While you expect classifications of the same source to agree, this is not always the case, since detections can differ in their properties.We therefore require a procedure to provide master classification of sources, that takes into account the information from the multiple detection classifications.
To obtain the master class for a given source, we use two different approaches: soft voting and hard voting classifiers.In the hard voting approach, also known as majority voting, we simply choose the most common classification among the different detections of a source.In the soft voting approach, we profit from the probabilistic nature of the classification: for each detection of a source, we weight each class according to their probability, and then average over all detections in order to get the class with the highest mean probability.Both soft and hard voting have been widely used in ensemble methods of machine learning (Zhou 2012).
We then arrange the results in two tables: a table of uniquely classified sources, and a table of ambiguous sources that yield different classifications with the two different methods.
To formalize the process of obtaining the master class for a given source, we introduce the following mathematical notation.
Let  be a source, and let D  be the set of all detections associated with .The hard voting classifier assigns the class C  that appears most frequently among the detections in D  .That is, where  is the set of all possible classifications, and [ → ] is an indicator function that equals 1 if detection  was classified as , and 0 otherwise.We have that   is the set of all possible classifications for detection   ∈ D  .For each detection   , let   (  ) be the probability that   belongs to class  in   .We can then compute the soft voting approach.We calculate the mean probability for each class  in  across all detections, i.e., where The full classification procedure that we have described here will be applied to the pre-processed CSC dataset, and their properties.
Results of this process are presented in Section 4. A summary of our full pipeline is presented in Algorithm 1.

A. Gaussian Mixture Model and Class Extraction
1. Clustering step: Assign a cluster to all data points using a Gaussian Mixture Model.2. Crossmatch step: Crossmatch data with other databases (in this case, SIMBAD) to extract classes for all data points, if available.

B. Cluster-wise classification algorithm
1. Distance step: For a data point (detection)   , compute the Mahalanobis distance to all possible classes.Store them as where   represents all possible classes for detection   .

Probability step:
Convert the distances to a probability distribution over classes with softmin((  ))

3.
Repeat for all registers in the data set.

C. Master classification 1. Hard voting:
For each source  in the data set, consider each of its detections and assign a class using the hard voting approach, which selects the most common classification among the detections: 2. Soft voting: For each source  in the data set, consider each of its detections and assign a class using the soft voting approach, which takes into account the probabilities of each detection belonging to each class:

Tables:
Aggregate the sources that agreed on their classification in both the hard and soft voting methods into a table of uniquely classified sources.Aggregate the sources that were classified differently in each method in a table of ambiguous sources.
Algorithm 1: Our proposed classification method.This technique takes advantage of a previous clustering output and use it to bound the classification to particular cluster spaces.

Hyperparameter tuning
The most important hyper-parameter to determine is the number of components  (n_components in the GaussianMixture scikitlearn architecture).In order to select the number of components, we preliminarily use the Bayesian Information Criterion (BIC).
The BIC was introduced in Schwarz (1978) and is defined as: where L is the maximum likelihood of the model,  are the degrees of freedom or number of parameters estimated by the model, and  is the number of samples or data points.The BIC measures how much variance in the data is explained by the incorporation of additional components.Additional components can help increase the likelihood, but beyond a certain complexity no additional information is gained, as we are overfitting.The BIC balances the goodness of fit of the model with the number of parameters.As the number of components in a GMM increases, the fit to the data improves, but the number of parameters also increases, leading to a more complex model.The optimal number of components is the one that provides a good balance between the fit to the data and the model complexity.
Usually a lower BIC is preferred.However, BIC decreases as the number of components increases, suggesting that the model is fitting the data better, but it also means that the model is becoming more complex and may be overfitting.We need to choose the number of components with the lowest BIC score with a reasonable number of parameters.We explore this selection using the elbow method.
The elbow method is a well-known heuristic technique for determining the optimal number of components in a model.This method involves plotting the BIC (or originally, the dispersion) as a function of the number of components, and selecting the number of components where the decrease in BIC slows down, i.e., the number of components for which the slope of the tangent changes drastically.This indicates that adding additional components is not significantly improving the fit to the data.The elbow method was first introduced in Thorndike (1953).Recent studies have shown that the BIC-based component selection performs better than the traditional dispersionbased elbow method (Schubert 2022).
We plot the BIC as a function of the number of components to discern where further addition does not substantially improve data fitting.Figure 1 shows that the BIC plot lacks a clear elbow; therefore, we employ the gradient function to determine the optimal number of components.At  = 6, we note a considerable change in the gradient's slope, which subsequently converges to 0 for  ≥ 7.This observation indicates that augmenting the components beyond  = 6 does not notably enhance the data fit.
In our analysis, we have used the BIC-based elbow method as an initial step guide.We have also relied on domain knowledge and a graphical analysis of the data distribution to confirm the number of components.We performed this analysis by experimenting with several number of components over the same dataset.This multifaceted approach allows us to make informed decisions about the number of clusters and to ensure the reliability of our results.We provide an analysis on the stability of this selection, and its impact on the final classification in online Appendix C.
As a result of this process, we selected the number of components to be 6.Further graphical and distributional analysis revealed that the formation of the clusters was mainly influenced by a subset of the selected features (Table 1), specifically variability probabilities and hardness ratios.These findings are discussed in detail in Section 4.
Another crucial hyperparameter is the type of covariance parameters to use.These are defined by the covariance_type parameter in the GaussianMixture method of scikit-learn.It is challenging to determine the position and shape features of the different components a priori, i.e., we do not know how the covariance matrices of each Gaussian need to be in order to optimally fit the data points.Therefore, we chose the 'full' configuration, which allows each component to define its own covariance matrix (position and shape) indepen- dently (Pedregosa et al. 2011).To ensure reproducibility, we fixed the random_seed of our Gaussian Mixture Model to 42 3 .

Validation
Validating unsupervised learning techniques requires a departure from typical methods predominantly used in supervised learning research.The validation of our proposed algorithm not only confirmed its suitability for our classification problem but also guided us in assigning a 'master class' to unique sources.We do this by using a 70%/30% stratified random split over the set of SIMBAD labeled CSC detections.We treat the 30% split as a benchmark set.These detections are treated as if they were unlabeled during the classification process, but we then compare the model's predictions to the known SIMBAD labels to assess its performance.This operation is performed using the train_test_split function from the scikitlearn library.Its effectiveness is evaluated through the analysis of the resulting confusion matrices, which contrast the ground truth classes (SIMBAD) with our predicted classifications.These results are presented in § 4.2.2.

RESULTS
In this section we present the results of applying the proposed pipeline to the CSC dataset described in § 2. Using the GMM results, as well as the class associations using the Mahalanobis distance within each 3 "The answer to the ultimate question of life, the universe, and everything is 42." (Adams 1995) cluster, we assign probabilistic classes for each object detection in the dataset.We summarize classes for each individual object using hard and soft voting classification.We provide final uniquely classified and ambiguous classification tables for individual sources.We validate our results by constructing confusion matrices for a set of control detections, and by comparing our results to compiled catalogs of X-ray binaries, AGNs, etc.We also quantify the uncertainty in our probabilistic classifications by evaluating the shape of the probability distribution over classes, and by evaluating consistency in the assigned class for bona-fide CSC sources and regions.We provide the code, data and results in the GitHub repository https://github.com/samuelperezdi/umlcaxs.

Clustering
We apply the Gaussian Mixture Model with the hyperparameters described in Section 3.3 to the selected CSC features ( § 2) in order to identify an initial set of clusters.Members of a given cluster have relatively similar features among them, such as hardness ratios and variabilities.Features tend to be distinct from cluster to cluster, but this does not exclude significant overlap, or implies that all points within a cluster belong to the same class.We applied the clustering to a total of 29655 individual CSC source detections.
Figure 2 shows a Mollweide projection of the sky, where the spatial distribution of the clustered detections is presented.We note a clear differentiation between sources located on the galactic plane, versus off-the-plane sources, as they tend to respectively be assigned to different clusters.For instance, clusters 3 and 4 predominantly reside in the galactic plane, with 34.3% and 28.8% of detections at galactic latitudes || < 5 • , respectively.Conversely, clusters 2 and 1 predominantly consist of extragalactic source detections, with only 12.6% and 9.8% of detections at galactic latitudes || < 5 • , respectively.To a lesser extent, cluster 5 is also associated with offplane sources, while cluster 0 exhibits a stronger connection to the galactic plane.
In Figure 3 we show the hardness-to-hardness diagrams for objects in different clusters, color-coded by variability on the hard band.Overall, cluster membership relates to hardness and variability properties.Cluster 4 is made up almost entirely of highly variable objects irrespective of their spectral shape, whereas cluster 3 is predominately made of very low variability objects with a hard spectrum.Other clusters show more dispersion in variability, although with discernible differences in both the overall variability and the spectral shapes.Clusters 1 and 2 occupy an intermediate region in hardness.
No single feature by itself can separate sources of different properties, but combined they can isolate specific behaviors.

Classification
So far, our clusters contain a mix of SIMBAD-classified objects and a majority of unclassified objects.We first investigate if there are dominant classes within each cluster In Table 3 we rank the 10 most common classes within each cluster.Unclassified detections are marked as NaN, and are the majority of data points in all clusters, except for cluster 4. Some of the clusters are clearly dominated by object detections of similar astrophysical types.Such is the case of cluster 4, for which over 90% of the classified objects are of types associated with young stars.Clusters 0 and 5 also appear to be dominated by young stars, whereas for cluster 2 practically the entire set of SIMBAD-classified objects corresponds to either supermassive black holes (QSO, AGN, Seyfert), or stellar size black holes (HMXB, XB).The X-ray properties of accreting compact systems are expected to be similar, but we find clusters (e.g.cluster 1) where AGNs and QSOs are significantly more common than X-ray binaries.Despite these trends, clusters are not generally populated by objects of a given class.

Classes
Some of the classes represented in the SIMBAD database are extremely rare, having less than 10 examples each.Our initial BIC test, on the other hand, indicates that the X-ray properties alone contain only a limited amount of information to separate the objects, rep-resented by the optimal number of GMM clusters selected (6).It would therefore be meaningless to attempt a classification using all the 125 classes that appear in the SIMBAD database for our sources.
For training, we exclude sources in these underrepresented classes.
• Extended Large accretors with optical lines (Seyfert): Seyfert_1 and Seyfert_2 • Small accretors (XB): HMXB, LMXB, and XB • Young stars (YSO): YSO, TTau*, and Orion_V* In the context of this study, "accretors" refer to accreting compact objects.In "Large accretors", the central object is a super-massive black hole, which results in a high luminosity."Small accretors" involve a stellar remnant (a neutron star or black hole) accreting matter from a companion star.
For a significant number of X-ray sources, SIMBAD only provides a general class that is not informative as of specific physical properties ("Star", "X", "Radio", "IR", "Blue", "UV", "gamma", "PartofG").We assume objects in those classes to be unclassified, and assign them classes using our algorithm.
All in all, the total number of CSC detections for which we will provide a label is 14, 507.

Detection level classification
Using the probabilistic approach described in § 3.2 we have assigned a probabilistic class to all 14,507 CSC detections.
Figure 4 shows the True vs. Predicted confusion matrix for the benchmark set previously described.Each element of the matrix indicates the fraction of test objects with a given True class that have been correctly classified by the algorithm.At the detection level, our X-ray based classification method is successful at distinguishing between large accretors (QSO, AGN, and Seyferts), small accretors (X-ray binaries), and young stars (YSOs, T-Tauri stars, and Orion V* variables).However, a fair amount of confusion persists between sub-classes of each major class.
For example, A QSO is correctly classified as a such about 40% of the times, whereas about 30% of the times, the algorithm thinks it is a Seyfert 1 galaxy.In fact, QSO X-ray spectra can be similar to Seyfert 1 spectra (Dadina et al. 2016).Lacking a distance or a redshift, the algorithm cannot distinguish between the two in the detection level.Similarly, a HMXB gets correctly classified over 50% of the time, and incorrectly classified as a LMXB in about 8% of the cases.On the other hand, LMXBs are classified as HMXBs in 44% of the test examples.Among young stars, most confusion occurs for T-Tauri stars, which are classified as such only a third of the times, with the remaining examples classified in equal amounts as either YSOs or Orion V*.
Large accretors are most likely to be confused with HMXB, but that only happens about 20% of the time, and small accretors are almost never assigned a large accretor label, with the exception perhaps of the XB class, which is confused with QSOs in 13% of the examples.The most robust classification is provided for young stars, which are almost never classified in other groups.This is illustrated in Figure 5, where we have grouped individual classes into major classes A significant part of the confusion is likely due to differences in the spectral states of the different objects during different epochs.More robust classifications result from considering multiple detections of a source, but we note that single-epoch, X-ray-only classification of CSC sources is possible at some significance level.
How do X-ray properties distribute between classes? Figure 6 shows a density map of the distribution of hardness ratios for the classified detections.We observe relatively distinct distributions in the spectral shape for each of the different classes, even if there are also similarities and degeneracies.For example, QSOs, AGNs, and Seyfert 1 galaxies all appear to have a similar mean in their hardness ratio distributions, but they differ in the standard deviations of both hardness ratios.Seyfert 2 galaxies, on the other hand, have a slightly elevated hard-to-soft ratio with respect to the other three groups.X-ray binaries show wider distributions of their hardness ratios, with low-mass X-ray binaries even showing signs of multimodality.This is relevant in the distinction between large and small accretors.Young stars show overall softer spectra with respect to compact objects, although in the case of YSO there seems to be at least two separate populations with different spectral shapes.T-Tauri stars have consistently lower hardness ratios, whereas Orion variables, which tend to show eruptive behavior, appear to occupy a region in between the two populations of YSO types.
Figure 7 shows a similar density map, but this time for the variabilities in the broad (0.5keV-7.0keV) and soft (0.5keV-1.2keV) bands.We note significant variations in the level of soft-band variability across classes, with less variability in the integrated broad band.For most of the classes, variability is not correlated between bands, indicating spectral variability.Among large accretors, QSOs appear less variable in the soft band with respect to AGN, and Seyfert 2 galaxies appear more spectrally variable than their Seyfert 1 counterparts, which also show more variability in the broad band.For small accretors, no significant changes in variability is discernible between HMXBs and LMXBs, whereas sources classified plainly as XBs show a smaller degree of soft band variability.Finally, for young stars, T-Tauri stars are clearly differentiated by a marked bimodal distribution, and the majority of objects having a significant variability probability (var_prob_b ~1).The latter might be related to flaring events in that stage of their evolution.

Master classification
So far we have obtained probabilistic classifications for every detection of a source in our dataset.In order to obtain more accurate classes for a given CSC source, we can now operate on the probabilistic labels assigned to each detection of a source.Spectral variability, different exposure times or signal to noise level might translate into different classifications in different epochs.From the 8756 unique sources in the data, 6802 (77.7%) have only 1 detection (for these sources, the detection level classification is our best guess), 1000 (11.4%) have 2 detections, 882 (10.1%) have between 3 and 10 detections, and finally just 72 (0.8%) sources have more than 10 detections.For each source we now assign a master class using the procedure described in § 3.2.2.
Overall, both the hard and soft voting systems agree with each other: for a given source, the highest average probability across classes usually aligns with the primary type assigned to the majority of a source's detections.Only 485 sources (6% of the total) result in different classifications.The Cohen's Kappa coefficient (Cohen 1960), a measurement of agreement between the two methods, was calculated to be 0.94.We compile a table of uniquely classified and a table of ambiguous class sources depending on whether the class assigned agree between the two methods.We also provide the mean classification probabilities and their standard deviations.In the uniquely classified sources table, master_class refers to the class assigned from both the soft and hard classifier, and agg_master_class denotes its corresponding aggregated class.In the ambiguous sources table, soft_master_class and hard_master_class, refer respectively to each of the methods.We provide the number of detections per source as detection_count.Sample extracts of the uniquely classified and ambiguous classification tables, sorted by detection count, are presented in Table A1

Astrophysical validation of the classification
Certain types of CSC sources are expected to occupy specific regions of the sky.For example, we expect YSOs to be associated with starforming regions along the galactic plane, and AGNs and quasars to be mostly located off the galactic plane.As a way to validate our classification, we investigate whether this expectations are satisfied -in a statistical sense-for the objects that we have classified.Figure 8 shows density maps for each aggregated class over Mollweide projections.We note that objects assigned to the YSO class clearly concentrate along the galactic plane, while AGNs are mostly confined to the extragalactic area of the sky.Objects that we have classified as X-ray binaries have a tendency to be located along the galactic plane, but not as starkly as in the case of YSOs.Seyfert galacies are mostly extragalactic, but less markedly than AGNs.These maps provide an early indication that our classifier is associating the right objects to the right environments, despite not having access to any information about the object's location in the sky.It is important to note that Chandra is not an all-sky survey, and therefore these maps should not be interpreted as proxies for the actual surface densities of each class.This association between assigned classes and environment holds even at smaller scales.Figure 9 shows previously unclassified sources within the Orion Nebula region (M42) to which we have assigned a master class.Most of the sources in this area are classified as YSO, as expected in this well known star-forming region (O'dell 2001).We also observe a few sources classified differently.We now discuss two examples of these: • 2CXO J053516.0-052353 has 4 CSC detections and is classified as an AGN. Figure 10 shows its probability distribution across aggregated classes.There is significant confusion between classes, with the AGN and YSO being equally probable.In observation ID 4374, the source was classified as a YSO with a high probability of 0.99, and in observation ID 3744, it was classified as Orion_V* with a probability of 0.75.In the remaining two observations, the source was classified as AGN with probabilities around 0.7, which resulted in the final hard-voting class.Furthermore, it is worth noting that the standard deviation is higher for the YSO class (0.41) compared to the AGN class (0.38).Such high standard deviations are not uncommon in our classification scheme.In summary, this is a close call between YSO and AGN class, which would still be classified as a Young Star if only the aggregated classes defined in § 4.2.1 were used.
• 2CXO J053615.0-051530 has 5 CSC detections and is classified as Seyfert_1 galaxy, with relatively high probabilities also assigned to the QSO and AGN classes.The closest match in the SIMBAD database, source V* V828 Ori, is as an Orion Variable located about 3 arcsec to the southwest of the CSC source.In Szegedi-Elek et al. ( 2013), the source is included in the TTau* group, and other general variable and proper motion stars catalogs.We cannot rule out the possibility that the CSC detection might be in fact a background AGN.
We also validate membership for objects classified as X-ray binaries.Source 2CXO J193309.5+185902, for example, is classified as a HMXB by our method with just one detection and a probability of 95%.This is consistent with other works such as Yang et al. (2022), in which the source is classified as a HMXB with a probability of 69 ± 14%.More generally, we can look for XB in regions where such type of X-ray source is expected, such as the spiral arms of galaxies.We explored source aggregated classifications over the region of the M33 galaxy, in the outskirts of which sources 2CXO J013301  machine-learning based classifications, such as Tranin et al. (2022).Some sources are also missclassified as X-ray binaries, based on comparison with independent classification.Source 2CXO J053747.4-691019, for example, is confidently classified as a HMXB by our algorithm.However, it has been classified as pulsar PSR J0537-6910 located in the Large Magellanic Cloud (Lin et al. 2012).We note, however, that in some X-ray binaries the compact object can be a pulsar (Wĳnands & Van Der Klis 1998).
Source 2CXO J031948.1+413046 is located at the heart of the Perseus galaxy cluster.It is classified very confidently (>0.9 probability in all three detections) as a Seyfert_2 galaxy candidate.Other sources in this core region have been studied and suggested as either Type 1.5 or Type 2 Seyferts, such as 3C 84 Rani et al. (2018); Véron-Cetty, M.-P.& Véron, P. ( 2006).These specific galaxies typically exhibit a brightly illuminated core, predominantly at infrared wavelengths.However, the number of Seyfert 2 galaxies in the Perseus cluster is limited, with (Matt, G. et al. 2012) being one of the previously identified examples.
A different test astrophysical environment is provided by NGC 4649 and its surroundings.This is an elliptical galaxy located in the Virgo cluster with a population of X-ray sources detected by Chandra, the vast majority of which has been classified as Low-Mass X-ray Binaries (Luo et al. 2013;D'Abrusco et al. 2014).This classification is based on their overall X-ray properties, density, and spatial coincidence with a rich population of Globular Clusters detected in Hubble-Space Telescope images (Strader et al. 2012).All of the previously unclassified CSC sources near NGC 4649 for which we are providing a master class are assigned a XB label by our pipeline, as shown in Figure 11.The CSC source closest to the X-ray source identified by Luo et al. (2013) as the AGN in the center of NGC4649 is 2CXO J124339.9+113309.As a previously classified source, it was not assigned a label in our classification, but it is located in the cluster 0, under the GinPair (Galaxy in Pair of Galaxies) category from SIMBAD.The next closest source (2CXO J124339.8+113307) is classified as an X-ray Binary in each of the two detections included in the CSC.

DISCUSSION
The advent of all-sky, time-domain X-ray surveying missions such as eROSITA, and the proposed capabilities of upcoming probe and flagship facilities such as the Advanced X-ray Imaging Satellite (AXIS), Arcus, and eventually Lynx, underline the importance of a robust method for the classification of X-ray sources.This is important not only because of the increasing volume of high-energy sources detected, but also for the design of multi-wavelength follow-up campaigns of particular objects of interest, in particular those that inform models of transient behavior and extreme accretion.Classification attempts are usually limited by two factors: the lack of a robust and balanced training set across X-ray classes, and the fact that a significant number of X-ray sources lack optical or infrared counterparts that could provide useful insights for the categorization of sources.
Leveraging the information contained in the Chandra Source Catalog, we have designed an unsupervised classification algorithm that acknowledges those limitations, with the scope of investigating to what extent the X-ray properties alone can provide relevant insights for the labeling of X-ray sources.We now discuss our findings, focusing on four major aspects: i) the performance of our method with respect to existing classification approaches and curated catalogs; ii) the implications for astrophysics of the relationship between specific classes and the distribution of X-ray properties; iii) the impact of multiple detections of a source on classification, and iv) the limitations of our method.

Comparison to other methods and catalogs
We compare our classification results with the following independent methods and compiled catalogs: • The training dataset described in Yang et al. (2022) in the context of supervised classification, containing 2947 literature verified CSC X-ray sources classified as AGN, Cataclismic Variable (CV), highmass star, HMXB, low-mass star, LMXB, neutron star, and YSO.(MUWCLASS TD) • The catalog of confidently classified CSC sources presented from the Yang et al. (2022) pipeline, which contains 31000 CSC sources, with the same labels as the previous items.(MUWCLASS CCGCS) • The machine learning classification approach of XMM-Newton 4XMM-DR10 catalog sources, presented in Tranin et al. (2022).This method assigns four different labels: AGN, stars, X-ray binaries (XRBs), and cataclysmic variables (CVs).(4XMM-DR10 Classification) • The catalog of 224168 galaxies presented in Toba et al. ( 2014) based on WISE and SDSS data.(WISE-SDSS Galaxies) • The catalog of 121 new redshift obscured AGNs of the Chandra Source Catalog presented in Sicilian et al. (2022) (Obscured AGNs).
To perform the comparison, we cross-matched our catalog of uniquely classified sources with each of the reference catalogs, either using a coordinate radius (≤1 ′′ ) or matching exact CSC names.Where appropriate, we merge similar classes into a single class prior to comparison.For example, we merge the AGN and Seyfert aggregated classes into a AGN+Seyfert class.We also merge our LMXB and HMXB classes into a single XB class.Table 4 summarizes the results of this comparison.The table lists for each case the number of matches, the precision, and the recall.The Precision column refers to the precision when all labels in the reference catalog are included, whereas the R Precision columns refers to the precision when only reference catalog labels that are in both catalogs are considered.
We note a significant increase in the precision when only the restricted set of labels is considered.In particular, 98% of the objects that we classify as AGN+Seyfert actually belong to that class according to (Yang et al. 2022), with precision in other classes ranging from 6% to about 80%.The lowest precision is obtained for X-ray binaries, for which only 6% of the objects that we classified as XBs actually belong to that class according to the reference catalog.The recall numbers tend to be more uniform across classes, with percentages ranging from 20% to 88%.For example, of all AGN+Seyfert objects in the (Yang et al. 2022) training set, we correctly classify as such 83% of them, the same proportion as for YSOs.The recall for the AGN+Seyfert class increases to 88% when the WISE-SDSS catalog is considered.We observe the worst recall performance when trying to correctly classify obscured AGNs.From the (Sicilian et al. 2022) catalog, we have only classified about 20% of their objects as AGNs or Seyfert.
In general, most classifications we labelled as AGN+Seyfert are accurate with respect to the MUWCLASS set, and we also recover a fair amount of objects actually belonging to that classes.X-ray binaries show comparatively lower precision and recall.We can recover about 50% of them from control catalogs, but we tend to assign XB labels to objects that are classified in other categories (mostly AGN) in control catalogs.Spectral similarities between these two types of accreting objects is a likely reason for this degeneracy, which can be broken if the distance to the sources is known (Volonteri et al. 2017).A significant fraction of the objects that we classify as YSOs are classified as low-mass stars in the MUWCLASS training set.62% .hard_hs vs. hard_hm distribution density estimation for each class group, visualized using a kernel density estimate (KDE).Darker hues represent areas of higher concentration within the distribution.Each corresponding aggregated class density is included, delineated by gray contours and referenced in the gray boxes.Note that the density estimation may appear to exceed the valid value boundaries due to the smoothing effect.This effect does not indicate actual data points outside the hardness ratio limits but results from the bandwidth choice in the estimation.var_prob_s distribution density estimation for each class group, visualized using a kernel density estimate (KDE).Darker hues represent areas of higher concentration within the distribution.Each corresponding aggregated class density is included, delineated by gray contours and referenced in the gray boxes.Note that the density estimation may appear to exceed the valid value boundaries due to the smoothing effect.This effect does not indicate actual data points outside the probability limits but results from the bandwidth choice in the estimation.
of the objects classified as stars in MUWCLASS, on the other hand, have YSO labels in SIMBAD.
We find 2006 matching sources with the 4XMM-DR10 classification catalog.Almost 80% of the sources that we classify as AGN or Seyfert galaxies agree with the Tranin et al. (2022) probabilistic classification, and of all AGNs and Seyfert galaxies in the 4XMM-DR10 catalog, we recover 62% of them with the correct class.Again, X-ray binaries show lower precision and recall, likely due to confusion with supermassive accretors.Of the 2006 matched sources, 435 were categorized as YSO in our classification scheme.The majority of these sources were classified either as Star or XRB by the 4XMM-DR10 classification.We obtain high precision and recall for optically or IR-selected galaxies in Toba et al. (2014), although these are based on a very small number of matches.

Accretion-powered sources
One remarkable result is that when we limit ourselves to the classes assigned by our algorithm, the confusion between small accretors and large accretors is not widespread.Such confusion would be expected on the basis of similar properties due to the accretion nature of the X-ray emission (Padovani et al. 2017).With certain degree of confidence, however, in our analysis we can distinguish X-ray binaries from AGNs, QSOs and Seyferts when only the X-ray properties are considered.As Figure 5 shows, large accretors are correctly classified as such in 75% of the cases, and incorrectly classified as X-ray binaries only on 19% of the cases.On the other hand, small accretors are correctly classified as X-ray binaries in 66% of the cases, and incorrectly classified as large accretors in 20% of the cases.
The feature distribution maps of Figure 6 and Figure 7 suggest  that the differentiation is possible based on a wider range of hardness ratios in X-ray binaries compared to AGNs and Seyferts, as well as a slightly increased average hardness for X-ray binaries, in particular HMXBs.The broader range of hardness ratios can be interpreted in terms of the different spectral states of X-ray binaries, which occur in much shorter timescales than in AGNs (Remillard & McClintock 2006).Within the large accretors, AGN and Seyfert galaxies are often confused, which is likely the result of a nomenclature difference, Seyferts being a particular type of AGN (Peterson 1997).We note, however, that Seyferts show more variability than AGNs, and Seyfert 2 galaxies in particular are more likely to be variable in the soft band, probably due to more obscuration due to their orientation with respect to the line of sight, in the context of the AGN unified model.QSOs in our sample, on the other hand, have tighter constrains in their hardness ratios and variabilities, consistent with the absence of spectral variability.In particular, they are less variable in the soft band with respect to AGNs and Seyferts.Stochastic variability in the thermal soft band of AGN spectra is usually related to disk instabilities, jet activity, or spottiness of infalling matter (Jovanović & Popović 2009), all of which might be more stable at the the high accretion rates seen in QSOs, provided that they do not have a blazarlike orientation.Hard band variability, on the other hand, relates to the hot corona (Ballantyne & Xiang 2020;Petrucci et al. 2001;Haardt & Maraschi 1991;Galeev et al. 1979), and appears more prominent in Seyfert 1 galaxies, for which the line of sight intersects the corona more directly (Soldi et al. 2014).Heavy absorption by the torus can also enhance both soft band variability the hardness ratios in Seyfert 2 galaxies due to preferential absorption of soft photons (Turner et al. 1997; Risaliti et al. 1999).The results in Figure 6 and Figure 7 show overall agreement with this unified model.Compared to QSOs and AGNs, objects classified as X-ray binaries show a larger range of variability in the broad band, resembling the patterns observed in Seyfert galaxies, with a slightly increased broad band variability.Resolved non-thermal coronal X-ray emission is likely to be at the root of this similarity between systems of significantly different masses and physical scales.Recent studies have provided evidence of similar spectral state evolution in the luminosity-hardness diagram, related to the accretion processes in both XBs and AGNs (Fernández-Ontiveros & Muñoz-Darias 2021).No dusty torus is present in X-ray binaries, however.This might explain why in Figure 7 both LMXBs and HMXBs tend to resemble more the variability pattern of Seyfert 1 galaxies, for which the line of sight intersects less obscuring material, and more of the corona.
We recognize that the incorporation of multi-wavelength counterparts can improve the separation between AGN-type sources and X-ray binaries, in particular if sources can be associated to optically derived redshifts, as has been shown for example in Yang et al. (2022).Association with host galaxies can also help in breaking degeneracies between these two types of accretors.We point out, however, that only a fraction of X-ray sources have optical and IR counterparts.In particular, less than half of the CSC sources have at least one optical or IR counterpart when a cross-match is performed using the Gaia DR3, Legacy Survey DR10, PanSTARRS-1, and 2MASS catalogs4 .This highlights the importance of assessing the quality of the classification with only X-ray information is used.

Young stars
Objects classified as young stars have consistently higher variabilities, in particular in the broad band.T-Tauri stars are the only type of object classified by our pipeline for which the majority of examples have high variability both in the broand an soft bands.High levels of variability may be tied to irregularities in the accretion processes occurring within these young stellar objects (Testa 2010;Preibisch et al. 2005).In addition, it could be linked to strong magnetic fields, which could lead to prominent magnetic activity such as flares and coronal mass ejections (Testa 2010).This high level of variability observed is consistent with our understanding of the the processes occurring at this stage of stellar evolution.While a similar level of magnetic activity is expected from Orion V* stars, we do not observe it in our sample.
Objects classified as YSOs also show a bimodal distribution in their hardness ratios.Hard X-ray emission during coronal events might be partly responsible for this behavior, but degeneracy with other types of hard-emitters, such as HMXBs, can also be associated to the bimodal distribution.In fact, missclassification of other types of objects as YSO (the class with the highest classification accuracy in validation) is not uncommon, although lower that you would expect from the limited amount of information encoded in the X-rays.For example, source 2CXO J033829.0-352701,located at the heart of the weakly active galaxy NGC 1399 within the Fornax Cluster, was classified very confidently as a YSO by our method, with a mean probability of 1.0 across six detections.It has been classified as an elliptical galaxy (De Vaucouleurs et al. 1991) with a Seyfert 2 AGN (Véron-Cetty & Véron 2006).The source is soft (hard_hs < −0.5) and shows remarkably low variability probability (var_prob_b < 0.2).While Seyfert 2 galaxies tend to be associated with harder spectra (the line of view intersects the obscuring torus that preferentially absorbs soft-X-ray photons), the specific orientation of the system with respect to the line of sight, as well as different contributions from scattering and reflection, can result in a softer spectrum.
We examined the properties of source detections with a YSO label in SIMBAD, as well as those labelled as YSOs in our classification, and found that the density distributions of properties are consistent for both the reference set and the classified set.There are, however, some difference in certain regions of the parameter space that could lead to misclassifications, because the reference set is not fully representative of the population.

Differences in the hard and soft voting classifiers
Our catalog of ambiguously classified objects comprises sources whose hard and soft classification disagree.That is, for these sources, the most common class among detections is not the same class with the highest average probability.These are cases in which there is ambiguity between two different classes of similar types, such as AGNs and Seyferts.We investigate some of these objects in order to gain additional astrophysical insight about how the algorithm asigns classes.
Source 2CXO J004228.2+411222 has 66 detections and is classified as a LMXB by the hard voting classifier, and as HMXB by the soft voting classifier.Independent work has classified it as an X-ray binary/black hole candidate in M31 (Arnason et al. 2020;Barnard et al. 2014Barnard et al. , 2013)).On the other hand, 2CXO J004246.9+411615, which has also been classified as an X-ray binary or black hole candidate in M31 (Barnard et al. 2014(Barnard et al. , 2013) ) and has 22 detections, is classified as a Seyfert_1 by the hard voting classifier, and as a LMXB by the soft voting classifier.In this latter case the mean probabilities are 0.22 for the Seyfert_1 class and 0.23 for the LMXB) class.So, whereas in the case of 2CXO J004228.2+411222 the ambiguity can be interpreted in terms of a difference of spectral hardness between two types of X-ray binaries and our method still assigns a general class that agrees with previous literature, in the latter case the difference in probability between the two candidate classes is negligible, and the ambiguity is unlikely to be resolved without additional information, such as the redshift of the intrinsic luminosity of the source.Such differences between soft and hard classification should be considered when interpreting the results of our classification catalog.
The X-ray pulsar magnetar 2CXO J010043.0-721133(Durant & van Kerkwĳk 2005, 2008) is classified as an Orion_V* by the hard voting classifier, and as a YSO by the soft voting classifier with 14 detections.Both assigned classes fall in the aggregated class YSO.While the detection could be consistent with a pulsar with a fading optical counterpart (Durant & van Kerkwĳk 2005, 2008), the reason that it got assigned a YSO class is that the class "pulsar" is not among the classes that we have selected to be assigned to the classified objects.Young stars and pulsars might share some of their properties in the space of X-ray features.
Finally, source 2CXO J020938.5-100847 is classified as a TTau* by the hard voting classifier, and as an AGN by the soft voting classifier.This source only has 2 detections.Mean probabilities are 0.19 for TTau* and 0.44 for AGN, which favor the   label.In fact, the source is located less than 2 arcseconds from the nuclear region of galaxy NGC 838 in the group HCG 16 (O'Sullivan et al. 2014), which has been spectroscopically classified as an active galaxy (LINER type).In this case, the mean probabilities are more informative of the true class of the objects.Overall, if the three aggregated classes (AGN+Seyfert, XB, YSO) are used, 46% of the sources in the list of ambiguously classified sources show agreement between the two voting systems.This suggests that even when there is disagreement between the two voting systems for the fine classification between the labels in Figure 4, a significant fraction of these ambiguous sources can still be successfully classified in the more general classes of Figure 5. Taken all together, the probabilistic elements of our classification can provide insight into the nature of these objects, even in the absence of optical counterparts.

Limitations and future work
We have demonstrated that a unsupervised approach to the classification of X-ray sources is possible even in the absence of an optical confirmation.Despite its applicability to X-ray catalogs, the GMM approach that we have adopted has a number of limitations that primarily stem from two factors: a) the intrinsic degeneracy in the observational features of sources of different types, and b) uncertainties in the knowledge base used (i.e.SIMBAD) to associate cluster membership to specific classes.
Feature degeneracies can be mitigated by incorporating ancillary information about individual sources, such as environmental properties (galaxy host information, among others).While this is in principle possible for a fraction of X-ray sources, it remains unfeasible for X-ray sources without associated optical counterparts.Limitations on the knowledge base can be addressed, in principle, by constructing detailed training sets based on spectroscopic classification or other independent signatures.An example of this is the training set presented in Yang et al. (2022).However, it is extremely difficult to construct such training sets while ensuring that they remain representative, balanced, and unbiased.This is both due to class imbalance and lack of independent classifications for X-ray sources in the majority of cases.Furthermore, we refrain from using coordinates as features due the CSC's skewed and non-homogeneous sky coverage, which could introduce biases.We therefore embrace these limitations as a challenge to push the limits of X-ray only unsupervised classification.In our experiments, we found our method to be robust against class imbalance, outperforming traditional classification algorithms.To evaluate the impact of class imbalance and compare our approach's performance with standard supervised methods, we detail a classification exercise using Support Vector Machines in online Appendix B.
Gaussian Mixture Models are constrained by the assumption of Gaussian-distributed clusters, which does not hold for most distributions.Their parametric nature also restricts their ability to fit multidimensional data distributions in practical problems (Mallapragada et al. 2010).We used heuristic methods (such as the BIC and elbow methods) and identification of property distribution across multiple experiments to choose a number of Gaussian components that matches the overall behavior of the feature distributions.Alternative methods involve the use of Gaussian Mixtures with a Dirichlet process prior (Teh et al. 2010;Görür & Edward Rasmussen 2010), or other non-parametric proposals (Mallapragada et al. 2010).Also, spectral clustering is a generalization of classical clustering methods and is capable of identifying non-convex clusters (Hastie et al. 2009).
The pipeline proposed here is similar to the cluster-then-label semi-supervised method, wherein a clustering algorithm is first employed, followed by a supervised learner to label the unlabeled points in each cluster (Zhu & Goldberg 2009).In our case, the Mahalanobis distance acts as the supervised learner, although each cluster can contain multiple classes, which justifies our probabilistic approach.A recurrent version of this model allows for the rapid classification of new detections, which would allow for the incorporation of the model to data processing pipelines of Chandra or other observatories.Using our method, sources with diverging hard and soft voting classifications could be used to recognize objects with significant spectral or flux variability between observations, such as X-ray binaries transitioning between states.
It should be noted that comparisons of probabilities between source detections in different clusters may not be informative for the purposes of comparing those two observations.That is because the probabilities presented in our results are relative to each cluster's individual space.For instance, the assigned probability of 100% for source 2CXO J020938.4-100846(as presented in Section 4) as an HMXB implies that within its specific cluster, the reference HMXB distribution is the closest to the source detection.
We emphasize that our pipeline is not only easily reproducible but also adaptable to other catalog classification endeavors.We invite interested readers to explore the provided GitHub repository and implement proposed modifications.This could help enhance the performance of the work or enable exploration of other applications, such as the classification of nature-specific types of objects (e.g., YSO and non-YSO classification).

SUMMARY AND CONCLUSIONS
The automatic classification of X-ray sources is a necessary step in extracting astrophysical information from compiled source catalogs.Classification is useful for the study of individual objects, statistics for population studies, as well as for anomaly detection, i.e., the identification of new unexplored phenomena, including transients and spectrally extreme sources.This is true for Chandra data, and for other X-ray missions that have compiled catalogs of X-ray sources and their properties, including XMM-Newton, eROSITA, etc.Despite the importance of this task, classification of astronomical X-ray sources remains challenging.This is due to mainly two issues: i) a significant fraction of X-ray sources lack optical or infrared counterparts that could provide useful ancillary information, such as redshift; and ii) the construction of reliable training sets that would allow automatic supervised learning classification is not trivial.These training sets are incomplete and unbalanced at best, due to the lack of spectroscopic confirmation for a significant number of sources.
Despite these difficulties, supervised classification of X-ray source catalogs has been attempted (Yang et al. 2022;Chen et al. 2023;Kumaran et al. 2023).Such efforts rely on the ability to obtain multiwavelength properties or observables for the sources to be classified, which is not always possible.In this work, we have developed an alternative methodology that employs unsupervised machine learning to provide probabilistic classes to Chandra Source Catalog sources.The method relies on clustering of X-ray properties such as hardness ratios and variability using Gaussian Mixtures.It suffers less from the lack of a reliable training set, as classes are assigned by association with cluster membership, as well as some topological considerations.Here are our main findings: (i) As an alternative to directly learning from a representative set of labelled multi-wavelength training sources, we have demonstrated that it is possible to assign probabilistic classes to X-ray sources by first clustering them according to the distribution of X-ray properties, and then using a metric (the Mahalanobis distance) within each cluster to evaluate their probabilistic classes based on their distance to a much less representative set of independently classified sources.
(ii) We provide a catalog of probabilistic classes for 8,756 sources, comprising a total of 14,507 detections.For each source, we provide a probability distribution over classes, as well as a measure of the uncertainty in its assigned class.We have validated our classification by using an internal validation set, and by comparing our classification with those performed using supervised approaches.
(iii) Our methodology is particularly successful at identifying emission from young stellar objects, with an 88% overall accuracy.When compared with existing classification catalogs, our classification can reach a level agreement above 90% in identifying AGNs and Seyfert galaxies.
(iv) Our method is able to distinguish between small and large compact accretors (that is, X-ray binaries and active galactic nuclei) with more than 50% confidence, despite the widespread assumption that such accreting systems can be very similar in their X-ray properties, and difficult to separate in two different classes in the absence of a redshift or distance measurement.This distinction is possible due to a wider range of hardness ratios in X-ray binaries compared to AGNs and Seyferts, along with a slightly higher average hardness for X-ray binaries, in particular HMXBs.
(v) Compared to QSOs and AGNs, objects classified as X-ray binaries show a larger range of variability in the broad band, resembling the patterns observed in Seyfert galaxies, with a slightly increased broad band variability.Resolved non-thermal coronal X-ray emission is likely to be at the root of this similarity between systems of significantly different masses and physical scales.
(vi) Differences in spectral variability and hardness ratios between objects classified as QSOs, Seyfert 1 and 2 galaxies, and other AGNs, are consistent with the unified AGN model.
(vii) Broad band variability stands as a key feature to differentiate XBs from AGNs, when relying solely on X-ray data.
(viii) Stellar variability due to coronal ejections in young stars is manifested in a distinct distribution of variabilities and hardness ratios in these objects, that informs their classification.Short-lived flares of hard X-ray photons are a likely cause of this distribution of the X-ray features.
(ix) Our unsupervised method is similar to a semi-supervised method, with a probabilistic metric (e.g. the Mahalanobis distance) replacing a classical supervised classifier.
Our complete results are accessible in the online supplementary material.Reproducibility is facilitated through the code available in our GitHub repository: https://github.com/samuelperezdi/umlcaxs.Additionally, we provide a notebook for classifying individual X-ray sources5 .We offer a web app playground for direct data interaction: https://umlcaxs-playground.streamlit.app/.
Here, users can access visualizations of source positions in the sky and their properties, along with respective classifications.

Figure 1 .
Figure1.The Bayesian Information Criterion (BIC) (top) and the gradient of the BIC (bottom) as a function of the number of components .The number of components ranges from 2 to 20.The gray region delineates the confidence interval as determined by the standard deviation of each 's iteration results.The BIC function is smoothly decreasing, while its gradient shows a constant behaviour for values greater than  = 6.The red dashed vertical line highlights the function point at  = 6, which is the number of components that shows a better configuration in this technique.

Figure 2 .Figure 3 .
Figure 2. Source detections in galactic coordinates over a Mollweide projection, discriminated by their assigned cluster in colors and markers.We see a trend of extragalactic or galactic for points in particular clusters.

Figure 4 .
Figure 4. True vs. Predicted confusion matrix for the benchmark set using individual detections only, normalized by row.For each true class the proportion of source detections of that class assigned to each of the possible labels is shown.

Figure 5 .
Figure 5. True vs. Predicted confusion matrix for the benchmark set, normalized by row.Classes were replaced by the aggregated classes AGN, Seyfert, QSO, YSO.The matrix displays the proportion of source detections in a specific class that were correctly classified or misclassified.The data used in this plot is not a final result of our pipeline.

Figure 6
Figure6.hard_hs vs. hard_hm distribution density estimation for each class group, visualized using a kernel density estimate (KDE).Darker hues represent areas of higher concentration within the distribution.Each corresponding aggregated class density is included, delineated by gray contours and referenced in the gray boxes.Note that the density estimation may appear to exceed the valid value boundaries due to the smoothing effect.This effect does not indicate actual data points outside the hardness ratio limits but results from the bandwidth choice in the estimation.

Figure 7 .
Figure7.var_prob_b vs. var_prob_s distribution density estimation for each class group, visualized using a kernel density estimate (KDE).Darker hues represent areas of higher concentration within the distribution.Each corresponding aggregated class density is included, delineated by gray contours and referenced in the gray boxes.Note that the density estimation may appear to exceed the valid value boundaries due to the smoothing effect.This effect does not indicate actual data points outside the probability limits but results from the bandwidth choice in the estimation.

Figure 8 .
Figure 8. Density map for the uniquely classified CSC sources that lacked a label before our study, in galactic coordinates, over Mollweide projections for each aggregated master class.

Figure 9 .Figure 10 .Figure 11 .
Figure 9.A set of previously unclassified sources in M42, overlaid on an AllWISE W1 band image in grayscale.The colors indicate different classes assigned by our algorithm.

Table 1 .
Properties selected for our analysis.Properties denoted with '_*' represent three features for energy bands: b (broad), h (hard), s (soft).Adapted from CSC documentation webpage.
var_newq_bproportion of the average of minimum and maximum count rates (i.e., data points in the light curve) during an observation relative to the mean count rate.var_max_b + var_min_b 2 • var_mean_b

Table 2 .
Preprocessing approach for each property.Energy bands are denoted by * = b, h, s (broad, hard, and soft).
and TableA2respectively. Full versions of the tables are available in the online supplementary material and the provided GitHub repository.

Table 3 .
The 10 most prevalent source classes in each of the clusters, ranked by number of examples.Source detections with non matching types are labeled as NaN.
.0+304043 and 2CXO J013324.4+304402 are located.Both sources are classified as XB, in agreement with other

Table 4 .
A summary of the comparison of our classification with previous work and catalogs.The number of matches corresponds to the number of sources belonging to our assigned class resulting from crossmatching the Chandra Source Catalog names or a matching criteria of ≤1 ′′ .