© The Authors 2026

1. Introduction

Carbon stars are characterized by an atmospheric carbon-to-oxygen number ratio C/O higher than 1. The C/O can result either from extrinsic processes such as the accretion of C-enriched material or from physical processes originating in the star itself. The latter is the case for thermally pulsing asymptotic giant branch (TP-AGB) stars. During their evolution along the AGB phase, stars with a mass M ≳ 1 M_⊙ (Di Criscienzo et al. 2016) experience the third dredge-up. This event mixes the products of nuclear burning into the outer envelope, enriching the photosphere with C. However, stars with masses M ≳ 3.5 M_⊙ are prevented from becoming carbon stars, due to the bottom burning (Iben & Renzini 1983) process that takes place at the base of the convective envelope of these stars and activates the CN cycle that converts carbon into nitrogen. Stars with masses between 1 and 3.5 solar masses (Di Criscienzo et al. 2016) undergo recurrent thermal pulses and dredge-up episodes. These mass limits strictly depend on the initial chemical composition and vary according to the adopted stellar evolution models (Ventura et al. 2016; Dell’Agli et al. 2017; Pastorelli et al. 2019; Abia et al. 2020; Marigo et al. 2022; Straniero et al. 2023). A series of third dredge-up events mix the products of nuclear burning into the outer envelope, enriching the photosphere with C. At the relatively low temperatures of AGB atmospheres, molecules form efficiently. Due to the high stability of the carbon monoxide molecule (CO), which consumes the less abundant of the two elements, the C/O abundance ratio becomes the key parameter: C/O > 1 leads to carbon-rich (C-rich; Iben 1983) stars, while C/O < 1 results in oxygen-rich (O-rich) stars. The carbon enrichment increases the opacity of the outer layers, and together with pulsation-driven shocks creates favourable conditions for dust condensation and growth. Radiation pressure on dust grains then drives an intense stellar wind, ejecting gas and dust into the interstellar medium (Höfner & Olofsson 2018).

As a result, carbon-rich AGB stars are one of the most important contributors to the enrichment of the interstellar medium (e.g. Marigo 2002; Karakas & Lattanzio 2007). Their dust production is dominated by amorphous carbon (amC), with a minor contribution from silicon carbide (SiC), as confirmed by several studies (e.g. Groenewegen 2022). Thus, carbon-rich AGB stars play a significant role in the chemical evolution of galaxies and in shaping their observed properties. They account for a substantial fraction of the infrared (IR) luminosity and are key contributors to the integrated light of unresolved galaxies (e.g. Maraston et al. 2006; Marigo et al. 2009; Boyer et al. 2011; Melbourne & Boyer 2013). In this work we focus on intrinsic carbon-rich AGB stars (resulting from the third dredge-up) and, for the sake of brevity, we hereafter refer to them simply as carbon stars (or C stars).

The properties of carbon stars make them promising distance indicators. In particular, the J-region asymptotic giant branch (JAGB) method (Freedman & Madore 2020) uses a subset of these stars as standard candles. This approach has been shown to provide distance estimates consistent with those from the tip of the red giant branch (TRGB) method (Lee et al. 1993), with recent applications confirming its accuracy (Lee 2023). Moreover, thanks to their higher luminosity, JAGB stars can be detected at greater distances than TRGB stars (Lee et al. 2024), offering a potential extension of the extragalactic distance scale and contributing to the calibration of the Hubble constant.

The advent of large-scale surveys such as Gaia (Gaia Collaboration 2016), the Two Micron All Sky Survey (2MASS; Skrutskie et al. 2006), and the All-Sky Release Catalog (AllWISE; Cutri et al. 2021) of the Wide-field Infrared Survey Explorer (WISE; Wright et al. 2010), has enabled the identification and characterization of a large number of carbon stars across the Milky Way (MW) and in the Magellanic Clouds (MCs). Taken together, these surveys provide extensive and homogeneous photometric coverage, from 0.5 to 22 microns, allowing a more reliable determination of the properties of carbon stars across different stellar populations and galactic environments.

In the last decade, numerous studies have been conducted to determine the properties of carbon stars by using the spectral energy distribution (SED) fitting method through multi-band photometry data (e.g. Bressan et al. 1998; Groenewegen et al. 2009; Riebel et al. 2012; Boyer et al. 2013; Srinivasan et al. 2016; Rau et al. 2017; Nanni et al. 2018; Groenewegen & Sloan 2018; Nanni et al. 2019; Nanni 2019; Groenewegen 2022). By comparing the observed photometric data to the available synthetic models, these studies have been able to constrain key quantities such as mass-loss rates, dust properties, luminosities, and effective temperatures of stars using multi-band photometry in optical and IR bands.

A large variety of synthetic models are available for this purpose in the literature. Among the hydrostatic stellar atmosphere models used to describe stellar spectra are COMARCS (built specifically for AGB C-rich or O-rich stars, Aringer et al. 2016, 2019) or MARCS (Gustafsson et al. 2008). These models are subsequently combined with radiative transfer codes that incorporate dust features into the stellar spectral energy distribution (SED). Commonly used software include DUSTY (Ivezić & Elitzur 1997), often applied together with COMARCS models (e.g. Nanni et al. 2018, 2019; Groenewegen 2022), and 2Dust (Ueta & Meixner 2003), which provides the basis for the GRAMS models (Sargent et al. 2011; Srinivasan et al. 2011) constructed from MARCS and COMARCS atmospheres (e.g. Riebel et al. 2012; Boyer et al. 2013).

The great majority of the aforementioned studies investigated the properties of carbon stars in the Magellanic Clouds (e.g. Groenewegen et al. 2009, Riebel et al. 2012, Rau et al. 2017). One reason is due to the star formation history characterizing these two galaxies, which has produced numerous stars now observed as carbon stars stars (see Blanco & McCarthy 1983 for a detailed study of carbon stars as tracers of past star formation activity in the LMC). At the same time, the interpretation of the data is facilitated by the known distances and extinction, which are often much less well constrained within the Milky Way. Only a limited number of studies concentrate on determining stellar properties for the MW, with notable examples being Nanni (2019) and Groenewegen (2022). With the advent of Gaia, obtaining highly accurate distance measurements for Galactic stars has become possible, and has enabled more reliable extinction estimations and more accurate determinations of stellar parameters.

For this study we analysed the Gaia DR3 Golden Sample of Carbon Stars (Gaia Collaboration 2023a), and derived stellar and dust properties for its sources. This dataset includes carbon stars in the Milky Way and in the Magellanic Clouds, and allows a consistent treatment of stars formed in different galactic environments and with different metallicities. This approach enables us to investigate populations that differ in metallicity and evolutionary history within a single self-consistent framework, ensuring that the results obtained for the different systems can be directly compared without biases introduced by heterogeneous selection or analysis methods. In particular, we placed special emphasis on the Galactic sample, for which the availability of precise Gaia parallaxes and recent extinction estimates allows a reliable determination of stellar and dust parameters. Carbon stars can also be Mira variables, for which variability can strongly affect the construction of spectral energy distributions. To account for this, we modelled multi-band light curves for a subset of Miras, and we used the mean magnitudes derived from the best-fit models in place of single-epoch Gaia photometry, thus ensuring consistent and variability-corrected inputs for the analysis.

This paper is organized as follows. Section 2 describes the dataset and the photometric sources selected for this work. In Sect. 3 we detail the method applied to correct the photometry for variability effects in Mira variables. Section 4 introduces the grid of models adopted for the SED fitting, while Sect. 4.2 outlines the fitting procedure itself. We report our findings in Sect. 5 and compare them with the literature in Sect. 6, addressing the relevant caveats and physical implications. Finally, our conclusions are summarized in Sect. 7.

