Stellar infrared excesses can indicate various phenomena of interest from protoplanetary disks to debris disks or more speculatively technosignatures along the lines of Dyson spheres. In this paper we conduct a large search for extreme infrared excesses designed as a datadriven contextual anomaly detection pipeline. We focus our search on FGK stars close to the main sequence to favor nonyoung host stars. We look for excess in the midinfrared unlocking a large sample to search in while favoring extreme IR excess akin to the ones produced by extreme debris disks EDDs andor planetary collision events. We combine observations from ESA Gaia Data Release 3 the Two Micron AllSky Survey and the unWISE version of NASAs Widefield Infrared Survey Explorer WISE and create a catalog of 4898812 stars with G lt 16 mag. We consider a star to have an excess if it is substantially brighter in the W1 and W2 bands than what is predicted from an ensemble of machine learning models trained on the data taking optical and nearinfrared information as input features. We apply a set of additional cuts derived from the machine learning models and the objects astronomical features to avoid false positives and identify a set of 53 objects including one previously identified EDD candidate. The typical infraredexcess fractional luminosities we find are in the range 0.0050.1 consistent with previous EDD candidates and potential planetary collision events.
['Debris disks', 'Infrared excess', 'Outlier detection', 'Astronomy databases', '363', '788', '1934', '83', 'Astrophysics - Earth and Planetary Astrophysics', 'Astrophysics - Instrumentation and Methods for Astrophysics', 'Astrophysics - Solar and Stellar Astrophysics']
A Datadriven Search For Midinfrared Excesses Among Five Million Mainsequence FGK Stars
{'No Header': 'Typeset using L A T E X twocolumn style in AASTeX631', 'A Data-Driven Search For Mid-Infrared Excesses Among Five Million Main-Sequence FGK Stars': 'Gabriella Contardo 1 and David W. Hogg 2, 3, 4  \n1 Theoretical and Scientific Data Science, Scuola Internazionale Superiore di Studi Avanzati (SISSA), Trieste, Italy 2 Center for Cosmology and Particle Physics, Department of Physics, New York University, New York NY, USA 3 Center for Computational Astrophysics, Flatiron Institute, New York NY, USA 4 Max-Planck-Institut fur Astronomie, Heidelberg, Germany', 'ABSTRACT': "Stellar infrared excesses can indicate various phenomena of interest, from protoplanetary disks to debris disks, or (more speculatively) techno-signatures along the lines of Dyson spheres. In this paper, we conduct a large search for 'extreme' infrared excesses, designed as a data-driven contextual anomaly detection pipeline. We focus our search on FGK stars close to the main sequence to favor non-young host stars. We look for excess in the mid-infrared, unlocking a large sample to search in while favoring extreme IR excess akin to the ones produced by Extreme Debris Disks (EDD) and/or planetary collision events. We combine observations from ESA Gaia DR3, 2MASS , and the unWISE version of NASA WISE , and create a catalog of 4,898,812 stars with G < 16 mag. We consider a star to have an excess if it is substantially brighter in W 1 and W 2 bands than what is predicted from an ensemble of machinelearning models trained on the data, taking optical and near-infrared information as input features. We apply a set of additional cuts (derived from the ML models and the objects' astronomical features) to avoid false positives and identify a set of 53 objects, including one previously identified EDD candidate. Typical infrared-excess fractional luminosities we find are in the range 0 . 005 to 0 . 1, consistent with previous EDDs candidates and potential planetary collision events.", '1. INTRODUCTION': "The presence of infrared excess in stellar objects is, in general, a sign of circumstellar dust, protoplanetary disks, or debris disks. The strength of IR excesses is expected to decrease as a function of stellar age as a star's disk evolves (Wyatt et al. 2015), tying it to the evolutionary process in the star, its potential planets and other planetesimals or asteroid belts. As such, strong IR excesses are often associated with (relatively) young stellar objects. Cotten & Song (2016) detect IR excess in main-sequence stars, with a range of fractional luminosity L IR /L ∗ between 10 -5 to 10 -2 , averaging around 10 -4 , and with a single candidate having a fractional luminosity ≃ 0 . 1. However, several studies have also highlighted the existence of high IR excess (with fractional luminosities around 0.02 and up to 0.1) for stellar objects with wider age ranges (Mo'or  \ngabriella.contardo@sissa.it  \net al. 2021; Balog et al. 2009; Kennedy & Wyatt 2013; Song et al. 2005; Zuckerman et al. 2012; Melis et al. 2021; Rhee et al. 2008), up to Gyrs. These fractional luminosities are in between the ones expected for (regular) debris disks ( L IR /L ∗ < 10 -3 ) and protostars ( L IR /L ∗ > 0 . 1) (Uzpen et al. 2005). Those excesses are currently thought to come from planetary or planetesimal collisions and often referred to as Extreme Debris Disks (EDDs). However, the stellar ages of some of these candidates are in tension with the current models of rocky planet formation which predict that most collision events (leading to such disks) should occur within the first 100 Myr. Therefore, these objects might indicate that these processes last longer than thought or that other processes lead to such disks.  \nKenworthy et al. (2023) recently identified a Sun-like object of 300 Myr undergoing a sharp increase in infrared, reaching a fractional luminosity similar to the ones stated above (0 . 04), followed by an obscuration in the optical some 1,000 days later. The authors show that  \nthe observations are consistent with a planetary collision between two exoplanets of several to tens of Earth masses at 2-16 astronomical units from the central star. Another Sun-like star, of ∼ 600Myr, with a detected IR-excess of 0 . 01 fractional luminosity with an apparent dimming in the optical was identified in Melis et al. (2021). Interestingly, there is also one account of a fast disappearance of the IR-excess (Melis et al. 2012).  \nCurrently, there are very few known EDD candidates: for instance, there were 17 reported in Mo'or et al. (2021). Those events are a relatively rare occurrence: Mo'or et al. (2021) searched through 78 , 650 objects and found 8 candidates -6 new and 2 known, i.e. an 'occurrence' 1 rate of 0.01%. A similar occurence rate was suggested by Kennedy & Wyatt (2013). Those extreme excesses thus likely come from rare or potentially very short-lived events. It is critical to increase our sample of candidates to better understand the underlying physical causes of these extreme IR excesses.  \nInterestingly, another possible-albeit less likelycauses of extreme IR excess are techno-signatures such as Dyson Spheres (hypothetical megastructures built around stars to collect energy). For instance, Suazo et al. (2022) provided upper limits on the prevalence of partial Dyson Spheres combining data from Gaia DR2 and AllWISE , using models of the resulting IR-excess (and optical obscuration) of different types (stages and temperature) of Dyson Spheres. Following on this work, Suazo et al. (2024) recently presented a search targeting such IR-excess types relying on these Dyson spheres models, with 7 candidates identified. The search for (extreme) contextual outliers in the IR is a nice example where the search for extraterrestrial intelligence (SETI) field intersects and can challenge the boundaries of current astrophysical models (in this instance, of planetary formation and evolution). As one would need to investigate and eliminate all other explanations before concluding an extreme IR-excess source is a Dyson Sphere, searches for 'Dysonian' signatures will potentially provide interesting objects for 'mainstream' astrophysics. This aspect has been coined as the 'ancillary benefit', one of the axes of merits one can evaluate technosignature searches by, defined by Sheikh (2020).  \nWe propose here to conduct a large search for such 'extreme' infrared excesses. We focus on a population of photometrically main-sequence stars, to favor poten-  \ntially non-young or transitional candidates. We also restrict our search to FGK stars, focusing on Sun-like systems: as this work is partially motivated by the search for techno-signatures, as well as potential planetary collision of the like found in Melis et al. (2010); Kenworthy et al. (2023), we want to focus on systems similar to ours. Additionally, reducing the variety of stellar types is also motivated by the nature of our approach. A larger diversity of stellar types would require our models to learn this and could lead to biases in our detection in terms of stellar types. However, the following methodology could and should be explored on other stellar populations.  \nWe choose to conduct our search in the mid-IR, and more specifically we look for excess in the WISE W 1 and W 2 bands: It has been shown that IR-excess typical of EDDs can be hot enough to show in W 1 and W 2, and we expect (possibly extremely naively) to see techno-signatures near the water triple-point. Restricting to these two WISE bands also unlocks a larger sample to investigate, as the sensitivity and reliability of WISE is much higher in those bands than in the W 3 and W 4 bands: our curated dataset, presented in Section 2, has 4.9 million objects. On the other hand, looking for excess in W 1 and W 2 means that we target only extreme events leaking in those bands, which can nonetheless be subtle, especially in the W 1 band.  \nTo approach this problem in a way that scales to our catalog, we present a data-driven pipeline. Our approach alleviates the need for good (computationally expensive) stellar modeling or fitting: we are looking instead for deviations (or contextual anomalies) according to the data , i.e., our sample of main-sequence FGK stars. That is, instead of using stellar modeling or approximations, we train a machine learning model to predict the expected W 1 and W 2 magnitudes from other (optical and near-infrared) photometric features. We then use this model to find stars with observed magnitudes in W 1 and W 2 higher than predicted (similar to what one would do with regular stellar modeling). Additionally, we deploy a set of measures to focus our search on stars that have consistently wrong predictions when using models trained on different subparts of the data. Details on our methodology are provided in Section 3. This approach is also computationally interesting compared to 'proper' SED fitting to scale to our 4.9 million objects dataset.  \nFrom a methodological perspective, this paper presents a pipeline to conduct targeted, contextual anomaly detection (see Chandola et al. (2009) for an anomalydetection taxonomy overview), benefiting from robust supervised machine learning methods (here Random Forest). Anomaly detection (in astronomy and astro-  \nphysics) is often framed as finding rare objects, i.e., finding data examples with low probability (in low-density regions) according to the data. To tackle this problem with ML methods, one can (for instance) get a density estimate, and retrieve the objects with low scores. Some caveats of this generic approach are that (i) obtaining reliable density estimates is non-trivial, especially as the number of dimensions increases, and (ii) such an approach will retrieve all (potentially) rare/odd objects, but not necessarily the 'interesting ones' (which motivates the growing literature on integrating activelearning in anomaly detection pipeline in astronomy such as in Ishida et al. (2021) and Lochner & Bassett (2021)). However, it is a good way to potentially uncover the infamous 'unknown unknowns'. On the other hand, in many cases, we might be looking for 'known unknowns' (as in well-specified unknowns, like here, where we want to look for specific types of anomalies: stars with IR-excess). 'Knowing' our unknowns, however, does not necessarily mean that we can turn this problem into a supervised classification one (where our detector would tell us a binary answer or a score probability of our object being an anomaly): We might not have actual examples of these anomalies, or only very few examples.  \nHere we propose to benefit from the fact that we are looking for targeted anomalies that deviate in a specific way (excess), in a specific region of the feature space (mid-infrared), considering their 'context' (their optical and near-IR observations). We will be looking for objects that are confidently wrong in the prediction of that part of the feature space (the MIR-bands) using the rest of the feature space (mostly optical and near-IR photometry, see Section 3 for details). This process is indirectly doing a form of density estimation: as we assume our model is flexible enough and generalizes well, objects with a high, confident error, deviate from the training set distribution in some significant ways, that are not predictable from their 'context' (input features). We note however that if the objects we are looking for were very common in our sample and with an excess predictable from the input features, our pipeline would not flag them as anomalies.  \nOne amusing thing about taking a data-driven approach is that we require less, and less accurate, knowledge about how stars work. Most infrared excess projects will have to deal with dust reddening and have physically accurate models for the multi-wavelength emission from stars. We require neither; we require only that infrared excess events be around stars that are typical in the  \nsense that there are (many) other stars in the target sample that are similar but lacking an infrared excess.  \nThe rest of this paper is organized as follows: the data curation is described in Section 2. The method details are provided in Section 3 along with a list of subsequent cuts applied to obtain our final set of 53 candidates, which include one previously detected EDD candidate (see Appendix D). Further analysis of these candidates including checking for H α emission, variability, Black Body fitting and their parameters, and rate of recovery for those, can be found in Section 4. We provide concluding remarks in Section 5. We provide in the Appendix A a summary Table of our candidates, and in Appendix C a check for mid-IR deficit. The data products as well as the pipeline code and the analysis code can be found at https://github.com/contardog/ NotATechnosignatureSearch.", '2. DATA': "We combine data from Gaia DR3 (Gaia Collaboration et al. 2016, 2023; Babusiaux et al. 2023), the Two Micron All Sky Survey ( 2MASS ) (Skrutskie et al. 2006) (Skrutskie et al. 2006), the unWISE catalogue (Schlafly et al. 2019) and AllWISE catalogue (Cutri et al. 2021). We describe in this section our data curation and selection pipeline.  \nUsing Gaia DR3 observations, we select stars with a G -band magnitude brighter than G < 16 mag . We apply a temperature cut s.t. the effective temperature is between 4000 and 7000K (restricting the stellar types to FGK), using effective temperatures from Gaia DR3 ( teff gspphot , inferred by GSP-Phot Aeneas from BP/RP spectra, apparent G magnitude and parallax). We impose the parallax over error to be over 10, ruwe < 1 . 4 and require that ebpminrp gspphot (reddening value using BP/RP spectra) is not null. We rely on the cross-match between Gaia DR3 and 2MASS ( gaiadr3.tmass psc xsc best neighbour ), and between Gaia DR3 and AllWISE ( gaiadr3.allwise best neighbour ), provided by the ESA, to get the counter-part candidates in each catalog. Using the DPAC curated catalogs available on ESA Gaia archive, we impose that the contamination flags ( cc flags ) from AllWISE is '0000', i.e. there is no (flagged) contamination in any of the bands. We also restrict our sample to enforce that the photometric uncertainty in 2MASS is under 0 . 2 for all bands, and with a 2MASS photometric quality flag of 'AAA' (i.e. SNR > 10). Additional cuts to prevent potential mismatch or contamination are also applied: we enforce that the  \n'number of mates' in both 2MASS and AllWISE is 0 (i.e. there is no other Gaia sources that have this object as best neighbor as well: this should ensure (to some extent) that we do not get sources that were actually more than 1 object but unresolved in 2MASS or WISE ), and that the 'number of neighbors' is equal to 1 (i.e. there is no additional source in the external catalog that matches the Gaia source within position errors). Both columns number of mates/neighbors are provided in the 'best neighbour' cross-match tables with Gaia DR3 from ESA for both 2MASS and AllWISE .  \nNote that those cuts do not necessarily prevent potential contamination from close nearby sources if both sources were resolved in all catalogs (we suggest an additional crowding flag to address this in Section 3). Finally, we impose that the angular distance between the Gaia source and its best neighbors in the external catalogs is under 0 . 15 arcseconds. The corresponding ADQL query for this cross-match selection is provided in the Appendix B. It results in a preliminary dataset of 18,751,187 stars, or roughly about 1% of Gaia DR3. The left panel of Figure 1 shows the distribution of this sample in the color-magnitude diagram.  \nWe restrict this sample further by selecting stars along the main sequence, with the following color cut:  \n0 . 65 < B p -R p < 1 . 5 (1)  \nand a cut in the color-absolute magnitude plane:  \n3 . 83 × log( B p -R p )+5 . 05 < M G < 3 . 83 × log( B p -R p )+5 . 95 (2)  \nwith M G computed as  \nM G = G +5+5 × log 10 ( ϖ 1000 mas ) (3)  \nwhere ϖ is the measured parallax.  \nThis reduces our dataset to 5,979,478 stars. We additionally apply a dust cut, using the dustmaps library (Green 2018) with the SFD (Schlegel et al. 1998) dust reddening map. We discard sources with reddening values E ( B -V ) > 1.  \nWe then use the cross-match service provided by the NOIRLab 's Astro Data Lab to obtain the corresponding unWISE values for W 1 and W 2 bands. We remove sources matched with more than one counterpart in the unWISE cross-match table (where the original conesearch was of 1 . 5 arc-seconds), and we apply a distance cut between matched sources under 0 . 15 arc-seconds. We control for valid values in W1 and W2 magnitudes and  \nrestrict our sample to stars with unwise flags values of 0 for the bitmask for both W 1 and W 2 bands (i.e. sources are not deemed impacted by nearby bright sources or artifacts).  \nThis results in a final sample of 4,898,812 million stars. The right panel of Figure 1 shows the distribution of this sample in Gaia color-absolute magnitude diagram.", '3. METHOD': "We propose to approach the problem of identifying infrared-excess stars in a data-driven fashion. Namely, our goal is to fit (or learn) a model on the data to predict the expected values for W 1 and W 2 bands, and to use the error in prediction as an indicator for anomalies, and potential infrared-excesses. Specifically, we are going to look for confidently wrong predictions of the model.  \nWe use the following features as input to our model: photometric observations (magnitude) from Gaia DR3 and 2MASS ( G,B p , R p mean magnitudes, and J , H , K s ), colors computed from Gaia and 2MASS ( B p -R p , B p -G , G -R p , J -H , J -K s , H -K s , and G -J ), and absolute magnitude M G . We additionally use the ruwe , parallax and ebpminrp-gspphot features from Gaia .  \nWe fit independent Random Forest (RF) regressors to predict respectively the colors K s -W 1 and K s -W 2. We use the default setup from sklearn (Pedregosa et al. 2011a) using a mean-squared error criterion, 100 trees, and no maximum depth. Comparisons with other (more flexible) models showed no improvement in prediction quality on a held-out sample, and therefore we chose this simple yet effective method.  \nWe perform an 8-fold split, where different RFs are trained separately on each fold (12 . 5% of the data or about 600,000 stars). We thus get 7 'test' predictions (for which the star did not appear in the training sample) for each star in our sample, for each color. Figure 2 shows the feature importances for a Random Forest predicting K s -W 1 (light blue) and an RF predicting K s -W 2. Interestingly, while colors features have high importance for both (as expected, since we predict a color), predictions of W 2 seem to rely more on ebpminrp (reddening value) and parallax than the RF for W 1. The order of colors is also interesting: G -R p seem to have more importance than J -H (likely because other colors from 2MASS already have high importance).  \nFrom each of these 7 individual 'test' prediction for color K s -W 1 (resp. K s -W 2), we can compute back the predicted W 1 (resp. W 2) magnitude of a star i , using  \n- This work's candidates (Contardo et al, 2024)\n- Cotten et al, 2016 (Nearby IR-excess)\n- Cruz-Saenz de Miera et al, 2014 (Sun-like stars)\n- Moor et al, 2021 (Extreme Debris Disks)  \n<!-- image -->  \nFigure 1. Left panel shows our initial sample of 18M FGK-type stars in Gaia DR3 with good cross-matches in AllWISE and 2MASS , in the color-absolute-magnitude space. The absolute magnitudes were computed as shown in Eq. 3, and there were no corrections made for extinction or reddening by dust here. We plot IR-excess candidates from the literature that appear in this preliminary sample: from Cotten & Song (2016) (cyan, targeting nearby stars), Cruz-Saenz de Miera et al. (2014) (pink, targeting Sun-like stars), and Mo'or et al. (2021) (blue, targeting Extreme Debris Disks). Our final prime set of candidates is shown in red. The dashed black box illustrates our main-sequence cut. Right panel shows the final sample used for our search, after our main-sequence cut, dust cut and unWISE control cut as described in Section 2.  \n<!-- image -->  \nFigure 2. Feature Importance from two Random Forests fitted to predict W 1 (blue) and W 2 (red).  \n<!-- image -->  \nthe predicted color and the true observed K s magnitude of that star as:  \n˜ W 1 i,j = Ks i -RF 1 j ( x i ) (4)  \nWhere RF 1 j ( x i ) denotes the output (predicted color K s -W 1) of the Random Forest j for a star x i . We can similarly compute a predicted W 2 magnitude from the sets of Random Forest fitted to predict K s -W 2.  \nWe can also compute an estimated magnitude W 1 (resp. W 2) for a star i by combining the 7 test predictions, as the median prediction over the folds, which we denote:  \n̂ W i = Median( { ˜ W i,j } j ∈ F i ) (5)  \nWhere F i denotes the folds' indices where i is not in the training set. Using ̂ W i as our predictor, we obtain a Mean Absolute Error of 0 . 016 for W 1 and 0 . 023 for W 2.  \nThe left panel in Figure 3 shows our sample in the W -̂ W error-space. Stars with very high prediction errors can be regarded as anomalous, either in excess (negative values) or deficit (positive values), as it means that the observed magnitude significantly deviates from the model prediction. Interestingly, one can see that significant infrared deficit (positive values) in both bands (i.e. positive diagonal) are barely existent, as expected. We highlight as red points a set of preliminary candidates selected from a cut within that error space (dashed lines, corresponding to Equation 8 described below). We can also see in this Figure panel that there seem to be two 'types' of excess candidates: one kind lying on the 1to-1 line in terms of error in predicted magnitude, and another kind where the excess (error) is (roughly two times) stronger in W2, but less candidates lying in between those 'groups'.  \nWhile high prediction errors indicate potential anomalies, we also want to focus our search on highly confident incorrect predictions. Therefore, we propose to introduce additional criteria to flag candidates that (i) have a high prediction precision, i.e. the predictions across the folds have low variance (the error is not dependant on the training set), (ii) are in well predicted regions of the feature space (i.e. similar examples have a high accuracy), (iii) are in a well-populated region of the dataset (they are not outliers in the feature space).  \nPrecision metric -: We define the fold-MAD (Median Absolute Difference) metric for a data point i as the median of the absolute differences between each prediction from the RFs trained on different folds, and the global prediction ̂ W i :  \nfold-MAD i = Median( {| ̂ W i -˜ W i,j |} j ∈ F i ) (6)  \nThis metric can be seen as a precision metric, relying on the scatter in the predictions for the object i .  \nAccuracy of similar examples -: For each object in our sample, we take their K-nearest neighbors ( K = 30) in the following feature space, in each test-fold: H -K s , G -R,J -K s , G -J, B -G,M G , parallax, ebpminrp (reddening), and ruwe. (Those features have been identified as the 'most important' features for predictions by the Random Forests). The features are rescaled in a MinMax fashion before computing the neighbors. We refer to this KNN as NN-color. We compute the K-MAD of a star i in the test-set of fold j as the median of the absolute differences between the true and predicted magnitudes ˜ W k,j over the K neighbors in that test-set:  \nK-MAD i,j = Median( {| W k -˜ W k,j |} k ∈ NN-colour( i,j ) ) (7)  \nWhere NN( i, j ) denotes the K nearest neighbors indices for the object i in the test-set of fold j , and W denotes either W 1 or W 2 bands. In other words, this measures if similar examples to a given object are well predicted or not. We can use this to favor candidates that have high errors but a low K-MAD (i.e. candidates for which similar objects were correctly predicted), increasing our confidence in the anomalous nature of the object, or to reweight the significance of an error in prediction through the K-MAD (i.e. a small prediction error is 'more significant' if the neighbors are very accurately predicted). We note however that if we had a 'cluster' of, say, objects with an IR excess that were all very similar in the feature space (enough so that they'd form the  \nmajority of each others' neighborhood), this criterion could mechanically discard these objects. We consider however that, given the size of our sample, the feature space we look at and the rarity of the objects we are looking for, this is relatively unlikely to happen.  \nThe right panel of Figure 3 shows the distribution of our sample in the prediction error divided by the K-MAD (computed per fold), averaged across the test-folds, for W 1 and W 2. We also plot the standard deviation across the folds as error bars for a selection of infrared excess candidates (negative values) under this criterion, shown in orange. The red candidates (high prediction errors) from the left panel of the same Figure match the red points in the right panel, and vice-versa for orange points, showcasing how some high-error sources might be significantly down-weighted by the use of the K-MAD.  \nWell-populated region of the feature space -: Additionally, we can compute the distance between the source and its neighbors. We can use this metric as a criterion for flagging or removing candidates that might live in a lowdensity region of the input feature space (hence with possibly less reliable prediction, and for which the KMAD might not be indicative).  \nWe propose to further investigate a set of candidate stars that are selected through a combination of the prediction error criterion (red points in Figure 3), the error over K-MAD criterion (orange point in Figure 3, downweighting the points that are in a region of the featurespace that is badly predicted, favoring high-accuracy regions), and the fold-MAD, to ensure our candidates have small prediction uncertainties (low scatter). More specifically, we select the stars with an error in each band under the 0 . 008-th percentile value (dashed light-grey diagonal line in the left panel of Figure 3), and with both errors under -0 . 1 mag (horizontal and vertical dashed light-grey line):  \n( W 1 i -̂ W 1 i ) + ( W 2 i -̂ W 2 i ) < -0 . 42 ( W 1 i -̂ W 1 i ) < -0 . 1 ( W 2 i -̂ W 2 i ) < -0 . 1 (8)  \nFor the error over k-MAD cut, we similarly select stars with a sum of the mean of error over k-MAD under the 0 . 007-th percentile value, and with both means of error  \nFigure 3. Left panel shows the difference between observed W 1 (resp. W 2 on y -axis) and the predicted ̂ W 1 (resp. ̂ W 2). We select stars with a high error in both bands as a preliminary sample of IR-excess candidates, depicted in red (cuts as dashed grey lines corresponding to Eq. 8). The right panel shows the mean across 7-folds of the prediction error, divided by the Median Absolute Error of the object's 30 nearest neighbors in the respective fold. Error bars are standard deviations across the 7-folds. We select another set of candidates in this space, depicted in orange. Red data points match the ones in the left panel, and respectively the orange points match the ones in the right panel.  \n<!-- image -->  \n<!-- image -->  \nKMADj  \nover K-MAD under -5:  \n1 7 ∑ j ∈ F i ( W 1 i -˜ W 1 i,j ) K-MAD1 i,j + 1 7 ∑ j ∈ F i ( W 2 i -˜ W 2 i,j ) K-MAD2 i,j < -23 . 1 7 ∑ j ∈ F i ( W 1 i -˜ W 1 i,j ) K-MAD1 i,j < -5 1 7 ∑ j ∈ F i ( W 2 i -˜ W 2 i,j ) K-MAD2 i,j < -5 (9)  \n59  \nWhich corresponds to the dashed grey lines in the right panel of Figure 3.  \nWe additionally select candidates that have a fold-MAD under 0 . 01 mag:  \nMedian( {| ̂ W 1 i -˜ W 1 i,j |} j ∈ F i ) < 0 . 01 Median( {| ̂ W 2 i -˜ W 2 i,j |} j ∈ F i ) < 0 . 01 (10)  \nThis results in a preliminary set of excess candidates of 127 stars.  \nWe note that this threshold cut is rather arbitrary: we naturally focus on a few numbers of extreme outliers and enforce a relationship in the errors in both W 1 and W 2. Other strategies could be used to decide where to put this threshold. For instance, one could use the rate  \nor contamination in the symmetric, infrared-deficit side: such events should not occur, and thus this gives us some information about the 'scatter' in our data and predictions and possible false-detection rate in our candidates (see Section 4.7 and Appendix C for more details).", '3.1. Additional cuts to reduce false-detection by contamination and other artifacts': "On this preliminary set of 127 candidates, we evaluate the following additional features that could lead to false anomalies:  \n- · Crowding: For each object, we perform a conesearch in Gaia , 2MASS , and unWISE catalogues of 5 arcseconds. If another object is found within that cone for any of the catalogs, we flag it as potentially being contaminated, or in a crowded area.\n- · Figure of Merit (FoM) from Gaia : Following the recommendations in Dennihy et al. (2020) to alleviate source confusion when looking for IRexcess in WISE , we check the figure of merit ( score ) provided by ESA Gaia cross-match with AllWISE . This score assesses the quality or 'likeliness' of the match (the higher the better), however, it is not normalized. Dennihy et al. (2020) suggest as a rule of thumb to discard sources with  \n2  \nW  \n2  \nW  \n- a high SNR but a low FoM (they suggest that excesses with SNR > 10 but FoM < 4 are likely the result of source confusion). We do not have objects with a FoM under 4 in our preliminary set of candidates.\n- · Proper-Motions Comparison : Dennihy et al. (2020) also suggest to check proper motions agreements between catalogs to remove potential source confusion. They indicate that the proper motion accuracy in CatWise depends on the magnitude of the source (10 mas yr -1 for bright sources, 30 mas yr -1 for W 1 ≈ 15 . 5 mag). We compare the proper motions from Gaia and from CatWise and apply a cut when the proper motions in either RA or DEC disagree more than the approximated accuracy for the source brightness.\n- · Disagreement between unWISE and AllWISE : To ensure confidence in the observed magnitude of our targets, we add a flag checking for agreement between unWISE and AllWISE . unWISE documentation recommends subtracting 4 mmag from unWISE W1 and 32 mmag from unWISE W2 to improve agreement between both catalogues. After doing so, we add an 'agreementflag' such that | W AllWISE -W unWISE | < 0 . 15 for both bands.\n- · Distance to galactic plane : We remove candidates that are close to the galactic plane s.t. abs ( b ) < 10.\n- · Binary and duplicated sources : Candidates that are indicated as 'non single stars' by Gaia DR3 are flagged. We also check each candidate with Simbad (using Gaia DR3 ID) and flag each star that is indicated to be eclipsing binaries or eclipsing binaries candidates. Gaia DR3 also provides a duplicated source flag (which may indicate observational, cross-matching or processing problems, stellar multiplicity, and possible unreliable astrometry or photometry). We therefore discard candidates with this flag.  \nWe combine these different flags to obtain a curated final set of 53 MIR-excess candidates. We provide in Table 1 the number of candidates remaining after applying each cut. The Gaia DR3 and 2MASS ID of those 53 candidates are available in Table 2 in the Appendix A, and as a csv file at https://github.com/contardog/ NotATechnosignatureSearch along with the code for this pipeline and the following analysis.  \nFigure 4 shows where our candidates lie in Galactic coordinates l, b , where the dashed lines indicate the cut at abs ( b ) < 10. We provide in Figure 5 histograms of our candidates' sample and of our complete sample in B p -R p , [ M/H ] (obtained from Gaia DR3 mh gspphot ), Teff, and M G . We observe a gap in our candidates' distribution around 1 in B p -R p , 5.5 in M G and 5400K in Teff. There also seems to be an over-density in our candidates' distribution compared to the full sample of stars with color B p -R p < 0 . 9 and Teff ∼ 5900K. Because of the data-driven nature of our approach, and the subsequent cuts used to create our candidates' sample, as well as the small size of our final sample, it is hard to give clear interpretations and conclusions about these trends. We note however that the dip observed here at B p -R p ∼ 1 is not strongly visible when looking at the distribution of the candidates selected using the Prediction Error only (Eq. 8) or the k -MAD cut only (Eq. 9), and only start appearing when combining Eqs 8, 9 and 10 altogether, and is further depleted by subsequent cuts (at a somewhat similar rate than other color-bins). We also observe an offset trend in M/H towards lower values than the full sample, with a few outliers standing out.  \n<!-- image -->  \nl  \nFigure 4. Distribution of our full sample and our candidates in Galactic coordinates. Dashed black lines indicate abs ( b ) = 10.", '4.1. Characterization of the hosts': "For the four observables in Figure 5, we conducted a Kolmogorov-Smirnoff test which indicates that the differences in distribution with the full sample are significant (with p-values < 1 e -3 ).  \nTable 1. Table of criterion for our candidates' selection, and the number of candidates remaining after each cut.  \nFigure 5. Histograms of distribution of our candidates' sample (dark red outline) and our full 4.9M stars sample (blue) in B p -R p , effective temperature, [ M/H ], and M G .  \n<!-- image -->", '4.2. H-alpha emission': "Our sample focuses on the main sequence, but stars might have only recently joined and could be in a 'transitional' phase between protoplanetary disks and debris disks (Wyatt et al. 2015). To properly understand the nature(s) of the excesses observed in our candidates, and to identify the most interesting candidates for follow-up observations (e.g. confidently old stars with an unusually high IR excess), it is crucial to constrain the hosts' ages. We investigate different markers that could indicate that the host is young or recently joined the main sequence. However, the analysis presented below is only a first step and we defer thorough age-estimations of the candidates to future work.  \nTo investigate the potential 'youngness' of our candidates, we start by evaluating the H α emission. We first use the measurement of the pseudo-equivalent width of the H α line provided by Gaia DR3 from the RP spectra, ew espels halpha . The value is expected to be negative when emission is present. The left panel of Figure 6 shows the distribution of this measurement for our candidates' sample (dark red line), as well as a randomly selected control sample of 1,600 stars from our 4.9M sample (blue) and a set of stars similar to our candidates called 'Neighbors' (dashed black line). This set of neighbors is computed in a different features-space than previously. We use the following features to favor stars of similar magnitude and spectral shape (while our previous KNN relied on color): M G , ebpminrp gspphot, parallax, ruwe , G , J , H , K s , B p , and R p . We refer to this KNN as NN-mag.  \nWe can see two clear outliers in our set of candidates with high H α emission, indicating young stellar objects. The rest of our candidates seem to follow the control distribution, perhaps with a slight over-density on the absorption side. The pseudo-equivalent width EWH α of all candidates is given in the summary Table 2 in Appendix A.  \nAs a sanity check, we perform a different measurement of potential emissions or absorptions of H α as follows: using the XP spectra and GaiaXPy library, we measure the flux value of the spectra at the H α wavelength (656.46nm). For an object i , we obtain the flux values of its 30 nearest neighbors (using NN-mag), and compute the difference between the object's flux and the median flux of its neighbors, divided by the median flux:  \nE Hα,i = F Hα,i -Median k ∈ NN-mag( i ) ( F Hα,j ) Median k ∈ NN-mag( i ) ( F Hα,j ) (11)  \nThis metric provides an indicator of how much the object flux diverges at the H α line from similar objects, where a positive value would indicate emission (respectively absorption for negative values). The middle panel of Figure 6 shows the distribution of E Hα for our set of candidates, the control sample and the neighbors sample, and the right panel shows E Hα against the pseudowidth provided by Gaia data product. We see agreement especially in the detection of the two high emissions objects, Gaia DR3 5260430146906778624 (top right panel in Figure 15) and Gaia DR3 3004530769758177280 (third row, left panel of Figure 16). Note that the second object corresponds to our candidate with the highest fractional luminosity ( ∼ 0 . 1).  \nVisualizing the XP spectra also highlights some odd objects with respect to their NN-mag neighbors (in grey), see for instance the top middle candidate in Figure 14 ( Gaia DR3 6705842630930258816).  \nFor the sake of completeness and using the data product available with the catalogs we used, we also provide in Figure 7 a visualization of the age estimates provided by the Final Luminosity Age Mass Estimator (FLAME) within Gaia (Creevey et al. 2023) for our candidates and control samples. The 'error bars' correspond to the lower and upper confidence levels of the FLAME's age estimates (i.e. the 16th and 84th percentile value of the 1D projected distribution from sampling in mass and age). The x -axis is the 'evolution stage' from FLAME, where 100 corresponds to Zero Age Main Sequence (ZAMS), 360 to the main sequence turnoff, and 490 would be the base of the Red Giant Branch (Hidalgo et al. 2018). However, Creevey et al. (2023) highlights that those age estimates assume some prior on the metallicity, and the authors advise caution in using those ages for stars outside of the -0 . 5 < [ M/H ] < +0 . 5 regime. Therefore, we plot only in Figure 7 the stars that have a metallicity estimate (provided by Gaia ) in this range for the control and neighbors subsample. For our candidates, we plot in red the stars within the correct metallicity range, and in green the ones outside this range. All the youngest ( < 2Gyrs) candidates fall outside of the metallicity range. However, some of the remaining candidates have large error bars. One of the candidates with a high H α emission has an age estimate of 5.79 Gyrs (the other candidate does not have an estimate). We therefore recommend cautiousness in making any conclusion from those estimates.  \n<!-- image -->  \nFigure 6. Distribution of measurements for H α emission for a random control subsample (blue), a set of stars similar to our candidates (black dashed line) and our candidates set (dark red line). Left panel shows the distribution of the pseudo-equivalent width of the H α line provided by Gaia DR3 (negative values indicate emission). Middle panel shows our measurement using Equation 11 (positive values indicate emission). Right panel shows those measurements against each other.  \n<!-- image -->  \nFigure 7. Age estimates and evolution stage provided by the Final Luminosity Age Mass Estimator (FLAME) in Gaia DR3, for our candidates (dark red), a set of similar stars (grey) and a control sample (blue), restricted to the stars that have -0 . 5 < [ M/H ] < +0 . 5. The error bars depict the lower (16th percentile) and upper (84th percentile) ages to cover the 68% confidence interval. The x -axis Evolution Stage goes from 100 (Zero Age Main Sequence) to 490 (Red Giant Branch). 360 would be the main sequence turn-off.  \n<!-- image -->", '4.3. Check for Radio sources': "Following Ren et al. (2024), we search for cross-matching radio sources with our candidates: radio sources in the close vicinity of a candidate might indicate possible contamination with background galaxies (e.g. dust obscured). We search for radio sources within 15' of our candidates in catalogs from the Very Large Array Sky Survey (VLASS, Gordon et al. (2021)), Rapid ASKAP Continuum Survey (RACS, Hale et al. (2021)), the FIRST survey (Helfand et al. 2015), the NRAO VLA  \nSky Survey (NVSS, Condon et al. (1998)), the TIFR GMRT Sky Survey (TGSS, Intema et al. (2017)), and we add to this list the LOFAR Two-metre Sky Survey (LoTSS) DR2 (Shimwell et al. 2022).  \nWe found radio sources near 6 of our candidates. Four of those have a radio source within 1'; one within 2', and one at 10.'65. A flag indicating candidates with a nearby radio source is included in the summary Table 2 provided in the Appendix A.", '4.4. Variability of hosts': "Another matter of interest in sources of IR excess is the variability of the object. Following Mo'or et al. (2021), we use the NeoWISE time series to investigate first the variability of our candidates in the mid-infrared. We use the correlation-based Stetson index (Stetson 1996) S J , which measures the correlation of the variability in two (or more) bands and can be computed for a given object as follows:  \nS J = ∑ n j =1 w k sgn( P k ) √ | P k | ∑ n j =1 w k (12)  \nWhere n is the number of paired observations across the different bands, w k is a weight per epoch (uniform in our case) and :  \nP k = ( √ n n -1 W 1 ,k -¯ W 1 σ W 1 ,k )( √ n n -1 W 2 ,k -¯ W 2 σ W 2 ,k )  \nWhere σ W i,k is the uncertainty for the k -th observation in the i -th band, and ¯ W i is the mean of the observed points in band i .  \nWe compare the Stetson index obtained by our candidates against a random control sample and a set of  \nneighbors (computed with NN-color), shown on the left panel of Figure 8. We observe a small overdensity for higher Stetson indexes in our candidates' set, highlighting these candidates as potentially (more) variable objects. The object with the highest Stetson index also corresponds to an object flagged as a variable from AllWISE 's data var flag . However, the rest of our candidates have low IR variability scores according to this flag. To further get a sense of the variability of our candidates, we also measure their Stetson index percentile score (right panel Figure 8) with respect to the Stetson indexes of the control sample (in blue), and with respect to their 30 neighbors (dashed black line).  \nFigure 9 shows the NeoWISE time-series in both W 1 and W 2 for two of our candidates, one with a high percentile score (left, 'more variable') and one with a low percentile score (right).  \nWe did not identify any rising pattern as the one in Kenworthy et al. (2023) in our candidates' NeoWISE timeseries. Those types of objects might not get detected by our current pipeline, depending on how AllWISE and unWISE treat the aggregation of NeoWISE data, and when the rise would occur: we might get underestimated magnitudes from either catalogs, or pre-rise magnitudes (from AllWISE ). Additionally, we removed candidates with a disagreement between the two catalogs: as unWISE used a longer observation span from NeoWISE , a strong change in an object magnitude could lead to disagreement with AllWISE as well. A dedicated search for these specific 'rising IR emissions' would therefore be relevant, using either NeoWISE or unTimely (Meisner et al. 2023).  \nOn the optical side, we check the photometric variable flag provided by Gaia (computed in the G band). 10 out of our 53 candidates are flagged as variable, the remaining have not been processed (flag not available, hence they can not be considered not-variable either). This is a significantly larger occurrence of variable objects ( ∼ 19%) compared to the occurrence rate in our control sample ( ∼ 7%), the 'neighbors' sample ( ∼ 2 . 5%), and our entire sample ( ∼ 5 . 3%). Since variability is an indicator of potential youngness, this could hint that our sample tends to have younger stars than our control set. However, dimming and transit patterns have been observed in the optical for Sun-like stars with similar MIRexcess with ages of 300Myr (Kenworthy et al. 2023) and 600Myr (Melis et al. 2021). A closer investigation of the variability and periodicity of our candidates would thus be necessary.  \nAccording to the variable classification pipeline (Rimoldini et al. 2023) provided by Gaia , one of our candidates  \nis an eclipsing binary ( Gaia DR3 6705842630930258816). Two objects are (tentatively) classified as RS Canum Venaticorum variable ( Gaia DR3 3183786280737192704 and 2256421065353454464), one as a δ Scuti or γ Doradus or SX Phoenicis star. Four remaining are flagged as 'solar-like'.  \nUsing the cross-match of Gaia DR3 with known variable sources in the literature, provided by Gavras et al. (2023), we identify two candidates classified as Eclipsing, one that was previously identified as such above by the Gaia variable flag and Gaia DR3 5260430146906778624.  \nLonger light curves in the optical can be obtained for some of these objects (30) from TESS mission. However, a deeper analysis of e.g. periodicity and conducting gyrochronology will require a careful pre-processing and detrending of those data, which we defer to future work. Observations from the Zwicky Transient Facility (ZTF) could also be leveraged (17 candidates) through the SNAD viewer (Malanchev et al. 2023) which also provides cross-matches across multiple catalogs and additional variability detection from various pipelines.", '4.5. Characterization of the Excesses through SED-fitting': "We now further investigate our curated set of 53 candidates to derive properties of the possible source of excess.  \nFirst, we use Ariadne package (Vines & Jenkins 2022) to perform Spectral Energy Distribution (SED)-fitting on each of our candidates. Given RA,DEC positions and optionally a Gaia DR3 identifier, Ariadne automatically retrieves the object's photometric observations from Gaia DR2, AllWISE , APASS , Pan-STARRS1 , SDSS , 2MASS , ASCC , IRAC/GLIMPSE , TESS , SKYMAPPER , GALEX , and Tycho-2 . Ariadne also incorporates priors for interstellar extinction parameters through dust maps. We use the default setup relying on SFD dustmaps (Schlegel et al. 1998). The goal of Ariadne library is originally to improve SED-fitting through bayesian model averaging of different models. However, in the following of our analysis, we will use a single SEDmodel, the Castelli-Kurucz (Castelli & Kurucz 2003). We also use the 'fitzpatrick' extinction-law in Ariadne setup. Ariadne then returns the best-fitted SED as well as estimated stellar parameters (Teff, logg, distance, A V , metallicity and radius). Ariadne uses by default all observed photometry up to W 2 for the fit, but not W 3 and W 4.  \nFigure 8. Left panel shows the distribution of the Stetson index (Equation 12) of our control sample (blue), of a set of stars similar to our candidates (black dashed line) and of our candidates (dark red line). Right panel shows the distribution of the percentile score of our candidates' Stetson index with respect to the control sample (blue) and with respect to their set of similar stars (black dashed line).  \n<!-- image -->  \nWe follow this pipeline, using Gaia DR3 values as priors for Teff and logg . However, we force the removal of the WISE values. The motivation is two-fold: (i) we put the SED-modeling in a similar context as our model, (ii) since we will use the fitted SED to compute properties of the potential source of the excess, we do not want to bias the predicted values in the WISE bands.  \nWe obtain, for each of our candidates, SED-predicted values for all the WISE bands. Note that these values use potentially more observations at different wavelengths than our model, and rely on stellar model, extinction and other priors (contrary to our entirely data-driven predictions).  \nUsing the SED-predicted fluxes and the observed ones, we can estimate the properties of a blackbody fitted on the residuals. To do so, we build a grid of blackbody spectra with a range of temperatures from 100K  \nto 1500K with a step of 5K, and a range of angular sizes from 10 -23 to 10 -16 steradians with 200 bins spaced evenly on log-scale 2 . For each candidate, we select the blackbody that minimizes the sum of the squared errors as: argmin BB ∑ j ( F BB ( λ j ) -( F obs ( λ j ) -F SED ( λ j ))) 2 where j denotes an index on wavelengths. We fit the blackbody using four wavelengths corresponding to K s , W 1 , W 2, and W 3 bands. Because we optimize the sum of the squared errors in the flux units, the brighter bands are weighted more than the fainter ones in the fit. While we do not use W 3 for our detection pipeline because it is less reliable, we choose to use it for the blackbody fitting to extend our wavelengths baseline. But it is important to note that some of our candidates have a low SNR flag in W 3, which might lead to incorrect blackbody estimates.  \nWe show in Figure 10 the flux value observed in W 2 band versus the flux predicted by the SED-fit as red points, confirming excess for all of our candidates (lying under the 1-to-1 dashed line). The blue dots show the flux in W 2 of the SED with the additional blackbody, against the observed flux. Our fitting procedure tends to slightly over-correct (some points lie above the 1-to-1 dashed line) in the W 2 band, but not in W 1 and W 3. Figure 11 shows the SED from Ariadne (dashed red line), the black body fitted (dashed light grey line), and the combination of the two (black line), for one of our candidates.  \nFigure 12 shows the different black bodies (BB) properties we obtain (bearing in mind that those are approximate fits) with our candidates' sample. The fractional luminosities f BB = L BB L ∗ (computed as the ratio of the integrated flux of the sole black body and the integrated flux of the predicted SED) range between 0 . 005 and 0 . 1, with an average of 0 . 02. These fractional luminosities align with what has been observed and expected for Extreme Debris Disks candidates.  \nThe SED and black-body fits for all of our candidates are shown in Figures 14-17, along with their XP spectra (H α line indicated in blue) and imaging visualization of the objects in W 1 and W 2 bands using AllWISE , and in the optical using the Digitized Sky Survey (DSS), on 3 ' square. We observe in these Figures a diversity of 'excess shape' and strength: we see candidates with some extreme 'bumps' in W 3 - which are often not well fitted by our black body estimates -, similar to EDDs candidates. However, as we also display the photometric  \nGaia DR3 6727198140355736448  \n<!-- image -->", 'Gaia DR3 4420716452299884288': "Figure 9. NeoWISE time-series for two of our candidates, in W 1 (blue) and W 2 (red) bands. Left panel shows a candidate with a high percentile score with respect to the control sample. Right panel shows a candidate with a low percentile score.  \n<!-- image -->  \nFigure 10. Observed flux (in W 2 band) vs predicted flux from Ariadne -SED only (red points) and Ariadne -SED + a BlackBody (blue points), fitted on the flux residuals in bands K s , W 1 , W 2 and W 3. The 1-to-1 line is depicted as a dashed grey line.  \n<!-- image -->  \nquality flag of W 3 and W 4 from AllWISE 3 , we can see that for some of those objects, the W 3 band might not be reliable. This might also impact the relevance of our BB parameters fit. We also want to point out that a blackbody, or a single blackbody, might not be the best way to represent the source of the excess. Some works use two blackbodies e.g. Chen et al. (2014). We defer the investigation of more complex models to future work.  \nThese candidates might interestingly fall in a region of debris disks evolution and parameter space that was scarcely populated previously and provide an interesting continuum of disks' evolution: Wyatt et al. (2015) proposed a five steps classification in the evolution of disks (from protoplanetary to debris), based on a number of a factor including flux ratio at different wavelengths (see Fig. 1 in Wyatt et al. 2015). Our candidates would fall around HD 141569 for the flux ratio at 12 µm . However, the lack of reliable observations at higher wavelengths, dust mass estimates, and, as mentioned previously, reliable age estimates, prevent us from placing our candidates in this classification scheme accurately. Additionally, the work conducted in Wyatt et al. (2015) focused on A stars, while our candidates are FGK stars, which might make direct comparison difficult. This motivates nonetheless further observations on this sample to better connect and understand those potential objects as potentially 'transitional' or undergoing extreme and rare events.", '4.6. Recovery Rate of Injections': "As a sanity check, we perform a set of injections to evaluate the recovery rate of our pipeline depending on BB's parameters. We select as a sample of stars the set of neighbors (in color space using NN-color) of our candidate, ∼ 1500 objects. We create a grid of BB parameters, where the temperature ranges from 200 to 1575 K by step of 25 K , and the angular size ranges from 10 -23 to 10 -16 with 50 bins in log-scale. Each star in the sample is injected with all BBs, creating an artificially injected dataset of 4.3M objects. We only consider here the recovery as passing the cut of the first three masks used in this paper relying on prediction error (Eq 8),", 'Gaia DR3 1673109634053612416': "Figure 11. SED and BlackBody fitting for our excess candidates. Yellow triangles are the observations gathered and used by Ariadne to compute the SED (dashed dark red line). Red points are the unWISE (for W 1 and W 2) and AllWISE (for W 3 and W 4) values. Blue points are the prediction of our model (averaged over the 7-folds). The best-fit blackbody is depicted as a grey dashed line. The final SED computed as Ariadne SED + BlackBody is depicted as solid black line. The letters above W 3 and W 4 are AllWISE photometric quality flag.  \n<!-- image -->  \nthe confidence of error compared to similar objects (Eq 9) and the precision criterion (Eq 10). Figure 13 shows the resulting recovery rate for each temperature-angular size bin, and our candidates as the cyan points. Some of our candidates' estimated BB parameters lie in a low or no-recovery region: this is very likely explained by the fact that our BB's parameters estimation is done with four observations, and often a rather approximated fit (as can be seen in Figures 14-17).  \nIt is also important to note that our selection pipeline relies in the first cut on thresholds (Equations 8-10): relaxing these choices and threshold might lead to more detections and higher recovery rates of fainter MIR-excess, but also potentially to an increase in the false detection rate.", '4.7. Occurence rate': "Our injection analysis combined with our candidates does not however directly allow to make statements on the occurrence rate of MIR-excess events in FGK main sequence stars, due, among other things, to the threshold steps of our pipeline and our data-driven approach. However, analyzing the rate of detection of in-  \nfrared deficit can help interpret our detection rate in the excess direction. As we assume that MIR deficits are not physically possible, investigating if and when we reach the detection of MIR-deficit gives us an idea of when we reach the 'noise' level of our model, and the possible false-positive or contamination rate.  \nWith our threshold values, looking for infrared deficits of similar strength (that is, reversing the cuts and values) in the opposite direction as done with Eq. 8 and 9, leads to a preliminary set of 13 'deficit' candidates (compared to 127 in the excess direction). This could be interpreted, to some extent, as a possible false-positive rate in our preliminary excess sample of ∼ 10%.  \nWhen we apply the additional cuts of Table 1, all the 13 candidates are eliminated. This, besides confirming the relevance of these checks to remove false detection, also tends to indicate a relatively high purity in our excess final sample.  \nConversely, replicating the pipeline using upper percentiles to compute the thresholds for infrared-deficit detection, leads to a final set of 7 candidates (compared to our 53 candidates), 3 of which can easily be discarded from imaging as being contaminated, the others being  \nFigure 12. Black bodies (BB) and stellar properties: Top left panel shows the fractional luminosity against the angular size, colored by the BB's temperature. Top right panel shows the BB's temperature against the angular size, colored by the effective temperature of the star provided by Gaia . Bottom panels show the same space colored by the absolute magnitude M G (left) and G (right).  \n<!-- image -->  \nextremely small deviation (see Appendix C for the Figures).  \nThis combined with the injection analysis seems to indicate that our percentile approach for our threshold choice allows good coverage for recovery while not leading to too much false-positive contamination. Our (uncorrected) 'occurrence' rate (53 candidates out of ∼ 5 million stars, or ∼ 0 . 001%) is much lower than the ones suggested in e.g. Kennedy & Wyatt (2013) (1 in 10,000 for old ( > Gyr) systems -e.g. like BD+20 307-, and 1% for younger ( < 120Myr) systems). However, their estimates are for IR excess at 12 µ m, while our detection pipeline uses observations at ∼ 2 -3 µ m.", '5. DISCUSSION': "We proposed and executed a data-driven approach for the discovery of infrared (IR) excesses around stellar objects, more specifically among main-sequence FGK stars. We used robust machine learning methods to bypass stellar modeling to estimate the 'expected' midIR (MIR) magnitudes, and we designed a set of curated statistics derived from the predictions and the predictions' quality to define a set of criterion cuts. We combined those with additional cuts to prevent as much as  \npossible false detection and artificial excess due, e.g., to contamination. We obtained a set of 53 candidates that display interesting (M)IR excess. A few of those objects appear to be young stars (showing high Hα emission). A significant fraction of our candidates seem to have variability in the optical, and some in the mid-infrared. This can also indicate youth, but a proper age estimation using gyrochronology and a variability analysis is required. Using SED-fitting and fitting black bodies on the flux residuals in 4 bands ( K s , W 1 , W 2 , W 3), we compute black body properties and obtain estimated fractional luminosities coherent with the ones associated with Extreme Debris Disks candidates and planetary collision candidates. These objects form an interesting set of candidates to investigate further, as they show a variety in their excess' morphologies, some of which might require more refined models than a simple single blackbody.  \nIt is natural to follow this work by gathering follow-up observations of the candidates, to confirm their nature (or discard potential confounding factors), to obtain better estimates of the properties of those disks, their compositions and their underlying processes. For example, Mo'or et al. (2024) recently investigated the presence of solid-state features around 10 µm in eight EDDs, and divided them into silica or silicate-dominated groups.  \nFigure 13. Map of the recovery rate of Black Body injections in a sample of stars similar to our candidates ( ∼ 1,500 stars). Our candidates' BB estimates are shown as cyan dots.  \n<!-- image -->  \nAnother key aspect will be to perform a deeper analysis to obtain reliable age estimates: the most interesting objects (as in: the ones that would challenge current models the most or would be the least expected) would be non-young. Additionally, identifying potential trends (if any) between those excesses and stellar ages might help better understand and constrain the underlying physical processes leading to such excesses.  \nAs we mentioned in the introduction, technosignatures such as Dyson spheres or swarms would lead to infrared excesses, of potentially different morphologies depending on the level of completion of the sphere/swarm. Suazo et al. (2024) recently conducted a search where they looked for specific types of IR-excess matching the Dyson spheres' models they defined. It would be interesting to follow a similar analysis on our set of candidates -which was found through a more flexible search, and thus all might not necessarily match the IR-excess morphologies of DS- to investigate if some of our candidates are potential matches and at which level of completion.  \nAnother natural extension of this work is to investigate other stellar types. These kinds of IR-excess have also been observed around A-stars (e.g. in Melis et al. (2013)), or in M-dwarfs (e.g. in the techno signature search of Suazo et al. (2024)). Replicating this work for  \ndifferent stellar types separately could lead to interesting observations, in terms of event rates and, further down the line, disk types.  \nWhile such a pipeline is relevant to identify a set of outliers candidates, it can only find (contextual) outliers according to the data . Making population statements (e.g. 'true' occurrence rate of planetary collisions) is not trivially doable with this approach, as it relies heavily on a series of threshold cuts to build the candidate sample. Additionally, since our model computing the 'expected (M)IR magnitude' is data-driven, we cannot directly use the error in prediction as the true IR excess (or deficit) for the full sample to investigate for instance the fraction of different levels of IR excess. We also note that, in principle, the prediction errors can be due mainly to variance in the data (and hence the model not being perfectly accurate) and might not necessarily indicate 'interesting anomalousness'. Here we can benefit from the fact that we are looking for the largest deviations, we expect (and have) a small sample of outliers to investigate, and we have other methods to validate an excess after detection. We also benefit from the fact that in this application, symmetric opposite anomalies (in the deficit side) are not physically possible. Thus we can use those as a probe for false-positive detection and/or reaching  \n'noise' level with our model. This would not always be available in other applications or general cases.", 'ACKNOWLEDGMENTS': "It is a pleasure to thank Jason Wright, Phil Armitage, and Erwan Allys for valuable discussions, as well as our anonymous reviewer and Carl Melis for their useful suggestions.  \nThis work has made use of data from the European Space Agency (ESA) mission Gaia (https://www.cosmos.esa.int/gaia), processed by the Gaia Data Processing and Analysis Consortium (DPAC, https://www.cosmos.esa.int/web/gaia/dpac/ consortium). Funding for the DPAC has been provided by national institutions, in particular the institutions participating in the Gaia Multilateral Agreement.  \nThis publication makes use of data products from the Two Micron All Sky Survey (doi: 10.26131/IRSA2 (Skrutskie et al. 2003)), which is a joint project of the University of Massachusetts and the Infrared Processing and Analysis Center/California Institute of Technology, funded by the National Aeronautics and Space Administration and the National Science Foundation.  \nThis publication makes use of data products from the Wide-field Infrared Survey Explorer, which is a joint project of the University of California, Los Angeles, and the Jet Propulsion Laboratory/California Institute of Technology, and NEOWISE, which is a project of the Jet Propulsion Laboratory/California Institute of Technology. WISE and NEOWISE are funded by the National Aeronautics and Space Administration.  \nThis publication also makes use of unWISE data product (doi:10.26131/IRSA525 (unWISE team 2019; Schlafly et al. 2019)) and AllWISE data products (doi: 10.26131/IRSA1 (Wright et al. 2019)).  \nThis research uses services or data provided by the Astro Data Lab at NSF's National Optical-Infrared Astronomy Research Laboratory. NOIRLab is operated by the Association of Universities for Research in Astronomy (AURA), Inc. under a cooperative agreement with the National Science Foundation.  \nThe Flatiron Institute is a division of the Simons Foundation.  \nSoftware: Astropy (Astropy Collaboration et al. 2013, 2018, 2022), numpy (Harris et al. 2020), ARIADNE (Vines & Jenkins 2022), scikit-learn (Pedregosa et al. 2011b), astroquery (Ginsburg et al. 2019), matplotlib  \n(Hunter 2007), pyvo (Graham et al. 2014), scipy (Virtanen et al. 2020), pandas (pandas development team 2020), dustmaps (Green 2018), gaiaXPy (GaiaXPy colab 2023).  \nFigure 14. SED and BlackBody fitting for our excess candidates. Yellow triangles are the observations gathered and used by Ariadne to compute the SED (dashed dark red line). Red points are the unWISE (for W 1 and W 2) and AllWISE (for W 3 and W 4) values. Blue points are the prediction of our model. The best-fit blackbody is depicted as grey dashed line. The final SED computed as Ariadne SED + BlackBody is depicted as a solid black line. Each plot also displays imaging from AllWISE ( W 1 , W 2) and DSS of a 3 ' square centred on the source candidate. The bottom part of each plot is the XP absolute spectra from Gaia DR3 computed using Gaiaxpy of the object (in red). The spectras of the 30 nearest-neighbors (using NN-mag) are in grey. H α line is shown in blue. The letters above W 3 and W 4 are AllWISE photometric quality flag.  \n<!-- image -->  \nFigure 15. SED and BlackBody fitting for our excess candidates. Yellow triangles are the observations gathered and used by Ariadne to compute the SED (dashed dark red line). Red points are the unWISE (for W 1 and W 2) and AllWISE (for W 3 and W 4) values. Blue points are the prediction of our model. The best-fit blackbody is depicted as grey dashed line. The final SED computed as Ariadne SED + BlackBody is depicted as a solid black line. Each plot also displays imaging from AllWISE ( W 1 , W 2) and DSS of a 3 ' square centred on the source candidate. The bottom part of each plot is the XP absolute spectra from Gaia DR3 computed using Gaiaxpy of the object (in red). The spectras of the 30 nearest-neighbors (using NN-mag) are in grey. H α line is shown in blue. The letters above W 3 and W 4 are AllWISE photometric quality flag.  \n<!-- image -->  \nFigure 16. SED and BlackBody fitting for our excess candidates. Yellow triangles are the observations gathered and used by Ariadne to compute the SED (dashed dark red line). Red points are the unWISE (for W 1 and W 2) and AllWISE (for W 3 and W 4) values. Blue points are the prediction of our model. The best-fit blackbody is depicted as grey dashed line. The final SED computed as Ariadne SED + BlackBody is depicted as a solid black line. Each plot also displays imaging from AllWISE ( W 1 , W 2) and DSS of a 3 ' square centred on the source candidate. The bottom part of each plot is the XP absolute spectra from Gaia DR3 computed using Gaiaxpy of the object (in red). The spectras of the 30 nearest-neighbors (using NN-mag) are in grey. H α line is shown in blue. The letters above W 3 and W 4 are AllWISE photometric quality flag.  \n<!-- image -->  \nFigure 17. SED and BlackBody fitting for our excess candidates. Yellow triangles are the observations gathered and used by Ariadne to compute the SED (dashed dark red line). Red points are the unWISE (for W 1 and W 2) and AllWISE (for W 3 and W 4) values. Blue points are the prediction of our model. The best-fit blackbody is depicted as a grey dashed line. The final SED computed as Ariadne SED + BlackBody is depicted as a solid black line. Each plot also displays imaging from AllWISE W 1 , W 2 and DSS of a 3 ' square centered on the source candidate. The bottom part of each plot is the XP absolute spectra from Gaia DR3 computed using Gaiaxpy of the object (in red). The spectra of the 30 nearest neighbors (using NN-mag) are in grey. H α line is shown in blue. The letters above W 3 and W 4 are AllWISE photometric quality flag.  \n<!-- image -->  \nWavelength [nm]", 'A. CANDIDATES TABLE': "Table 2 provides a summary view of our candidates with Gaia DR3 source ID, magnitudes in G and W 2 bands, variable flags from AllWISE and Gaia (when available, NaN does not mean a candidate is not variable), Stetson index S J , effective temperatures of our blackbody fit, fractional luminosities according to blackbody fit, pseudo-equivalent width of the H α line from Gaia (negative value corresponds to emission), Gaia 's age estimate, and the distance of the closest radio source within 15' if any.  \nTable 2 . Summary table of Candidates  \nTable 2 continued  \nTable 2 (continued)", 'B. ADQL QUERY FOR GAIA DR3 DATA SELECTION AND CROSS-MATCH': "Gaia ADQL query for data selection and cross-match with AllWISE and 2MASS of FGK stars. This query has to be looped over ranges of start index and end index using Gaia DR3 random index provided, to cover the entire Gaia data.  \n```\nSELECT dr3.*, xmatch.*, allw.*, tmassmatch.*, xjoin.*, tmass.*, dr3\\_all.duplicated\\_source FROM gaiadr3.gaia\\_source\\_lite AS dr3 JOIN gaiadr3.gaia\\_source AS dr3\\_all USING (source\\_id) JOIN gaiadr3.allwise\\_best\\_neighbour AS xmatch USING (source\\_id) JOIN gaiadr1.allwise\\_original\\_valid AS allw ON xmatch.allwise\\_oid = allw.allwise\\_oid JOIN gaiadr3.tmass\\_psc\\_xsc\\_best\\_neighbour AS tmassmatch USING (source\\_id) JOIN gaiadr3.tmass\\_psc\\_xsc\\_join AS xjoin USING (clean\\_tmass\\_psc\\_xsc\\_oid)\n```  \n```\nJOIN gaiadr1.tmass\\_original\\_valid AS tmass ON xjoin.original\\_psc\\_source\\_id = tmass.designation WHERE (dr3.random\\_index BETWEEN start\\_index AND end\\_index AND dr3.parallax\\_over\\_error > 10.0 AND dr3.teff\\_gspphot BETWEEN 4000 AND 7000 AND dr3.phot\\_g\\_mean\\_mag < 16 AND dr3.ruwe < 1.4 AND dr3.ebpminrp\\_gspphot IS NOT NULL AND allw.cc\\_flags = '0000' AND tmass.ks\\_msigcom < .2 AND tmass.h\\_msigcom < .2 AND tmass.j\\_msigcom < .2 AND tmass.ph\\_qual = 'AAA' AND tmassmatch.number\\_of\\_neighbours = 1 AND xmatch.number\\_of\\_mates = 0 AND xmatch.number\\_of\\_neighbours = 1 AND tmassmatch.number\\_of\\_mates = 0 AND xmatch.angular\\_distance < .15 AND tmassmatch.angular\\_distance < .15)\n```", 'C. INFRARED DEFICIT': "As noted in Section 4.7, using the same values used for Eq. 8 and 9 in a symmetric way to look for infrared deficit , i.e. looking for deficit of similar 'strength', leads to only 13 deficit candidates in the preliminary cut (as opposed to 127 in the excess direction). All of them then get discarded either by the secondary cuts (Table 1).  \nThus, we have no IR-deficit candidates that are symmetrically similar in strength to our IR-excess candidates.  \nApplying a similar methodology where Eqs 8 and 9 using upper percentiles, i.e. looking at the 'most deficit' objects (but not necessarily as strong in deficit as our excess candidates), we obtain 83 preliminary candidates ( vs 127 for excess). Applying all the secondary cuts leaves only 7 candidates (shown Figure 18). This highlights the efficiency of those secondary cuts to potentially remove false candidates.  \nThree of those have likely contamination from their environment shown by the imaging. They appear as a relatively small deficit in W 1 and W 2 according to the SED-fit too. For all of these examples, the W 3 and W 4 land above the SED-fit. However, they also all have a quality flag U (meaning that the source measurement had an SNR < 2, and the profile-fit magnitude is a 95% confidence upper limit).", 'D. OVERLAP WITH PREVIOUSLY DETECTED IR-EXCESS': "We check our pipeline selection and list of candidates against existing IR-excess candidates from the literature. We summarize this search in Table 3. For each reference investigated, we retrieve catalogs from either SIMBAD or Vizier, retrieve the Gaia DR3 match from SIMBAD if there is one (second column), and check if those objects are in our full 4.9M sample (third column). We indicate which of those objects make it through our preliminary cut (combining Equations 8,9 and 10) (fourth column) and through the complete pipeline (fifth column). Note that some objects in those catalogs might not be matched yet with a DR3 identifier and could thus actually appear in our sample.  \nThe excesses most similar to our candidates are likely to be Extreme Debris Disks candidates, and we indeed find that these objects are the ones that most often make it through our preliminary cut. Out of 23 objects matched from the catalog available on SIMBAD 4 , 5 appear in our full sample, and they all make it through our preliminary cut, confirming the ability of our pipeline to detect such IR-excesses. Two of those are then discarded by the crowding cut, one is discarded by the proper-motion disagreement cut, and one by the k-NN 'density' cut. One remains in our final set of candidates, Gaia DR3 2942418533272205312 (bottom left plot in Figure 17).  \nOut of the few other IR-excess candidates appearing in our 4.9M sample (42 in total), only one passes through the preliminary cut (from the Reserved catalog by Cotten & Song (2016)), but is discarded by the crowding criterion.  \nThis encourages us to believe that we explored a large sample of stars relatively disjointed from previous searches, both in the population sample analyzed (little overlap in the main catalog) and in the type of excess detected. It is relevant to note that the IR-excesses in the literature were detected in W 3 and W 4 bands, and hence might not produce an  \nFigure 18. SED, imaging and XP spectra of IR-deficit candidates .  \n<!-- image -->  \nWavelength [nm]  \nexcess in the MIR detectable by our pipeline. It also requires stringent cuts on the quality and SNR of W 4, which eliminates a large number of objects. Our detection pipeline allows for the discovery of extreme outliers (similar to potential EDDs) by relying on W 1 and W 2 bands only, while being rather conservative in further removing potential false detection.", 'REFERENCES': 'Arun, R., Mathew, B., Manoj, P., et al. 2019, AJ, 157, 159, doi: 10.3847/1538-3881/ab0ca1 Astropy Collaboration, Robitaille, T. P., Tollerud, E. J., et al. 2013, A&A, 558, A33, doi: 10.1051/0004-6361/201322068  \nAstropy Collaboration, Price-Whelan, A. M., Sip"ocz, B. M., et al. 2018, AJ, 156, 123, doi: 10.3847/1538-3881/aabc4f  \nAstropy Collaboration, Price-Whelan, A. M., Lim, P. L., et al. 2022, ApJ, 935, 167, doi: 10.3847/1538-4357/ac7c74 Babusiaux, C., Fabricius, C., Khanna, S., et al. 2023, A&A,  \n674, A32, doi: 10.1051/0004-6361/202243790 Ballering, N. P., Rieke, G. H., Su, K. Y. L., & Montiel, E. 2013, ApJ, 775, 55, doi: 10.1088/0004-637X/775/1/55  \nTable 3. Table of overlap between existing IR-excess catalogs and our sample and set of candidates. We retrieve the DR3 identifiers counterparts of the candidates using SIMBAD and indicate the number of objects in our full sample, the number of objects passing the preliminary cut and the full pipeline selection.  \nBalog, Z., Kiss, L. L., Vink\'o, J., et al. 2009, ApJ, 698, 1989, doi: 10.1088/0004-637X/698/2/1989  \nCastelli, F., & Kurucz, R. L. 2003, in Modelling of Stellar Atmospheres, ed. N. Piskunov, W. W. Weiss, & D. F. Gray, Vol. 210, A20.  \nhttps://arxiv.org/abs/astro-ph/0405087  \n- Chandola, V., Banerjee, A., & Kumar, V. 2009, ACM computing surveys (CSUR), 41, 1\n- Chen, C. H., Mittal, T., Kuchner, M., et al. 2014, ApJS, 211, 25, doi: 10.1088/0067-0049/211/2/25  \nCondon, J. J., Cotton, W. D., Greisen, E. W., et al. 1998, AJ, 115, 1693, doi: 10.1086/300337 Cotten, T. H., & Song, I. 2016, ApJS, 225, 15, doi: 10.3847/0067-0049/225/1/15 Creevey, O. L., Sordo, R., Pailler, F., et al. 2023, A&A, 674, A26, doi: 10.1051/0004-6361/202243688 Cruz-Saenz de Miera, F., Chavez, M., Bertone, E., & Vega, O. 2014, MNRAS, 437, 391, doi: 10.1093/mnras/stt1888 Cutri, R. M., Wright, E. L., Conrow, T., et al. 2021, VizieR Online Data Catalog, II/328 Da Costa, A. D., Canto Martins, B. L., Le˜ao, I. C., et al. 2017, ApJ, 837, 15, doi: 10.3847/1538-4357/837/1/15 Dennihy, E., Farihi, J., Gentile Fusillo, N. P., & Debes, J. H. 2020, ApJ, 891, 97, doi: 10.3847/1538-4357/ab7249 Gaia Collaboration, Prusti, T., de Bruijne, J. H. J., et al. 2016, A&A, 595, A1, doi: 10.1051/0004-6361/201629272 Gaia Collaboration, Vallenari, A., Brown, A. G. A., et al. 2023, A&A, 674, A1, doi: 10.1051/0004-6361/202243940 GaiaXPy colab. 2023, GaiaXPy, doi: 10.5281/zenodo.8239995  \nGavras, P., Rimoldini, L., Nienartowicz, K., et al. 2023, A&A, 674, A22, doi: 10.1051/0004-6361/202244367 Ginsburg, A., Sip"ocz, B. M., Brasseur, C. E., et al. 2019, AJ, 157, 98, doi: 10.3847/1538-3881/aafc33 Gordon, Y. A., Boyce, M. M., O\'Dea, C. P., et al. 2021, ApJS, 255, 30, doi: 10.3847/1538-4365/ac05c0 Graham, M., Plante, R., Tody, D., & Fitzpatrick, M. 2014, PyVO: Python access to the Virtual Observatory, Astrophysics Source Code Library, record ascl:1402.004 Green, G. 2018, The Journal of Open Source Software, 3, 695, doi: 10.21105/joss.00695 Hale, C. L., McConnell, D., Thomson, A. J. M., et al. 2021, PASA, 38, e058, doi: 10.1017/pasa.2021.47 Harris, C. R., Millman, K. J., van der Walt, S. J., et al. 2020, Nature, 585, 357, doi: 10.1038/s41586-020-2649-2 Helfand, D. J., White, R. L., & Becker, R. H. 2015, ApJ, 801, 26, doi: 10.1088/0004-637X/801/1/26 Hidalgo, S. L., Pietrinferni, A., Cassisi, S., et al. 2018, ApJ, 856, 125, doi: 10.3847/1538-4357/aab158 Hunter, J. D. 2007, Computing in Science & Engineering, 9, 90, doi: 10.1109/MCSE.2007.55 Intema, H. T., Jagannathan, P., Mooley, K. P., & Frail, D. A. 2017, A&A, 598, A78, doi: 10.1051/0004-6361/201628536 Ishida, E. E., Kornilov, M. V., Malanchev, K. L., et al. 2021, Astronomy & Astrophysics, 650, A195 Kennedy, G. M., & Wyatt, M. C. 2013, MNRAS, 433, 2334, doi: 10.1093/mnras/stt900 Kenworthy, M., Lock, S., Kennedy, G., et al. 2023, Nature, 622, 251, doi: 10.1038/s41586-023-06573-9  \nLochner, M., & Bassett, B. A. 2021, Astronomy and  \nComputing, 36, 100481 Malanchev, K., Kornilov, M. V., Pruzhinskaya, M. V., et al. 2023, PASP, 135, 024503, doi: 10.1088/1538-3873/acb292 Manzo-Mart\'ınez, E., Calvet, N., Hern\'andez, J., et al. 2020, ApJ, 893, 56, doi: 10.3847/1538-4357/ab7ead McDonald, I., Zijlstra, A. A., & Watson, R. A. 2017, MNRAS, 471, 770, doi: 10.1093/mnras/stx1433 Meisner, A. M., Caselden, D., Schlafly, E. F., & Kiwy, F. 2023, AJ, 165, 36, doi: 10.3847/1538-3881/aca2ab Melis, C., Olofsson, J., Song, I., et al. 2021, ApJ, 923, 90, doi: 10.3847/1538-4357/ac2603 Melis, C., Zuckerman, B., Rhee, J. H., & Song, I. 2010, ApJL, 717, L57, doi: 10.1088/2041-8205/717/1/L57 Melis, C., Zuckerman, B., Rhee, J. H., et al. 2012, Nature, 487, 74, doi: 10.1038/nature11210 -. 2013, ApJ, 778, 12, doi: 10.1088/0004-637X/778/1/12 Mo\'or, A., \' Abrah\'am, P., Szab\'o, G., et al. 2021, ApJ, 910, 27, doi: 10.3847/1538-4357/abdc26 Mo\'or, A., \' Abrah\'am, P., Su, K. Y. L., et al. 2024, MNRAS, 528, 4528, doi: 10.1093/mnras/stae155 pandas development team, T. 2020, pandas-dev/pandas: Pandas, latest, Zenodo, doi: 10.5281/zenodo.8092754 Pedregosa, F., Varoquaux, G., Gramfort, A., et al. 2011a, Journal of Machine Learning Research, 12, 2825 -. 2011b, Journal of Machine Learning Research, 12, 2825 Ren, T., Garrett, M. A., & Siemion, A. P. V. 2024, Research Notes of the American Astronomical Society, 8, 145, doi: 10.3847/2515-5172/ad5017 Rhee, J. H., Song, I., & Zuckerman, B. 2008, ApJ, 675, 777, doi: 10.1086/524935 Rimoldini, L., Holl, B., Gavras, P., et al. 2023, A&A, 674, A14, doi: 10.1051/0004-6361/202245591 Schlafly, E. F., Meisner, A. M., & Green, G. M. 2019, ApJS, 240, 30, doi: 10.3847/1538-4365/aafbea Schlegel, D. J., Finkbeiner, D. P., & Davis, M. 1998, ApJ, 500, 525, doi: 10.1086/305772  \nSchmidt, E. G. 2015, ApJ, 813, 29, doi: 10.1088/0004-637X/813/1/29 Sheikh, S. Z. 2020, International Journal of Astrobiology, 19, 237, doi: 10.1017/S1473550419000284 Shimwell, T. W., Hardcastle, M. J., Tasse, C., et al. 2022, A&A, 659, A1, doi: 10.1051/0004-6361/202142484 Skrutskie, M. F., Cutri, R. M., Stiening, R., et al. 2003, 2MASS Source Catalog, IRSA, doi: 10.26131/IRSA2 Skrutskie, M. F., Cutri, R. M., Stiening, R., et al. 2006, AJ, 131, 1163, doi: 10.1086/498708 Song, I., Zuckerman, B., Weinberger, A. J., & Becklin, E. E. 2005, Nature, 436, 363, doi: 10.1038/nature03853 Stetson, P. B. 1996, PASP, 108, 851, doi: 10.1086/133808 Suazo, M., Zackrisson, E., Mahto, P. K., et al. 2024, MNRAS, 531, 695, doi: 10.1093/mnras/stae1186 Suazo, M., Zackrisson, E., Wright, J. T., Korn, A. J., & Huston, M. 2022, Monthly Notices of the Royal Astronomical Society, 512, 2988, doi: 10.1093/mnras/stac280 unWISE team. 2019, unWISE catalog, IRSA, doi: 10.26131/IRSA525 Uzpen, B., Kobulnicky, H. A., Olsen, K. A. G., et al. 2005, ApJ, 629, 512, doi: 10.1086/431479 Vines, J. I., & Jenkins, J. S. 2022, MNRAS, doi: 10.1093/mnras/stac956 Virtanen, P., Gommers, R., Oliphant, T. E., et al. 2020, Nature Methods, 17, 261, doi: 10.1038/s41592-019-0686-2 Wright, E., Eisenhardt, P., Mainzer, A., Ressler, M., & Cutri, R. e. a. 2019, AllWISE Source Catalog, IRSA1, doi: 10.26131/IRSA1 Wyatt, M. C., Pani\'c, O., Kennedy, G. M., & Matr\'a, L. 2015, Ap&SS, 357, 103, doi: 10.1007/s10509-015-2315-6 Zuckerman, B., Melis, C., Rhee, J. H., Schneider, A., & Song, I. 2012, ApJ, 752, 58, doi: 10.1088/0004-637X/752/1/58'}