2. Dataset

Several catalogues of carbon stars have been compiled and published over the past 30 years, many of which are based on a combination of earlier studies. This has resulted in a complex landscape of partially overlapping datasets. Among the most recent efforts is the work by Suh (2024), which integrated multiple catalogues with stricter selection criteria to construct a refined sample of AGB carbon stars. Additional examples of recent C-star catalogues are, for instance, Abia et al. (2022), Chen & Yang (2012), and Li et al. (2018).

Despite efforts to standardize selection criteria, many of the catalogues discussed above still suffer from limitations in homogeneity. This study emphasizes homogeneity as a fundamental criterion by adopting the Gaia DR3 Golden Sample of Carbon Stars (GGSCS, hereafter) as a sourcelist, for which we infer both stellar and dust properties. In addition to being homogeneous, the GGSCS has extensive sky coverage, and high astrometric precision. The authors of this catalogue identified the stars based on the strength of C₂ and CN molecular bands in their Gaia XP spectra. Stars with colours bluer than G_BP − G_RP = 2 or fainter than G = 17.65 mag were excluded as detailed by the Gaia Collaboration (2023a). The catalogue includes a total 15 740 sources in the Milky Way and the Magellanic Clouds.

Figure 1 shows the sky distribution of the sources and their location on a JHK colour-magnitude diagram (CMD).

Fig. 1. Top panel: Sky distribution of the targets in the final catalogue. Bottom panel: Extinction-corrected colour-absolute magnitude diagram (J − Ks vs. MKs). The right subpanel shows the histograms of the absolute Ks magnitude probability density function (PDF) for the three populations. The colour-coding is as follows: SMC stars in green, LMC in red, and MW in blue.

2.1. Catalogue assembly

We complemented Gaia photometry with infrared data by matching the positions of the sources in the GGSCS with the 2MASS and AllWISE catalogues using a 1.0 arcsec search radius via the CDS Xmatch Service¹. We found 15 096 sources having a counterpart in both catalogues, and we retained 15 072 stars having all ten of the filter values available.

The three Gaia photometric filters are sensitive to stellar effective temperature (T_eff) and dust absorption, while the WISE bands are sensitive to dust emission and can trace chemical signatures, such as the silicon carbide feature at 11.3 μm in the W3 band, in the circumstellar envelopes. The 2MASS bands are sensitive to both the stellar T_eff and to dust emission. We note that Gaia and WISE magnitudes partially account for photometric variability through multi-epoch averaging, whereas this is not the case for 2MASS observations that may not accurately represent the average magnitude of variable sources. Section 3 discusses this aspect.

A number of studies similar to ours (e.g. Riebel et al. 2012; Groenewegen & Sloan 2018) also employed photometric data from Spitzer (Werner et al. 2004, filters I1, I2, I3, I4). We preferred not to use these data in order not to compromise the homogeneity of our sample, as only 41% of the stars have at least one Spitzer filter available.

We verified our positional crossmatch against the catalogue by Gavras et al. (2023), which uses a more advanced algorithm incorporating observational and astrophysical parameters beyond positional coincidence. A total of 12 679 out of 15 072 (84%) stars in our sample share identical counterparts with those reported by Gavras et al. (2023), confirming the reliability of our associations. The remaining sources are primarily located in high-density regions. Our final catalogue explicitly flags the sources found by Gavras et al. (2023).

2.2. Distances and extinction

Distances are a key parameter in our analysis, as they play a crucial role in both the extinction estimation and the SED fitting (see Sect. 4.2). We adopted the distances from Bailer-Jones et al. (2021) for stars not located in the MCs, acknowledging the inherent limitations of this approach, such as the adoption of Galactic priors and the single-star assumption, as discussed in detail in the original work. For stars in the LMC, we adopted a distance of 49.97 kpc (Pietrzyński et al. 2013), while for the SMC we used 62.44 kpc (Graczyk et al. 2020).

For all stars in our catalogue, we estimated integrated extinction values using the 3D dust extinction map of Lallement et al. (2022) combined with distance estimates from Bailer-Jones et al. (2021). Unlike asymptotic extinction, this approach traces the dust density along each line of sight, yielding a a more accurate determination of the foreground extinction, A_V, rather than relying on the total asymptotic values. The resulting A_V is then converted to other photometric bands (e.g. A_BP, A_K, A_W1) using the tabulated extinction coefficients provided by Wang & Chen (2019). These authors caution that these extinction coefficients were derived from red clump stars, and may not be appropriate for stars with a significantly different spectral energy distribution. We compared the extinction we derived with the values given by Abia et al. (2022), and found them to be compatible, suggesting that systematic uncertainties should be small enough not to impact the results and conclusions of the present analysis. The same extinction prescription was also adopted for the Magellanic Cloud sources. Although dedicated reddening maps for the Clouds are available (e.g. Skowron et al. 2021), we verified that the mean extinction values derived from the Lallement et al. (2022) map along the corresponding lines of sight are fully consistent with those reported by Skowron et al. (2021) for both the MCs. Therefore, the use of a homogeneous extinction treatment for all sources does not introduce significant systematic differences in our analysis.

For Galactic sources with larger inferred distances and the sources in the MCs, we therefore adopted the maximum integrated extinction along the line of sight provided by the map. As a consequence, the extinction for these objects may be slightly underestimated; however, this effect is expected to be small and does not qualitatively affect the results presented in this work.

2.3. Membership

To determine the membership of stars in the MCs, we applied a clustering algorithm to Gaia DR3 coordinates and proper-motions. To make it more resilient against outliers, we used the RobustScaler approach for scaling the data, following a similar application to open clusters by Hunt & Reffert (2023). We then applied the HDBSCAN algorithm (Ester et al. 1996), identifying 5235 LMC stars and 979 SMC stars. A total of 326 stars were found that do not belong to the MCs, but lack Bailer-Jones et al. (2021) distance estimates. Thus, we excluded these stars from our analysis. We also excluded 2 sources lacking extinction estimates. Our final sample consists of 14 747 sources, each with complete photometric coverage across all Gaia, 2MASS, and AllWISE filters, as well as reliable distance and extinction estimates.

2.4. Selection biases

The construction of the sample relies on Gaia detections and therefore unavoidably introduces selection effects. Due to its optical sensitivity, Gaia is biased against carbon stars affected by strong circumstellar extinction, such as sources undergoing very high mass-loss rates and surrounded by thick dust shells, which may fall below the Gaia detection threshold. At the same time, intrinsically faint carbon stars with low mass-loss rates may also remain undetected, particularly at the distance of the Magellanic Clouds, even in the absence of significant circumstellar extinction. The impact of these distance-dependent selection effects is clearly illustrated in Fig. 1. The three samples cover different ranges of of K-band absolute magnitude, with the SMC distribution showing a more pronounced truncation at the faint end compared to the MW, confirming that Gaia primarily detects only the brightest sources at larger distances. As a result, the sample may be incomplete at both the low and high ends of the mass-loss distribution, with these effects being substantially more severe in the Magellanic Clouds than in the Milky Way.

Additional selection effects are introduced by the use of WISE photometry. While WISE is essential for identifying dusty AGB stars, it is subject to sensitivity limits, source confusion, and saturation, especially at longer wavelengths and in crowded regions. These effects can affect both the faintest sources and the brightest objects with extreme infrared emission.

Consequently, the final catalogue should be regarded as representative of the carbon-star population within the dynamic range jointly accessible to optical and mid-infrared surveys. Because these selection effects differ between the Milky Way and the Magellanic Clouds, it is not straightforward to directly compare the samples. Therefore, to properly investigate metallicity-driven effects, it is necessary to account for differences in luminosity and evolutionary stage. We plan to address this specific analysis in a follow-up paper.

3. Average magnitude correction for Mira variables

Carbon stars are long-period variables (LPVs), often identified as Miras. For such objects, single-epoch photometry may not provide a reliable estimate of the mean magnitude, potentially introducing systematic biases in the derived stellar parameters.

To account for the variability of these stars, we analysed light curves from two main sources: Gaia DR3 epoch photometry (in the G, BP, and RP bands) and the OGLE III (Soszyński et al. 2009) and OGLE-IV (Iwanek et al. 2022) survey (in the V and I bands)². Although Gaia DR3 provides mean magnitudes based on multi-epoch observations, these are not derived from dedicated light curve fitting and may still be affected by variability-related biases. The OGLE data provide densely sampled optical light curves, making them particularly well suited for detailed variability analyses.

For all stars, we performed multi-band light curve fitting using a custom Python code developed for this work. The procedure builds on the LombScargleMultiband implementation (VanderPlas & Ivezić 2015), combined with additional publicly available tools for multi-band variability analysis. It simultaneously fits all available filters for each star, ensuring a consistent period, phase, and amplitude across bands. The resulting models provide an empirical, self-consistent description of the variability as a function of time. From the best-fit light curve models, we derived the mean magnitudes in the three Gaia bands (G, BP, and RP). These mean magnitudes are then adopted as input for the SED fitting procedure, described in Sect. 4.2. By construction, they are more robust than simple averages of individual epoch measurements, since they incorporate the full variability pattern across filters. For stars with OGLE data, we also used the V and I mean magnitudes derived from the same procedure to increase the reliability of the Gaia light curves fit. In the fitting procedure, we did not take into account Gaia photometric data with the columns rejected_by_photometry and variability_flag_g_reject as true.

We applied light curve corrections to a subset of stars classified as Miras, representing 18% of the total sample. These stars are regular pulsators, and the light curve fit is more reliable for this type of star than for the non-regulars. The identification of Miras was carried out by combining information from the OGLE catalogue with our Gaia-based analysis. For stars included in both Gaia and OGLE, we directly adopted the OGLE classification of Mira-type variables. To extend this classification to stars with only Gaia data, we examined the distribution of OGLE Miras in the diagram of period (P) versus variability amplitude in RP (AmpRP), finding that most occupy a region defined by an AmpRP > 0.2 mag and by 100 < P < 650 days. These empirical boundaries were then applied to the Gaia-only sources to identify additional Mira candidates. The final Mira subsample thus comprises both OGLE-flagged Miras and Gaia-only stars selected according to these criteria.

The periods derived from our light curve fits are fully consistent with OGLE values, confirming the reliability of our selection. Appendix A discusses in more detail the technical aspects of the algorithm adopted, as well as some light curve fitting examples.

To illustrate the effect of the correction, Fig. 2 compares the mean magnitudes provided by Gaia DR3 (Gaia Collaboration 2023b) to the values obtained from our multi-band light curve fitting. The left panel of Fig. 2 (BP band) exhibits a clear tilt, indicating that the applied correction increases for brighter BP magnitudes, while the effect is much weaker in RP and intermediate in the G band. This behaviour may reflect the strong optical variability of Mira variables, which is more pronounced in the BP band than in RP, combined with the fact that Gaia magnitudes are flux averages and do not account for the light curve shape. At fainter magnitudes, increased photometric noise and sparser effective sampling tend to wash out this effect, reducing the apparent correction.

Fig. 2. Absolute difference between the un-reddened GaiaBP, RP, and G magnitudes and our mean magnitudes derived from the light curve fit. The colours are the same as in Fig. 1. The y-axis range is restricted to highlight the region of interest, excluding a few outliers. The dashed black line shows the linear fit, and the shaded region denotes the 1σ confidence interval of the regression. The resulting fit parameters (slope and intercept) and their uncertainties are reported within each panel.

The correction is generally small, with absolute differences below 0.2 mag for most stars, although a few hundred objects show larger corrections, occasionally exceeding 1 mag. Corrections exceeding 1 mag occur only in the BP band (and partly in G), while they are much smaller in RP and negligible in the near-infrared. The total number of Miras identified is 2494. While these cases represent a small fraction of the sample, applying the correction ensures greater reliability of the SED fitting–both statistically, by improving the overall homogeneity of the dataset, and individually, by refining the estimates for Mira-type variables. The results of the light curve fits, such as the magnitude and the period found, are available in electronic format at the CDS.

4. Methods

In order to characterize the catalogue, we constructed a grid of spectral energy distribution models representative of carbon stars, and performed a fit them to the photometric data. In the present section, we describe first the models computation, and then the fitting method.

4.1. Model grid

In general, the observed SED of carbon stars results from the photospheric stellar emission being reprocessed by circumstellar dust. To model this, we adopted synthetic stellar spectra from the COMARCS³ library (Aringer et al. 2016, 2019), that provides hydrostatic stellar atmospheres specifically developed to study the structure and spectral properties of AGB stars. It is based on the MARCS code (Gustafsson et al. 2008) and includes updated molecular opacities and line lists, including for C-rich compositions. In this work, we employed all COMARCS models that satisfy the following criteria: M < 5 M_⊙, C/O ≥ 1, T_eff ≤ 5000 K, and T_eff ≥ 2500 K.

We used COMARCS photospheric spectra as input for DUSTY (V4)⁴ (Ivezic et al. 1999), a 1D radiative transfer code for modelling the radiation from astrophysical sources reprocessed by a dust shell. It solves the radiative transfer problem for spherical geometry using scaling properties that minimize free parameters, with dust temperatures determined from local radiative equilibrium. In addition to the photospheric spectrum, DUSTY requires information on the dust density structure to predict the output SED. In particular, it can either be computed self-consistently by solving the wind hydrodynamic equations (mode RDW), or an analytical profile can be specified to approximate the density profile. In principle, the former is more appropriate for AGB stars (Elitzur & Ivezić 2001; Nikutta et al. 2014; Höfner & Olofsson 2018), but it is considerably time-consuming. Therefore, the latter is often adopted in the literature (e.g. Groenewegen & Sloan 2018) due to its computational efficiency. In order to consistently compare our results with the literature, we also adopted the analytical profile, namely a r⁻² power law.

For similar reasons, we described the composition and optical properties of the dust in accordance with the recent literature (e.g. Groenewegen et al. 2009; Groenewegen & Sloan 2018; Riebel et al. 2012). We adopted a single chemical mixture made of amorphous carbon (amC). The optical constants have been taken from the ACAR sample provided by Zubko et al. (1996). For simplicity, we assumed a single sublimation temperature T_sub = 1200 K. The grain size is assumed to be 0.1 μm. For each COMARCS spectrum, we computed multiple DUSTY SEDs by varying the dust properties, namely the dust temperature and optical depth τ_V of the dust envelope. While the sublimation temperature of amorphous carbon is fixed at 1200 K, the inner dust temperature (T_dust) is an input parameter that defines the temperature at the inner radius of the dust shell. The latter, in particular, is a key parameter for determining the influence of dust on the observed properties of the stars. Its value correlates with the amount of dust around the star, thus with the IR excess of the light coming from the SEDs of the stars. Following Groenewegen (2022, and references therein) we assumed a fiducial wavelength of 0.55 μm for τ_V. We computed models for five different dust temperatures, from T_dust = 400 K to 1200 K in steps of 200 K, and for 55 values of the optical depth of the dust, ranging from 10⁻⁴ to 10. These values, together with the properties of the adopted COMARCS spectra, define the grid nodes, that are summarized in Table 1, for a total of ∼51 000 models.

The mass-loss rate Ṁ is computed using the equation from Gullieuszik et al. (2012), originally described by Groenewegen et al. 1998):

$$\begin{array}{c}\hfill {\tau }_V=5.405\times {10}^8\frac{\dot{M}{r}_{\mathrm{gd}}{Q}_V/a}{r_d{R}_{★}{v}_{\exp }{\rho }_d}·\end{array}$$(1)

Here Ṁ is the mass-loss rate in M_⊙ yr⁻¹, v_exp is the shell expansion velocity in km s⁻¹, and R_★ is the stellar radius in solar radii (R_⊙). Furthermore, ρ_d is the grain density in g cm⁻³, Q_V is the absorption efficiency factor in the V band, a is the dust grain radius in cm, r_gd is the dust-to-gas mass ratio, and r_d is the distance to the dust shell in stellar radii. For the calculation, we fixed v_exp = 10 km s⁻¹ (e.g. Groenewegen et al. 2009; Groenewegen & Sloan 2018), we calculated Q_V using the Python package miepython⁵, and we assumed a ρd = 2.2 g cm⁻³. We derived the optical depth τ_V from the best-fit DUSTY model. We assumed a constant dust-to-gas ratio of r_gd = 1/200 for the MW, LMC, and SMC, consistent with the standard approach adopted in the majority of comparable studies in the literature. The inner radius R_in is derived from the best-fit DUSTY model and subsequently rescaled accordingly to the stellar luminosity, since DUSTY computes this quantity assuming a reference luminosity of 10⁴ L_⊙. DUSTY outputs the inner radius of the dust shell in centimetres, which corresponds to the product of the stellar radius R_★ and the dimensionless parameter r_d. Thus, in Eq. (1), we calculated $R_{★}{r}_d=R_{\mathrm{in}}=R_{\mathrm{DUSTY}}{\left(\frac{L_{★}}{{10}^{4}L\odot }\right)}^{1/2}$. The luminosity of each star has been calculated by integrating its SED, and then multiplying it by 4πd².

The synthetic grid comprises over 50 000 models, with input parameters summarized in Table 1.

4.2. SED Fitting

The stellar and circumstellar parameters of the sources in our sample were derived by fitting the observed photometric data to synthetic model SEDs, with the best-fit model for each star determined by minimizing the reduced χ² value. The fitting procedure follows the classical χ² formalism,

$$\begin{array}{c}\hfill {\chi }_{\mathrm{red}}^2=\frac{1}{N-p}{\displaystyle \sum \frac{{[logF(\text{Obs},\lambda )-logF(\text{Model},\lambda )]}^2}{{\left(\frac{{\sigma }_{F(\text{Obs},\lambda )}}{F(\text{Obs},\lambda )ln10}\right)}^2}},\end{array}$$(2)

where F(λ) is the flux at a specific wavelength, N is the number of photometric data points for each source, p is the number of free parameters, and σ represents the uncertainty on that measurement. Synthetic fluxes are obtained by convolving the model spectra with the transmission profiles of each observational filter⁶.

To estimate uncertainties for the derived properties, we stochastically resampled the photometric data and repeated the SED fit 150 times for each star. For each star and filter, the observed magnitude was modelled as a Gaussian distribution centred on the measured value with a 1σ width equal to the photometric error reported in the corresponding reference catalogue. Random draws from these distributions were used to perform χ² fits, yielding a distribution of stellar parameters.

The final parameter values correspond to the medians of these posterior distributions, while the median-absolute-deviation times 1.4826 (one sigma in a Gaussian distribution) quantify their uncertainties. Using the median minimizes the impact of outliers, and we verified that median and mean values closely agree. This approach provides robust parameter estimates and a reliable assessment of their uncertainties. Figure 3 shows six examples of fitted SEDs.

Fig. 3. Examples of SED fits for six stars with different stellar and dust parameters. The red points indicate the observed fluxes normalized to the stellar luminosity, the blue line shows the best-fitting model, and the small black crosses represent the model fluxes convolved with the observational filter transmission curves, which are the quantities used in the SED fitting.

To assess the goodness of fit, we considered both the minimum χ² value and the average offset between the observed photometric data and the best-fitting model. Approximately 95% of the stars have χ² < 2000; only these are retained for our analysis and will be included in the published dataset.

We note that the absolute value of the reduced χ² in SED fitting cannot be interpreted in the classical statistical sense, as it strongly depends on the number of photometric points, on the adopted uncertainties, on the model grid resolution, and on the intrinsic variability of AGB stars. For this reason, χ² values from different studies cannot be directly compared. In our work, χ² is used exclusively as a relative internal metric to rank the goodness of fit among models for the same source.

Assuming that a perfect agreement between observed and model fluxes (i.e. an average flux ratio of unity) gives a ratio equal to one, we find that 97% of the stars show an average ratio $0.9955<\overline{\left(\frac{\mathrm{Obs}}{\mathrm{Model}}\right)}<1.0045\mathrm{erg}/\mathrm{s}{\mathrm{cm}}^2$, indicating generally good fits across the sample. At the end of the chi-squared process, we find a best-fit model for 13 777 sources (93% of the total). These stars compose the dataset that we discuss and publish.

5. Results

This section presents the main results of our SED fitting procedure. To compare the best-fit parameters for the Milky Way (MW), the Large Magellanic Cloud (LMC), and the Small Magellanic Cloud (SMC), we select, for each galaxy, a subsample composed of those stars that satisfy the following conditions:

$$\begin{array}{c}\hfill X_{\mathrm{C}\text{-}\mathrm{rich}}<W<1.7,X_{\mathrm{C}\text{-}\mathrm{rich}}=0.7+0.15{(M_K+7.65)}^2.\end{array}$$(3)

Here W is the Wesenheit index, defined in this work as the difference between the optical and near-infrared indices, W = W_{RP, BP − RP} − W_{Ks, J − Ks}. The individual components are defined as W_{RP, BP − RP} = G_RP − 1.3(G_BP − G_RP) and W_{Ks, J − Ks} = K_s − 0.686(J − K_s). The quantity X_C-rich defines the carbon star region in the Gaia-2MASS diagram (Lebzelter et al. 2018). We further adopted a threshold of W > 1.7 mag to identify the so-called extreme C-rich stars (Lebzelter et al. 2018; Nanni 2019). This classification, based on optical and near-infrared Wesenheit indices, is distinct from the mid-infrared selection of Extreme AGB stars proposed by Blum et al. (2006). This selection ensures that we compare stars with equivalent photometric properties across the MW, LMC, and SMC, allowing a consistent assessment of their best-fit parameters and any intrinsic differences among the three environments. In the following section, we refer to this subset as the C-region stars. Figure 9 shows the location of the C-region stars on the Gaia-2MASS diagram (carbon stars region). Table 2 lists the median values of the best-fit parameters for only the selected C-region stars in the three galaxies, followed by the relative dispersion. Table 3 shows a preview of the final catalogue, which provides the derived parameters for the entire sample of all stars analysed in this study.

5.1. Effective Temperature

The median effective temperature of the entire sample is T_eff = 3100 K, with a dispersion of σ T_eff = 300 K, consistent with the dominant population of carbon stars, as found in many literature studies (e.g. Abia et al. 2020). Focusing exclusively on the C-region stars subsample within the three populations, Table 2 summarizes the corresponding statistical properties, while Fig. 4 shows the distributions for these stars. Note that the individual measurement uncertainties (also listed in Table 2) can be primarily affected by the density of the COMARCS grid, which is not equally spaced. Although the median effective temperatures of the LMC and SMC samples are comparable within the uncertainties, their distributions show an extended high-T_eff tail, which is more prominent at lower metallicity. This behaviour is qualitatively consistent with stellar evolution expectations, according to which carbon stars in metal-poor environments can reach slightly higher effective temperatures (e.g. Blanco & McCarthy 1983; Nanni et al. 2021). The comparison between the two Magellanic Clouds and the Milky Way is subject to selection effects, since Gaia samples the carbon star populations differently in these environments, as described in Sect. 2.4.

Fig. 4. Distribution of Teff for C-region stars in the MW, LMC, and SMC. The colours and shapes of the three subsamples are shown in the legend.

5.2. Optical depth and mass-loss rate

Figure 5 presents the distribution of the mass-loss rates. These values are fully compatible between both MW and the MCs. The statistical properties of the mass-loss rates and optical depths for the C-region stars are reported in Table 2. As shown, both parameters are largely comparable across the three populations, with median mass-loss rates ranging from log Ṁ ≈ −7.8 to −8.2 M_⊙ yr⁻¹ and optical depths τ_V between 0.1 and 0.4. Although the MW sources appear to exhibit higher optical depths than their MCs counterparts, these differences are consistent within the estimated uncertainties. Therefore, these trends should be interpreted with caution, given that the distributions also remain sensitive to Gaia selection biases. We plan to address this issue in a follow-up paper, which will perform a more refined comparison accounting for differences in luminosity, evolutionary stage, and stellar mass.

Fig. 5. Distribution of mass-loss rates for stars in the MW, LMC, and SMC. The colours and shapes of the three subsamples are shown in the legend.

5.3. Mass

The stellar mass has only a minor influence on the overall shape of the SED, and our fitting procedure is therefore relatively insensitive to this parameter. Nevertheless, the derived median values reported in Table 2 are consistent with the expectations for carbon stars (1.5–4 M_⊙, Höfner & Olofsson 2018). Given the discrete nature of the grid, these values should be regarded as indicative rather than quantitatively significant. However, the slightly higher average masses inferred for the LMC and SMC are consistent with expectations from stellar evolution at lower metallicities, as determined by Pastorelli et al. (2020).

5.4. Carbon-to-oxygen ratio

The carbon-to-oxygen ratio (C/O) is another crucial parameter of the models; the values adopted here are those provided by the COMARCS model grid. On average, we observe higher C/O ratios in the Magellanic Clouds than in the Milky Way, with the median values presented in Table 2. Such a behaviour is qualitatively consistent with expectations for lower-metallicity environments, where the initial oxygen abundance is smaller and fewer dredge-up episodes are required to reach C/O > 1. In addition, stellar evolution models predict more efficient carbon enrichment during thermal pulses at low metallicity (Nanni et al. 2013; Di Criscienzo et al. 2013; Boyer et al. 2013; Abia et al. 2020). However, given the selection effects discussed in Sect. 2.4, this result should not be interpreted as direct evidence of a metallicity-driven effect.

Regarding the dust temperature at the inner boundary, for the majority of our entire sample, which is characterized by optically thin envelopes (τ_V ≲ 0.1), we determine values around ≈1000 K. We note, however, that for such low optical depths, the concept of a single dust temperature is less physically meaningful as the SED is not dominated by the dust emission from a single layer. Consequently, these values should be treated primarily as boundary conditions for the fit rather than intrinsic properties of the sources.

6. Discussion

This work has characterized the stellar and dust properties of the Gaia golden sample of carbon stars. As detailed in Sect. 2.4, the sample is biased against both dust-obscured and intrinsically faint sources. Since these selection effects impact the MW and MCs differentially, caution is required when performing cross-galaxy comparisons. In the following, we discuss our findings in the context of analogous studies from the literature.

In this work, we have also explored an alternative approach to the SED-fitting procedure, which adopts a systematic normalization of all the photometry bands to the 2MASS K_s filter. This method can ensure an SED-fit without the dependence by the star’s distance. Appendix B discusses the results of this approach, as well as their comparison to the standard SED-fitting method.

6.1. Effective temperature

As Sect. 5.1 discusses, we derived and compared the effective temperature distributions of MW and MCs carbon stars. A comparison of our results with the effective temperatures derived by Andrae et al. (2023) reveals a significant discrepancy: their median temperature is ≈5000 K, significantly higher than ours, with some stars exceeding 10 000 K. However, Andrae et al. (2023) noted that the GSP-Phot module can provide unreliable parameters for AGB stars. Therefore, our temperature estimates, derived specifically for these cool evolved stars, should be considered more reliable.

We further compared our results with previous SED-fitting studies of carbon stars, namely Groenewegen et al. (2009), Groenewegen & Sloan (2018), Groenewegen (2022), and Riebel et al. (2012). Figure 6 directly compares our results and with those of the literature. These works were selected because they share stars with our sample, adopt similar physical assumptions, and cover both the MW and the MCs, while differing in their treatment of circumstellar dust and radiative-transfer methods. In total, we identified 3 stars in common with Groenewegen et al. (2009), 21 with Groenewegen & Sloan (2018) and 3839 with Riebel et al. (2012), enabling a robust comparison.

Fig. 6. Comparison of the effective temperature values for the common sources between our work and those by Groenewegen et al. (2009) and Groenewegen & Sloan (2018) (left diagram) and Riebel et al. (2012) (right diagram). The data are shown using hexagonal density binning with a grid size of 40 and a threshold count of 10.

Overall, we find good agreement with these SED-fitting studies of carbon stars. For all stars in common, the differences in effective temperature are within ±500 K (see the left panel of Fig. 6), and a mild correlation is observed between the two T_eff estimates, with the coolest stars in our sample tending towards higher temperatures in the literature values.

All the aforementioned SED-fitting studies of carbon stars adopt COMARCS atmospheric models for the stellar component, but differ in their treatment of circumstellar dust and in the radiative-transfer codes employed: Riebel et al. (2012) used 2DUST, whereas the others relied on DUSTY or DUSTY-based codes (e.g. MoD; Groenewegen 2012). Differences in T_eff may therefore arise from several factors, including the choice of radiative-transfer code, grid resolution, and photometric datasets.

Additional sources of uncertainty include extinction and distance. Inaccurate distances can bias the best-fit models, as the fitting procedure compensates to reproduce the observed photometry and can also affect the inferred extinction. This effect is particularly relevant for MW stars, whereas distances and extinction values for MCs stars are well constrained and largely consistent across comparable studies. Finally, differences in the adopted dust optical constants can significantly modify the emerging SED shape, further affecting the derived effective temperatures.

6.2. Mass-loss rate

The mass-loss rate distributions of the Milky Way and the Magellanic Clouds show consistent results (see Sect. 5 and Fig. 5). For carbon stars, the dependence of mass-loss on metallicity is not straightforward. Since carbon is produced internally through dredge-up episodes, theoretical models and observations suggest that carbon stars can drive strong winds even at low metallicity (e.g. van Loon 2000; Groenewegen et al. 2007; Höfner & Olofsson 2018). The intrinsic mass-loss rate of C stars is primarily driven by the subsequent formation of amorphous carbon (amC) dust, the key mechanism responsible for wind acceleration (Nanni et al. 2021). From a theoretical perspective, metal-poor environments are expected to favour higher surface C/O ratios, owing to the lower initial oxygen abundance and to the increased efficiency of the third dredge-up at low metallicity. This condition can, in principle, enhance the circumstellar opacity and facilitate mass loss. However, such expectations apply to stars at comparable evolutionary stages and luminosities, and do not necessarily translate into observable differences when heterogeneous samples are compared.

Two main differences emerge between the MW and MC samples: considering the whole catalogue, the MW contains a larger number of stars with lower mass-loss rate and a slightly higher fraction of stars with higher mass-loss rate compared to the MCs (see Fig. 5). These differences are almost certainly driven by selection effects rather than intrinsic physical variations, as described in Sect. 2.4.

Moreover, as mentioned above, a robust assessment of metallicity-related effects on mass loss would require a comparison between stars matched in fundamental parameters, such as luminosity and evolutionary stage. Such an analysis, while enabled by the large sample and the availability of Gaia parallaxes, is beyond the scope of the present work.

We also compared our estimated mass-loss rate with the same set of studies used for the T_eff comparison (see Sect. 6.1). The results, shown in Fig. 7, do not reveal any clear systematic offset across the full range of mass-loss regimes.

Fig. 7. Comparison of the gas mass-loss rate values for the common sources between our work and those by Groenewegen et al. (2009) and Groenewegen & Sloan (2018) (left diagram) and Riebel et al. (2012) (right diagram). The uncertainties plotted are only the ones derived in this study. The cross in the bottom right corner of the right panel represents the average uncertainty. We note that the average uncertainty depends on the mass-loss regime, as is clear in the left panel.

In particular, when compared with the values from Riebel et al. (2012), we find an excellent agreement, with almost all stars lying within 1 dex of the zero-offset line. A comparison with Groenewegen et al. (2009), Groenewegen & Sloan (2018) also shows a generally good consistency, with an average offset of about 0.5 dex and no evident trends. We note that the apparent trend with mass-loss rate visible in Fig. 7 likely reflects the increasing sensitivity of SED-based Ṁ determinations to modelling assumptions at both low and high optical depths, rather than a genuine physical discrepancy. While a 1 dex difference represents an order of magnitude, such a spread is not unexpected in the context of AGB mass-loss modelling, given the high sensitivity of the results to the adopted model parameters. For instance, five stars in common exhibit, in our study, mass-loss rates that are 2–3 dex higher than the values reported in those studies (see Fig. 7). For these five objects, we find differences in bolometric luminosity of approximately 30–50% relative to our values. It is worth noting that SED fitting becomes intrinsically uncertain at low optical depth. In this regime, even modest variations in the adopted luminosity (or in any parameter that affects the SED normalization) can shift the inferred optical depth by a significant amount, which in turn propagates into large changes–up to several dex–in the derived mass-loss rate. This sensitivity is an important limitation in modelling low-τ carbon stars, and likely contributes to the offsets observed for these five objects. Furthermore, as demonstrated by Groenewegen & Sloan (2018), the choice of different optical constants for amorphous carbon grains (e.g. Rouleau & Martin 1991 vs. Zubko et al. 1996) can independently shift the derived Ṁ up to a factor of several.

Nanni et al. (2019) provide a more systematic assessment of the possible caveats, highlighting how several modelling assumptions can easily bias the final estimates. Beyond the choice of the dust optical constants, other ingredients of the models may significantly affect the outcome, such as the adopted initial expansion velocity or the neglect of dust–gas drift, which can alter the relation between gas and dust densities. Another fundamental ingredient that can have a strong impact on the results, in particular on the derived mass-loss rates, is the treatment of the wind dynamics. Even the assumed coefficients for amorphous carbon, as mentioned before, or the prescription for the seed nuclei abundance, can play a role. Additional uncertainties come from possible deviations from spherical symmetry of the model, as well as from photometric variability, which introduces extra scatter in the SED fitting and is taken into account in this work (see Sect. 3 and Appendix A). Taken together, these factors imply that mass-loss rates and related quantities should always be considered within the context of the adopted modelling framework, rather than as absolute values. Possible discrepancies might also be attributed to differences in the SED fitting parameters or wavelength coverage.

In addition to the star-by-star comparison with literature values, we examined the global behaviour of our results in the luminosity–mass-loss rate (L–Ṁ) diagram, shown in Fig. 8. For the LMC, SMC, and MW, we find that our stars (black points) occupy the same range of luminosities and mass-loss rates as reported by previous studies, confirming the overall consistency of our determinations. A slight selection bias from the Gaia sample is evident, however, as both the MCs and MW samples do not reach the most extreme mass-loss regimes, as the stars studied by Groenewegen et al. (2009), Groenewegen & Sloan (2018) and Groenewegen (2022). Moreover, our mass-loss estimates are consistent with expectations from theoretical stellar evolution models (see Pastorelli et al. 2019 and Pastorelli et al. 2020).

Fig. 8. From left to right: MW, LMC, and SMC stars. Luminosity against mass-loss rate values for stars in our sample (black dots), Groenewegen (2022, red triangles), Groenewegen et al. (2009, blue pentagons), Groenewegen & Sloan (2018, green squares), and Riebel et al. (2012, purple crosses). The black point in the bottom right corner of each panel represents the average uncertainty of our estimated mass-loss rates for each galaxy.

6.3. Gaia-2MASS diagram

The Gaia–2MASS diagram (see Sect. 5 for its definition) provides further insight into the catalogue adopted and our findings. Introduced by Lebzelter et al. (2018), it combines optical and near-infrared photometry to provide an efficient diagnostic for the chemical type and evolutionary status of AGB stars. Figure 9 shows this diagram for our entire sample, colour-coded by the mass-loss rate derived from SED fitting. Magellanic Clouds carbon stars populate the expected region, with the dustiest and most evolved objects–characterized by higher Ṁ values–progressively shifting towards redder colours, and in general all occupying the same region in the diagram.

Fig. 9. From left to right: MW, LMC, and SMC stars. Gaia-2MASS diagram coloured by the optical depths τV values derived in our study. The dotted lines represent the observational-derived lines delimiting the different regions occupied by different types of AGB stars, and are defined as in Lebzelter et al. (2018). In the left panel, each region has been labelled accordingly with red text. The stars inside the C-rich region are the ones defined in this paper as C-region stars.

The MW sample behaves similarly. The majority of MW stars are located in the same carbon-rich region as their MC counterparts, confirming their C-rich classification. However, the MW population shows a broader spread, including a tail towards redder Wesenheit values with higher Ṁ, suggesting the presence of more extreme-carbon stars, which are characterized by an Ṁ > 10⁻⁶ M_⊙ yr⁻¹ (Nanni et al. 2019). Additionally, several MW sources occupy a low-luminosity region where we expect to see a mix of faint AGB and RGB stars. These sources likely correspond to less evolved, dust-poor carbon stars, whose photometric properties have yet to be significantly affected by mass loss. This overall distribution confirms that the Gaia–2MASS diagram not only distinguishes between O-rich and C-rich AGB populations but also traces evolutionary effects such as dust formation and mass-loss enhancement. In addition, we highlight that the large size of our sample allows us to statistically trace the general distribution of the dust and stellar parameters of these carbon stars. The grid-like appearance of our results in Fig. 8 reflects the discrete sampling of optical depth values adopted in our model grid. This does not affect the global trends, but it explains why our points appear more regularly spaced than those from other studies. A rigorous decontamination would require a dual-grid fitting approach (C-rich vs. O-rich), which is beyond the scope of this work.

6.4. Variability of the carbon stars

Since the pulsations of carbon stars are the main driver of their mass-loss processes, we expect a connection between their photometric variability and the circumstellar dust content. We find a clear correlation between the mass-loss rate derived from our models and the variability amplitude measured in the GaiaG band, in the sense that stars with larger amplitudes tend to show higher Ṁ. The fitted relation we derived is G_Amplitude = s · log Ṁ + h, with s = 1.17 ± 0.1 and h = −7.67 ± 0.1. In the linear fit, measurements with lower uncertainties were assigned higher weights. At the same time, we note that the dispersion in mass-loss rate at a given pulsation amplitude is substantial, reaching approximately two orders of magnitude, higher than the intrinsic uncertainty. This large spread indicates that pulsation amplitude alone is not sufficient to uniquely determine the mass-loss rate. Other stellar properties, such as luminosity, evolutionary stage (the progression along the TP-AGB), and chemical composition, are likely to contribute significantly to the observed scatter. Within these limitations, the observed trend likely reflects the fact that stronger pulsations can facilitate the injection of material into the stellar envelope, thereby favouring dust formation and enhanced mass loss. A similar qualitative behaviour has been reported for oxygen-rich Mira variables, where optical variability amplitudes correlate with the intensity of SiO maser emission, interpreted as a consequence of pulsation-driven radiative pumping (Alcolea et al. 1990). Figure 10 illustrates both the overall trend and the large intrinsic scatter in the relation between mass-loss rate and G-band amplitude.

Fig. 10. Scatter plot of Gaia G mean magnitude derived from our light curve fit method (see Sect. 3) and mass-loss rates from our SED fit. The red dotted line is a linear fit between amplitude and mass-loss.

6.5. Expectation from theoretical models

Aringer et al. (2019) showed that spectra and colours of C stars are most sensitive to the C/O ratio rather than to the overall metallicity, and they are also strongly affected by T_eff and molecular line saturation. Spectral markers of nitrogen and oxygen enrichment (e.g. CN, C₂, and CO features) also depend on T_eff, becoming more evident at higher T_eff values. In particular, at higher T_eff, variations in O and N abundances have a stronger and more observable impact on the spectra, as molecular absorption becomes less saturated. These trends broadly agree with expectations that C-star spectra are less sensitive to metallicity and T_eff than M star spectra, which are dominated by metallicity-linked molecular species like TiO.

7. Conclusions

In this work we have presented a homogeneous SED-fitting analysis of a large sample of carbon stars in the Milky Way and the Magellanic Clouds, combining Gaia DR3, 2MASS and WISE photometry with the radiative transfer code DUSTY. We derived effective temperatures, luminosities, optical depths and mass-loss rates for 14 747 sources. The results confirm the general picture of carbon-rich TP-AGB stars as efficient mass-losing objects, with typical average effective temperatures at 3150 K and a standard deviation of 430 K. Mass-loss rate values span the range 10⁻¹¹–10⁻⁴ M_⊙ yr⁻¹. Specifically regarding the LMC, we find a distribution of mass-loss rates consistent with the theoretical predictions calibrated for this galaxy by Pastorelli et al. (2020).

We find that stars with larger variability amplitudes generally show higher mass-loss rates, supporting the hypothesis that pulsations play a major role in levitating material into the dust-formation zone. However, as noted in the analysis (see Sect. 6.4), a large spread in Ṁ for a given amplitude is present, indicating that other factors also play a significant role in driving the mass loss. The distributions of effective temperature and mass loss differ across environments. In particular, the lower-metallicity samples exhibit an extended high-T_eff tail, rather than a strong systematic shift in the median values, in qualitative agreement with expectations from metallicity-dependent stellar evolution. However, selection biases may affect this result. The comparison with previous studies reveals broad agreement in the trends but also differences of up to some orders of magnitude in the absolute values of Ṁ. These discrepancies reflect the strong sensitivity of the derived quantities to the adopted physical assumptions.

As emphasized by Nanni et al. (2019), the determination of mass-loss rates through SED fitting inevitably depends on the models. The choice of optical constants for amorphous carbon grains, the treatment of wind dynamics, or the assumed gas-to-dust ratio and expansion velocity can significantly affect the derived Ṁ. Additional uncertainties arise from spherical symmetry and from photometric variability. A major role can also be played by the treatment of the dynamic interaction between gas and dust around the AGB star. Nevertheless, the homogeneous methodology adopted in this work ensures that all sources have been analysed within a consistent modelling framework. This enables reliable relative comparisons among stars in different environments, even if the absolute values depend on the adopted physical prescriptions.

Despite these caveats, our analysis delivers an extensive and homogeneous catalogue of stellar and circumstellar parameters for carbon stars in the Gaia DR3 Golden Sample. These results provide a valuable reference for future studies of carbon-star properties based on Gaia data. It may also serve to refine the criteria for the identification and classification of carbon stars in forthcoming Gaia data releases, including DR4 and DR5.

Future progress will benefit from extending the fitting procedure to larger and more diverse samples, ideally encompassing all available carbon-star catalogues within a consistent and homogeneous framework. It will also be essential to use machine-learning methods to create continuous model grids, in order to reduce the dependence of the results on the grid itself, and to adopt methods that account for the time-dependent behaviour of pulsating AGB stars. In this context, fitting SEDs at different phases of the pulsation cycle would enable a more complete characterization of the variability in stellar parameters. A particularly promising direction involves the development of time-dependent hydrodynamic models that include the effects of pulsation, together with the application of Bayesian statistical techniques to quantify uncertainties and correlations among the derived quantities. Moreover, improved distances and extinction values will play a crucial role in the determination of the parameters, especially for the MW sample. These advances will enable a more robust use of AGB carbon-star samples in studies of stellar evolution and dust production in galaxies.

Data availability

The results presented in Table 3 is available at the CDS via https://cdsarc.cds.unistra.fr/viz-bin/cat/J/A+A/708/A120. The ATHENA-C model grid, described in Sect. 4 (or Table 1), is openly available on Zenodo at https://doi.org/10.5281/zenodo.18418537

Acknowledgments

This project has received funding from the European Union’s Horizon 2020 research and innovation programme under the Marie Skłodowska-Curie grant agreement No 101072454, as part of the Milky Way-Gaia Doctoral Network (https://www.mwgaiadn.eu/). The authors, especially A.L., would like to thank S. Srinivasan for the important and useful conversations regarding the statistical treatment of the results and uncertainties; A. Bonanos for her guidance and support during the course of the project; and C. Kalup for the important discussion on the complex mosaic of available Galactic extinction maps. The contribution by MT and GP is funded by the European Union – NextGenerationEU and by the University of Padua under the 2023 STARS Grants@Unipd programme (“CONVERGENCE: CONstraining the Variability of Evolved Red Giants for ENhancing the Comprehension of Exoplanets”). GC acknowledges financial support from European Union–Next Generation EU, Mission 4, Component 2, CUP: C93C24004920006, project ‘FIRES’. This work has made use of data from the European Space Agency (ESA) mission Gaia (https://www.cosmos.esa.int/gaia), processed by the Gaia Data Processing and Analysis Consortium (DPAC, https://www.cosmos.esa.int/web/gaia/dpac/consortium). Funding for the DPAC has been provided by national institutions, in particular the institutions participating in the Gaia Multilateral Agreement. This research has made use of the VizieR catalogue access tool (Ochsenbein et al. 2000), CDS, Strasbourg Astronomical Observatory, France (DOI: 10.26093/cds/vizier). This research has made use of the Spanish Virtual Observatory (https://svo.cab.inta-csic.es) project funded by MCIN/AEI/10.13039/501100011033/ through grant PID2020-112949GB-I00 (Rodrigo & Solano 2020). This research has used data, tools and materials developed as part of the EXPLORE project that has received funding from the European Union’s Horizon 2020 research and innovation programme under grant agreement No 101004214. We wish to thank the “Summer School for Astrostatistics in Crete” for providing training on the statistical methods adopted in this work.

References

Appendix A: Light curve analysis

To analyse the variability of our sample stars, we developed a dedicated Python code that processes time-series photometry and reconstructs phase-folded light curves across multiple filters. The code follows a modular design, ensuring that each step of the analysis is fully reproducible and easily adaptable to different datasets. The input consists of Gaia DR3 epoch photometry and, when available, complementary data from OGLE. These measurements are assembled into a dedicated TimeSeries object that standardizes the data structure by storing observation times, magnitudes with their associated uncertainties, and the corresponding photometric bands. This approach enables the consistent treatment of sources with heterogeneous sampling and multi-band observations within a unified framework.

The first task of the code is the determination of the variability period. For this purpose, we implemented the multi-band LombScargle periodogram, which simultaneously fits all available photometric bands and maximizes the likelihood of a common variability period, while allowing for differences in amplitude and phase among the filters. The algorithm provides as output the best-fit period, its uncertainty, and the complete periodogram for diagnostic purposes. For those stars in common with the OGLE survey, we verified that the periods derived by our code are fully consistent with the OGLE values, thus confirming the reliability of our implementation.

Once the period is determined, the code reconstructs the light curve model. The data are folded with the derived period and fitted with a multi-harmonic Fourier series across all available bands. This step allows us to produce both a global model of the variability and individual-band reconstructions. Figure A.1 shows an example of a light curve multi-band reconstruction. The resulting models are used to investigate the variability properties of the stars and to provide homogeneous light curve products. We decided to fit with one single period to each star in order to be more accurate for the Miras light curves. The code also has the option for a multi-periodic.

Fig. A.1. Observed multi-band and multi-period photometry for the star Gaia DR3 2935152926080940544 with the best-fit light curve models. The points represent the observed magnitudes, and the curves represent the best-fit light curve for each photometry band. The horizontal dashed lines are the mean magnitudes for each best-fit model of each band. To create the multi-band light curve model for this star, we used OGLE bands I (purple) and V (yellow), and Gaia bands G (green), BP (blue), and RP (red).

The outputs include the best-fit periods with uncertainties, the model light curves for each band, phase-folded plots of the data and models, classification flags (e.g. Miras, other LPVs), and diagnostic statistics such as χ² and residual scatter. These products are collected in structured tables and figures, which form the basis of the final catalogue used in this work.

Appendix B: Testing the chi-squared methodology

During the SED-fitting analysis, we tested several alternative approaches to assess their impact on the derived stellar parameters. Among these, the most relevant test concerned the normalization of both observed and model fluxes to the K band. In this approach, the K band is fixed, and the fitting becomes sensitive primarily to the shape of the SED rather than to its absolute scaling. Yang et al. (2023) provide an example of an SED fitting study that adopts this methodology.

The results presented in Fig. B.1 show that the effective temperature distribution obtained with the K-band normalization is compatible, within the uncertainties, to the one derived with classical method. The only notable difference in the temperature distribution is the presence of hotter stars with the classical method. The same stars studied with the K-normalization method have a temperature compatible with the bulk of the sample. The normalization tends to yield lower optical depths and, consequently, lower mass-loss rates compared to the classical method, described in Sect. 5. In particular, as shown in Fig. B.2, it is clear that K-normalized results generally provide lower mass-loss values and are more sparse towards the low-end of the mass-loss distribution, while the results derived with the general classical method provide more self-consistent results at higher mass-loss rates. This is expected, since the mass-loss rate, which is derived from the optical depth, depends strongly on the IR excess of the stellar flux. When the SED is normalized, sensitivity to the absolute IR excess is reduced, while sensitivity to the spectral shape at shorter wavelengths increases. In contrast, the T_eff values are absolutely comparable.

Fig. B.1. Histogram of distribution of the effective temperature for only LMC stars, derived by the classical SED fitting that we applied in this work (red distribution) and the SED fitting using both observations and models normalized to the K filter (black).

Fig. B.2. Histogram of distribution of the mass-loss rates for only LMC stars, derived by the classical SED fitting that we applied in this work (red distribution) and the SED fitting using both observations and models normalized to the K filter (black).

All Tables

All Figures

	Fig. 1. Top panel: Sky distribution of the targets in the final catalogue. Bottom panel: Extinction-corrected colour-absolute magnitude diagram (J − Ks vs. MKs). The right subpanel shows the histograms of the absolute Ks magnitude probability density function (PDF) for the three populations. The colour-coding is as follows: SMC stars in green, LMC in red, and MW in blue.
In the text

Fig. 2. Absolute difference between the un-reddened GaiaBP, RP, and G magnitudes and our mean magnitudes derived from the light curve fit. The colours are the same as in Fig. 1. The y-axis range is restricted to highlight the region of interest, excluding a few outliers. The dashed black line shows the linear fit, and the shaded region denotes the 1σ confidence interval of the regression. The resulting fit parameters (slope and intercept) and their uncertainties are reported within each panel.

In the text

	Fig. 3. Examples of SED fits for six stars with different stellar and dust parameters. The red points indicate the observed fluxes normalized to the stellar luminosity, the blue line shows the best-fitting model, and the small black crosses represent the model fluxes convolved with the observational filter transmission curves, which are the quantities used in the SED fitting.
In the text

	Fig. 4. Distribution of Teff for C-region stars in the MW, LMC, and SMC. The colours and shapes of the three subsamples are shown in the legend.
In the text

	Fig. 5. Distribution of mass-loss rates for stars in the MW, LMC, and SMC. The colours and shapes of the three subsamples are shown in the legend.
In the text

	Fig. 6. Comparison of the effective temperature values for the common sources between our work and those by Groenewegen et al. (2009) and Groenewegen & Sloan (2018) (left diagram) and Riebel et al. (2012) (right diagram). The data are shown using hexagonal density binning with a grid size of 40 and a threshold count of 10.
In the text

Fig. 7. Comparison of the gas mass-loss rate values for the common sources between our work and those by Groenewegen et al. (2009) and Groenewegen & Sloan (2018) (left diagram) and Riebel et al. (2012) (right diagram). The uncertainties plotted are only the ones derived in this study. The cross in the bottom right corner of the right panel represents the average uncertainty. We note that the average uncertainty depends on the mass-loss regime, as is clear in the left panel.

In the text

Fig. 8. From left to right: MW, LMC, and SMC stars. Luminosity against mass-loss rate values for stars in our sample (black dots), Groenewegen (2022, red triangles), Groenewegen et al. (2009, blue pentagons), Groenewegen & Sloan (2018, green squares), and Riebel et al. (2012, purple crosses). The black point in the bottom right corner of each panel represents the average uncertainty of our estimated mass-loss rates for each galaxy.

In the text

Fig. 9. From left to right: MW, LMC, and SMC stars. Gaia-2MASS diagram coloured by the optical depths τV values derived in our study. The dotted lines represent the observational-derived lines delimiting the different regions occupied by different types of AGB stars, and are defined as in Lebzelter et al. (2018). In the left panel, each region has been labelled accordingly with red text. The stars inside the C-rich region are the ones defined in this paper as C-region stars.

In the text

	Fig. 10. Scatter plot of Gaia G mean magnitude derived from our light curve fit method (see Sect. 3) and mass-loss rates from our SED fit. The red dotted line is a linear fit between amplitude and mass-loss.
In the text

Fig. A.1. Observed multi-band and multi-period photometry for the star Gaia DR3 2935152926080940544 with the best-fit light curve models. The points represent the observed magnitudes, and the curves represent the best-fit light curve for each photometry band. The horizontal dashed lines are the mean magnitudes for each best-fit model of each band. To create the multi-band light curve model for this star, we used OGLE bands I (purple) and V (yellow), and Gaia bands G (green), BP (blue), and RP (red).

In the text

	Fig. B.1. Histogram of distribution of the effective temperature for only LMC stars, derived by the classical SED fitting that we applied in this work (red distribution) and the SED fitting using both observations and models normalized to the K filter (black).
In the text

	Fig. B.2. Histogram of distribution of the mass-loss rates for only LMC stars, derived by the classical SED fitting that we applied in this work (red distribution) and the SED fitting using both observations and models normalized to the K filter (black).
In the text