Mass and wind luminosity of young Galactic open clusters in Gaia DR2

Context. Star clusters constitute a significant part of the stellar population in our Galaxy. The feedback processes they exert on the interstellar medium impact multiple physical processes from the chemical to the dynamical evolution of the Galaxy. In addition, young and massive stellar clusters might act as efficient particle accelerators and contribute to the production of cosmic rays. Aims. We aim at evaluating the wind luminosity driven by the young (<30 Myr) Galactic open stellar clusters observed by the Gaia space mission. This is crucial for determining the energy channeled into accelerated particles. Methods. To do this, we developed a method relying on the number, magnitude, and line-of-sight extinction of the stars observed per cluster. Assuming that the stellar mass function follows a Kroupa mass distribution and accounting for the maximum stellar mass allowed by the age and mass of the parent cluster, we conservatively estimated the mass and wind luminosity of 387 local clusters within the second data release of Gaia. Results. We compared the results of our computation with recent estimates of young cluster masses. With respect to these, our sample is three times more abundant, particularly above a few thousand solar masses. This is of the utmost relevance for predicting the gamma-ray emission resulting from the interaction of accelerated particles. The cluster wind luminosity distribution we obtained extends up to 3x10^38 erg/s. This is a promising feature in terms of potential particle acceleration scenarios.


Introduction
Stellar clusters belonging to the disk of the Galaxy are usually referred to as open clusters.They consist of homogeneous groups of stars with the same age and initial chemical composition.As star formation occurs primarily in clusters, they constitute ideal laboratories for tracing the star formation and evolution in our Galaxy (see Krumholz et al. 2019, for a recent review).The dynamical evolution of clusters is a result of several processes, including stellar evolution, two-body relaxation, tidal interactions, and shocks.During their lifetime, clusters tidally interact with their parent giant molecular clouds and with the Galactic structure.They eventually dissolve after a few relaxation times.
Young clusters, in particular, are key to understanding star formation processes because they enable the initial mass function (IMF) to be traced.The most massive of these objects, so-called young and massive stellar clusters (YMSCs), have recently received attention in the high-energy astroparticle physics community because they have been detected in very high-energy (VHE, 100 GeV) gamma rays (Aharonian et al. 2019) with the High Energy Stereoscopic System (H.E.S.S.).
Table C.1 is available at the CDS via anonymous ftp to cdsarc.cds.unistra.fr (130.79.128.5) or via https: //cdsarc.cds.unistra.fr/viz-bin/cat/J/A+A/686/A118Morphological measurements of this emission have enabled the derivation of the spatial distribution of accelerated particles in the surroundings of several Galactic clusters (including Westerlund 1 and the Cygnus Cocoon).The gamma-ray profiles and intensities suggest a continuous injection that may be powered by the stellar winds of massive stars provided that ∼1% of the wind power is converted into accelerated particles.The Cygnus Cocoon has even been observed at the highest energies ever achieved in gamma-ray astronomy: two photons with an energy above 1 PeV (1 PeV = 10 15 eV) have been detected from the Large High Altitude Air Shower Observatory (LHAASO; LHAASO Collaboration 2024).These observations probed the capabilities of stellar clusters to operate as PeV-particle accelerators (i.e., PeVatrons) and beyond.
Explaining the origin of cosmic rays (CRs) in the PeV domain remains one of the challenges in astroparticle physics today.While supernova remnants (SNRs) are often addressed to explain their energetics, the maximum energy that can be produced by these sources remains unclear (Celli et al. 2019;Cristofari et al. 2020;Diesing 2023).On the other hand, the power provided by stellar winds can be 10% of the power provided by SNRs (Cesarsky & Montmerle 1983;Morlino et al. 2023).Moreover, YMSCs have favorable conditions in terms of shock size and duration, so that it is plausible that energies up to ∼PeV can be obtained (Morlino et al. 2021).For these reasons, YMSCs were suggested as candidate sources of Galactic CRs up the knee region, complementary to SNRs.In addition to energetic considerations, acceleration of particles from the wind of massive stars appears to be a necessary component to explain the 22 Ne/ 20 Ne anomaly in CR composition (Casse & Paul 1982;Prantzos 2012;Tatischeff et al. 2021): the enhancement observed with respect to solar abundances requires an acceleration from a carbon-enriched medium rather than from the standard interstellar medium (ISM), such as from the winds of massive stars.The relative contribution of these different source populations, SNRs and YMSCs, to the observed CRs is yet to be clarified, as is their role in energy, and in particular, in the region of about a few PeVs (Vieu & Reville 2023;Morlino et al. 2023).
Furthermore, the exact location where particle acceleration would take place in YMSCs is a topic of active debate in the community.Several proposals have been put forward, based either on the acceleration in wind-wind collision regions (Bykov et al. 2020), on stochastic second-order acceleration in the highly turbulent medium of the cluster core (Klepach et al. 2000), or on acceleration at the wind termination shock (WTS; Morlino et al. 2021).When a cluster is older than a few million years, supernova (SN) explosions start to dominate the dynamics (Vieu et al. 2022), but for younger clusters, the energy input available for particle acceleration can only come from stellar winds.Understanding their exact kinetic luminosity is therefore of paramount importance.It is therefore necessary to characterize the cluster physical properties, in particular, their mass and wind luminosity.This is the objective of this work, while in two accompanying papers, we will explore the acceleration model and its radiative signatures.In all of the three works, we focus on young clusters for the reasons explained above.
The Gaia satellite has recently performed the most accurate astrometric and photometric survey of the Milky Way.This extended our knowledge about the Galaxy structure and evolution.The application of advanced classification algorithms based on machine-learning techniques to its second data release (DR2) has resulted in the identification of an extended sample of clusters (Cantat-Gaudin et al. 2020).Even though this catalog does not correspond to the latest available data release, it nevertheless contains detailed information concerning cluster member stars and their physical parameters, which are crucial (in the context of the WTS particle acceleration model) to evaluate the luminosity of the collective wind.
This paper is organized as follows.In Sect.2, we introduce and characterize the properties of the open cluster data sample we adopted for this study, which is devoted to systems younger than 30 Myr.We aim to provide all the necessary ingredients for the computation of the particle acceleration spectrum produced at the cluster WTS.We therefore describe the method we developed for estimating the mass and wind power of each cluster, for which we relied on the properties of the observed cluster star members reported by Gaia.In particular, we obtain in Sect. 3 a lower limit to the values of cluster masses by assuming the same stellar mass distribution function for all young clusters and normalizing it to reproduce the Gaia observations in terms of the number of stars within a given magnitude range.We further account for the optical extinction in the cluster direction.We compare this result with existing cluster mass estimates in Sect. 4. A fraction of the cluster sample investigated here was also reported in the Milky Way Stellar Cluster (MWSC) catalog and mass values derived by Just et al. (2023) through their tidal parameters.Moreover, the recent study by Almeida et al. (2023) yielded mass values for some of the clusters belonging to the same Gaia data release we investigated from an approach based on a Monte Carlo method.We hence discuss the compatibility of the different results.In Sect.5, we present a procedure to evaluate the related mass uncertainties based on the number of observed WR stars in Galactic open clusters.By including the contribution of observed WRs, we further derive the cluster wind luminosities in Sect.6.Finally, we conclude in Sect.7 and additionally provide three appendices that discuss technical aspects of our cluster mass estimation approach, the energetic contribution of SNe occurring in the cluster sample under investigation, and details about the CDS online supplementary data file.

Data sample and methods
In its DR2, 2017 open clusters have been identified in the Gaia data (Cantat-Gaudin et al. 2020).For 1867 of these, cluster parameters such as the distance modulus, extinction, and age were obtained through artificial neural network techniques.Compared to previous surveys, such as the MWSC catalog (Kharchenko et al. 2013), the sample of open clusters from Gaia DR2 reports detailed 3D information that is only allowed by astrometric measurements (proper motions and parallaxes).Therefore, we decided to adopt the Gaia DR2 open cluster catalog as a reference for our study, and to complete missing information about clusters that is required by the physical modeling from the MWSC survey.
An important feature of the Gaia release is that it provides values for the cluster age.This is a crucial parameter for the application of acceleration models because the main source of energy in these systems changes as they evolve.The energetics of star clusters is dominated by stellar winds for the first few million years and subsequently by SNR shocks (Vieu et al. 2022).Therefore, we applied a selection criterion to the initial cluster sample and only retained those younger than 30 Myr, which roughly corresponds to the main-sequence (MS) lifetime of stars with M 8 M , that is, the minimum mass value for stellar evolution to terminate with an SN explosion.This reduced the number of considered clusters to 390.
For each cluster in the selected sample, we derived the wind luminosity, which strongly depends on the cluster evolutionary stage and composition in terms of its member stars.However, because of the limitations due to instrumental sensitivity and light extinction, the actual stellar population of each cluster had to be inferred from the observed one.
Our approach is based on the assumption of a stellar mass distribution inside the cluster, ξ(M) ≡ dN/dM, for which we derived the correct normalization to reproduce the number of detected stars per cluster, N * , in the given magnitude range of the observations (for further details, see Sect.3).To do this, (i) we first extracted the number and magnitude of the observed member stars for every young cluster reported by Cantat-Gaudin et al. (2020), as well as the extinction toward the corresponding direction; (ii) we then obtained the observed minimum and maximum G-band magnitude of its member stars through a kernel density fit of the stellar magnitude distribution of each cluster; (iii) we converted these values into the maximum and minimum intrinsic luminosity of member stars, respectively, through both bolometric and extinction corrections; (iv) we later assumed the mass-luminosity relation of MS stars to derive the expected maximum and minimum stellar mass observed per cluster, which we labeled M * max and M * min , respectively; A118, page 2 of 13 (v) we finally derived the normalization of the stellar mass function for each cluster to reproduce the observed number of stars N * in the derived mass range [M * min ; M * max ].Since the correction for G-band extinction is based on the measured stellar effective temperature, the lack of this parameter in three clusters of the sample (i.e., IC 2948, NGC 1333, and Pismis 16) reduced the size of the final cluster sample to 387.The distribution of their age and distance from Earth are shown in Fig. 1, with median values of 14.5 Myr and 2.3 kpc, respectively.The spatial distribution of our sample (panel b) is almost flat up to ∼3 kpc with large fluctuations, while it decreases at larger distances.This likely indicates that the sample is complete up to this distance.The fluctuations are probably due to the local inhomogeneities.The actual distance up to which the Gaia sample can be considered complete is not known with accuracy: The analysis of the extinction maps (Lallement et al. 2022) suggests that clusters should be detectable up to 3 kpc, even though there are directions in the Galactic plane where the extinction is high even at closer distances.
Figure 1 also reports the number of stars associated with each cluster, from which we derived the correct normalization of the mass distribution.After this was obtained, we derived the current stellar mass function of each cluster by accounting for stellar luminosities beyond the Gaia sensitivity following theoretical considerations.In particular, we fixed an absolute minimum stellar mass M min to 0.08 M , corresponding to the limiting mass value in brown dwarfs for which nuclear burning is inhibited (see Richer et al. 2006).The determination of the current maximum stellar mass was more involved.We fixed the absolute maximum stellar mass in the Milky Way to ∼150 M as inferred from former analyses of the Arches cluster (Figer 2005;Kroupa 2005).More recent works investigating the cluster R 136 inferred a maximum stellar mass as high as ∼300 M (Crowther et al. 2010).However, we here adopted the conservative value of 150 M , also because higher values do not change our results significantly given the very low number of young and massive star clusters in the sample that might in principle host very massive stars.However, it is reasonable to assume that the maximum stellar mass depends on the total mass of the parent cluster.Even though this connection is still a matter of debate, we relied on the relation found by Weidner & Kroupa (2006), where the maximum stellar mass roughly corresponds to ∼10% of the cluster mass (see also Yan et al. 2023 for a more recent study of the correlation between the maximum stellar mass and the parent cluster mass).
Stellar evolution affects the high-mass tail of the distribution, which becomes progressively depleted because of SN explosions.In other words, the older the cluster, the lower the mass of the currently most massive member star (Buzzoni 2002).For the sake of completeness, we mention that the number of massive stars may also be reduced because of dynamical effects that may result in the ejection of stars from the cluster.(Oh et al. 2015) have estimated that clusters more massive than ∼10 3 M may eject a fraction of O-type stars up to ∼20% with a velocity 10 km s −1 .We neglected this complication, which may lower the final estimate of the cluster kinetic luminosity by a few tens of percent.For the purposes of our computations, we constrained M max to be the lower of these values.Further details about the specific values obtained per cluster are given in Sect.3.Because the mass estimation procedure defined here strongly relies on the capabilities of a single star identification, our result should rather be interpreted as a lower bound to the actual cluster mass, and hence mass-loss rate and wind power, as further discussed in the following sections.

Estimating cluster masses
The details of the steps undertaken in the mass estimation procedure are explained here.
(i) First, we evaluated individual cluster masses by assuming a broken power law functional form with four components for the stellar mass function ξ(M).Following Weidner & Kroupa ( 2004), A118, page 3 of 13 Celli, S., et al.: A&A, 686, A118 (2024) we set and M 1 = 1.00 M , M being the mass of the Sun.
(ii) Next, we applied to each cluster a kernel fit procedure to the G-band magnitude distribution of all its stellar members.The fit was applied to extract both the minimum M G min and the maximum M G max magnitudes of stars detected per cluster.In place of the maximum magnitude, we considered the mode of the stellar magnitude distribution per cluster, as it corresponds to the magnitude value up to which the Gaia survey data can be considered complete, as detailed in Appendix A.
Then, we derived the corresponding bolometric magnitudes M bol via the bolometric correction factor BC G .For a stellar Gband magnitude M G , this conversion is where BC G is obtained considering the observed effective temperature T eff of the stars within the cluster with the parameters α i given by Gaia FLAME 1 .
(iii) We converted the apparent bolometric magnitudes into intrinsic magnitudes by extracting line-of-sight extinction values A G from the v6 Gaia Starhorse catalog 2 .From these, the stellar luminosities were computed as (iv) These luminosities were used to derive the corresponding stellar mass M s (L s ) by inverting the following massluminosity relation, in the form of a broken power law, This relation was derived in Menchiari ( 2023) from merging existing parameterizations from Yungelson et al. (2008) and Eker et al. (2018), and it can be applied to the MS evolution in a broad mass range (extending even beyond 100 M ).In Eq. ( 5) 15, and = 0.817.We further assumed that all member stars were located at the same distance as the cluster: Hence, the minimum magnitude member star per cluster allowed us to infer the maximum mass M * max , while the mode magnitude determines the minimum stellar mass M * min .The distributions of these values for the clusters we analyzed are also provided in Appendix A.
1 https://gea.esac.esa.int/archive/documentation/GDR2/Data_analysis/chap_cu8par/sec_cu8par_process/ssec_ cu8par_process_flame.html 2 https://doi.org/10.17876/gaia/dr.2/52 The number of member stars N * between the minimum inferred stellar mass M * min and the maximum inferred stellar mass M * max is known from observations and shown in Fig. 1c.The distribution of the number of stars per cluster has a median value of 58, which is in line with typical values of open clusters.Finally, the normalization k of Eq. ( 1) was obtained by inverting the following expression: After normalizing the mass distribution of each cluster to match the Gaia number of observed stars, we obtained the total number of stars expected per cluster according to theoretical considerations, that is, we extended the following integral to compute the cluster mass in between a minimum M min and maximum M max stellar mass, in accordance with both stellar and cluster evolutionary theory.In particular, M min was assumed to be an absolute minimum mass of ∼0.08 M , related to the formation of brown dwarfs, while M max was determined by the following considerations.On one hand, the age of the cluster affects the highest value of the stellar mass that has not yet exploded into SN (Buzzoni 2002), under the hypothesis that the stellar formation process is instantaneous, such that the cluster age corresponds to the stellar age.In order to derive the relation of the maximum stellar mass allowed by its age, we adopted data from Fig. 3 of Buzzoni (2002), collecting the turn-over time of a sample of MS stars with solar metallicity (representing the reference time for SN explosion) as a function of their mass.We parameterized and inverted the analytical relation describing the data of Buzzoni (2002) according to the following expression: where the parameter values resulting from the fit are provided in Appendix A.
On the other hand, the maximum mass of a member star is determined in the first instance by the parent cluster mass, namely since the time of its formation (see Fig. 1 in Weidner & Kroupa 2004), most likely as a result of the interplay between the stellar feedback and the binding energy of the cluster-forming molecular cloud core.The latter effect naturally embeds the highest stellar mass observed.However, in order to determine the maximum stellar mass expected per cluster with the Weidner & Kroupa (2006) approach, knowledge of the cluster mass is required, which (by definition) is unknown.We hence proceeded iteratively by first computing a zeroth-order estimate of the cluster mass M c,0 , which we call seed cluster mass, which is obtained through Eq. ( 7) taking M max as the minimum between 150 M and the highest stellar mass allowed by stellar evolution given the cluster age T age (according to Eq. ( 8)).We hence proceeded by iterations aiming at refining the prediction about M max by accounting for the relation with the parental cluster mass, as derived by Weidner & Kroupa (2004).
Namely, for the clusters where M max (T age ) > M max (M c,0 ), we implemented an iterative approach based on the seed cluster mass values: We first derived the maximum stellar mass expected at the time of cluster formation according to Weidner & Kroupa (2006), and then adopted this value to recompute an updated cluster mass through Eq. ( 7), repeating the calculation to obtain A118, page 4 of 13 Fig. 2. Maximum expected stellar mass as a function of cluster age.The gray markers refer to the sample of selected young clusters from the Gaia catalog, whose maximum masses are evaluated from the iterative procedure described in the text.The solid blue line shows the result from Eq. ( 8), and the dashed red line represents the highest stellar mass value observed in the Milky Way at 150 M .the corresponding M max and so on, until after a few iterations, the procedure converged with an accuracy below 5% with respect to the cluster masses evaluated in the previous step.The values of M max obtained with this procedure are shown in Fig. 2 as filled markers, in comparison with the age-limited model from Buzzoni (2002).In Fig. 3a, we compare these values to those we inferred from Gaia observations (M * max ).The latter are clearly always lower than what is expected from theory, except for a few cases in which the adopted cluster ages are doubtful, as we discuss further in the next section.We additionally compared the maximum stellar mass values determined for each cluster with the mass values of observed star members that are available from the StarHorse routine.We found systematically lower values in the latter case, as expected.
The resulting cluster mass distribution is shown in Fig. 3b, with a median value of ∼413 M and a standard deviation of ∼1472 M .The highest value estimated in the sample is M max c 2.2 × 10 4 M and corresponds to Westerlund 1.We note that the distribution of the logarithm of young cluster masses has a belllike shape, which is expected to be induced by the completeness limit of the survey.The distribution can be fit with a Gaussian function with a mean value of log 10 [M µ /M ] = 2.6 and a standard deviation of σ log[M µ /M ] = 0.45.In the next section, we compare our mass estimates with other results in the literature, while in Sect.5, we describe how the mass uncertainty of our method was quantified based on the observed number of WR stars per cluster.

Comparison with existing cluster mass calculations
In this section, we compare our estimate of cluster masses with recent results from other works.In particular, Just et al. (2023) explored the latest MWSC data release that provided the socalled tidal radius r t , a spatial parameter first introduced by King (1962) to empirically describe the radial profile of the projected surface density of a star cluster, specifically the point where the projected density drops to zero.The set of three King's parameters was later found to fully describe the density profile of a spherical system in a quasi-equilibrium configuration.The tidal radius is therefore commonly adopted to compute the mass of virialized systems, such as the long-lived globular clusters.The application to open clusters is in turn more uncertain because their dynamics is dominated by the tidal forces of the Galaxy rather than by their internal potential.Nonetheless, Just et al. (2023) have computed tidal masses of 2227 clusters, regardless of their nature.Of these, 505 are younger than 30 Myr, but only 149 star clusters are in common with the Gaia sample we investigated.By comparing the masses of the common clusters, we found that the tidal approach yields masses that are systematically lower than those based on star counts used here.This behavior might be due to mass segregation effects (Portegies Zwart et al. 2010;Kirk & Myers 2011), inducing smaller estimates of the cluster radius that might strongly bias the results toward lower masses (Ernst et al. 2010).
A more proper comparison can be performed with a recent analysis by Almeida et al. (2023) that included a cluster mass evaluation for 773 open clusters from Gaia DR2.The authors developed a novel strategy based on the simulation of a synthetic cluster population followed by an isochrone fitting procedure that enabled them to constrain cluster ages, distances, masses, and binary fractions.This time, the common cluster sample (younger than 30 Myr) contains 102 objects, for which we investigated the differences not only in mass, but also in the other parameters derived by the authors.With respect to the parameters provided by the Gaia/MWSC catalogs, the distances found by Almeida et al. (2023)  The masses of single clusters agree overall within the uncertainties, but the average mass estimated by Almeida et al. (2023) is higher by ∼40% than ours, as shown in Fig. 4.This is expected because the two methods often rely on different cluster physical parameters to determine their masses.In particular, the method developed by Almeida et al. (2023) also accounts for binary stars, which we neglected.However, our sample is larger than that used by Almeida et al. (2023), and even more importantly, it includes stellar clusters with higher masses.This is crucial for the computation of gamma-ray emission from these objects, providing us with a more realistic estimate of their role as particle accelerators at the highest energies.
Finally, in Table 1 we compare our results for a few of the most massive star clusters of our sample with the mass and age values found in other works.The cited works usually estimated the cluster masses by evaluating their stellar content with a joint analysis of optical and near-infrared imaging.Our estimates are in general lower (by a factor ∼2 in three out of seven cases and a factor ∼15 in one case), with the exception of one case (Danks 2).We show in the next section that the number of WR stars also suggests that we underestimate the masses of the most massive stellar clusters by a factor ∼3 at least, as demonstrated below.

Estimating the cluster mass uncertainty using WR stars
The cluster masses estimated in Sect. 3 should be taken as a lower limit for two reasons.First, we neglected the stellar com- Notes.The cluster ages are also reported, namely the values we considered from Gaia/MWSC catalog (third column) and the values usually adopted in the literature (fifth column).
References ponent in binary systems, which may cause us to underestimate the cluster mass by a few tens of percent.For instance, the recent estimate by Hunt & Reffert (2024, see also Hunt 2023) suggested that the fraction of unresolved binaries can increase the cluster mass by 10%−30%.Second, the identification of cluster members by Cantat-Gaudin et al. ( 2020) was achieved through a specific neural network algorithm.However, techniques based on machine learning strongly depend on the initial set used to train the neural network.Other works that were based on different neural network architectures and had different training sets identified more stars (even up to a factor of 2; see, e.g., Mahmudunnobe et al. 2021;Van Groeningen et al. 2023), which clearly would impact any mass estimation method based on the star count.Moreover, Buckner et al. (2024) applied single star identification criteria to mock clusters, showing that for distances 1 kpc, they missed up a fraction >50% of the cluster members.Therefore, it would be highly desirable to have an independent estimate of the cluster masses.This can be obtained by comparing the number of WR stars associated with star clusters with those expected from theory, as explained below.
To provide a theoretical estimate for the number of WR stars, we assumed that at the end of the MS phase, all stars with a mass above M WR,min experience the WR phase.This minimum mass was estimated from stellar evolution codes to be between 22 and 37 M (Eldridge & Vink 2006).The average duration of the WR phase was estimated from stellar evolution models to be ≈[0.1−1]Myr, depending on several parameters such as the initial mass, metallicity, and stellar rotation (Meynet et al. 1994;Meynet & Maeder 2005).Because of the large theoretical uncertainties, we adopted the range inferred from the observed WR stars, that is, ∆t WR = [0.25−0.40]Myr (Rosslowe & Crowther 2015).Hence, under the assumption that the star formation process in the cluster is instantaneous and that the duration of the MS phase is regulated by Eq. ( 8), the number of WR stars is given by all stars whose MS phase lasted between t − ∆t WR and t (t being the cluster age).This translates into a mass interval between M WR,min (t) = min[M WR,min , M(t MS = t)] and M WR,max (t) = M(t MS = t − ∆t WR ).In reality, the star formation (SF) process is not instantaneous, but lasts for a time ∆T SF that mainly depends on the structure of the parent molecular cloud from which the cluster is born.Numerical simulations showed that this time can last up to ∼3 Myr (Krause et al. 2020) and that A118, page 6 of 13 the star formation rate is not constant during this time.Here, we accounted for the delayed SF process, while assuming that the star formation rate (SFR) remained constant during the entire period.Under these assumptions, the total number of WR stars as a function of cluster mass and age is where we explicitly introduce the dependence of ξ on the total cluster mass M c . Figure 5 shows the result of solving Eq. ( 9) for different cluster masses, fixing M WR,min = 25 M and using two different values for ∆T SF , namely 0 (instantaneous star formation), and 1 Myr.In the former case, N WR (t) shows a peak close to the threshold age when the most massive star in the cluster enters the WR phase, while in the latter case, this peak is smoothed out.The cutoff at t 6.5 Myr corresponds to the MS duration of stars with a mass of 25 M .
At the time of writing, the Galactic WR catalog (Rosslowe & Crowther 2015) 3 counts 669 WRs, 49 of which are assumed to be associated with 14 of the young clusters in our sample, namely Ruprecht 44, Westerlund 1 and 2 (Wd 1 and Wd 2), Trumpler 16, Sher 1, NGC 3603, Hogg 15, Danks 1 and 2, NGC 6231, Dolidze 3, Markarian 50, and Berkley 86 and 87.In descending order, Wd 1 is associated with 24 WRs (second only to the Galactic center), followed by Danks 1 with 6 WR stars, NGC 3603 with 5 WRs, and Wd 2 with 2 WRs.The remaining clusters are associated with a single WR.
We compared the predicted number of WRs by applying Eq. ( 9) to the cluster sample using the lower limits of the masses estimated in Sect. 3 and the associated ages.This procedure resulted in a predicted number of WRs equal to zero, regardless of the value of ∆T SF .This result may be due to both an underestimated value of the cluster masses and/or an incorrect estimate of their ages.The latter, in particular, is supported by several pieces of evidence, including Fig. 5 and the discussion in the previous section.As shown in Fig. 5, the clusters would be expected to host WR stars in the age range between ∼2.5 and 6.5 Myr, depending on their total mass.By assuming the catalog ages, we should not observe WRs in some cases in which they are reported.The only solution to this problem is a revision of their age estimate.However, an incorrect age estimate cannot entirely justify the discrepancy.Even when we assume that all clusters have an age of ∼4 Myr (to maximize the number of WR stars), the theoretical number of WRs expected in our cluster sample would be 11 (3 of which lie in Wd 1) compared to 49 (24 of which lie in Wd 1).We interpret this result as further evidence that the cluster masses evaluated in Sect. 3 are underestimated.
While confirming that the total cluster masses are underestimated, the observed number of WR stars provides us with a method for quantifying the mass uncertainty.By assuming that the actual cluster masses are a factor χ higher than the value estimated in Sect.3, we aim at reproducing the observed WR number.However, before proceeding with a quantitative estimate of χ, we note that neither the star cluster catalog nor the WR catalog are complete.We evaluated the completeness of the former in Sect. 2 to be about 2−3 kpc.The WR catalog can in turn be considered complete up to magnitude 8 in K s band (Rosslowe & Crowther 2015, see Sect. 4.3.1 and references therein).This roughly corresponds to a distance of 3 kpc in the Galactic plane (see, e.g., Kanarek et al. 2015).To be conservative, we restricted our analysis to WR/cluster associations closer 3 See the online catalog at http://pacrowther.staff.shef.ac.uk/WRcat/, last update from June 2023.6.Comparison of the observed number of WR stars in our cluster sample and the theoretical predictions where the cluster mass is multiplied by an arbitrary factor.The different line styles refer to different values of the duration of the star formation process, as indicated in the legend, and for three different cases: A (green), B (orange), and C (blue).The horizontal line corresponds to the number of observed WRs in the cluster sample, restricted to systems closer than 2.5 kpc, where both the Gaia and WR catalogs are assumed to be complete.See text for details.than 2.5 kpc.Within this distance range there are five associations with a total of seven WRs.
Figure 6 reports the predicted number of WR stars in clusters within 2.5 kpc as a function of χ and for different values of ∆T SF .To evaluate the impact of different assumptions, we show three sets of curves corresponding to different assumptions, -Case A: ∆t WR = 0.25 Myr and M WR,min = 37 M ; -Case B: ∆t WR = 0.40 Myr and M WR,min = 22 M ; -Case C: as case B, but with all clusters with estimated ages younger than 10 Myr set to be equal to 4 Myr.Case C was included to account for the uncertainty in the age estimate for very young clusters to maximize the number of WRs with respect to the cluster age.It should therefore be considered an extreme case.The result for case A suggests that our method underestimates the cluster masses by a factor 5, while for case B, we have χ 2.8−3.3,depending on the value of ∆T SF .
A118, page 7 of 13 Celli, S., et al.: A&A, 686, A118 (2024) The most optimistic of all cases, case C, requires χ 1.6−2.0.In addition to the total number of WR stars, we also compared the mocked distribution of clusters with at least one WR with the observed one by performing a Kolmogorov-Smirnov test.The estimated p-value is 0.4 for both cases A and C, and it is 0.6 for case B, suggesting that the latter is more compatible with the observed distribution than cases A or C. Finally, we note that when we neglected the completeness issue and applied the same analysis to the whole star cluster and WR samples, we obtained similar conclusions to those shown in Fig. 6.
A few additional comments are in order.The analysis discussed above may be affected by several additional uncertainties.From an observational point of view, WR stars might erroneously be associated with star clusters, as discussed, for example, in Rate et al. (2020).However, it seems unlikely that these misassociations apply to the majority of clusters, especially to those closer than 2.5 kpc where the parallax estimate is reliable.On the other hand, the theoretical approach we developed is quite simplified, especially concerning the physics of WR stars.The stars that undergo a WR phase and how long this phase lasts remain highly debated questions.However, the assumed intervals for M WR,min and ∆T WR cover the ranges usually adopted in the literature.A further increase in the number of the WR stars would require either M WR,min < 22 M or ∆T WR > 0.4 Myr.Another possibility to obtain a larger number of WRs is to change the IMF into a top-heavy distribution for high masses (i.e., harder than the Salpeter scaling M −2.3 assumed here), which would result in a higher number of massive stars.This may be the case for very massive star clusters such as those observed in the Galactic center, for instance, Arches (Hosek et al. 2019) and even Wd 1 (Lim et al. 2013), but not in general because the IMF of local clusters agrees well with a Salpeter distribution (Bastian et al. 2010).Most likely, the less well-constrained parameter that also affects our mass estimate is the value of cluster ages, which is difficult to obtain by isochrone fitting alone, particularly for young clusters.As we showed in Sect.4, the age values available in the literature differ significantly from the values quoted by Cantat-Gaudin et al. (2020) for at least some massive clusters.Nonetheless, even when this uncertainty was completely removed (as for case C), the cluster masses would still need to be increased by at least 60%.
Overall, we estimate an uncertainty for the cluster mass lower limits derived above of at least a factor ∼3 towards higher values.

A lower limit to cluster wind luminosities
Based on the above considerations, we proceeded to compute the cluster mass-loss rate by summing the mass-loss rates of all its member stars Ṁs (M) as For the stellar mass-loss rate, we relied on the analytical and empirical prescription given by Nieuwenhuijzen & de Jager (1990) relative to MS stars, depending on stellar luminosity L s (M) and radius R s (M), while the dependence of the stellar radius on mass was taken from Demircan & Kahraman (1991), resulting in The prescription adopted for the stellar mass-loss rate in Eq. ( 11) was chosen because its empirically driven nature fits the purposes of a statistical study, however it does not account for the possible clumpiness of stellar winds.Clumps can enhance the strength of the emission lines, which depends on the density squared (e.g., Hα lines), implying a lower value of the inferred mass-loss rate (Renzo et al. 2017).However, if the clumps are dense enough, they could become optically thick to the emitted radiation, suppressing the emission (see discussion in Vink 2022).The dominant effect remains unclear so far, and more investigation of the topic is required.For instance, in a combined analysis of the Hα and P v doublet lines detected from ζ Puppis, Oskinova et al. (2007) found that the two effects compensate almost perfectly and that the inferred mass-loss rate is very close to the rate predicted without accounting for clumpiness.Hence, in the following discussion, we simply rely on the empirical prescription for the mass-loss rate given in Eq. ( 11) and do not account for clumpiness.
The other important parameter is the wind speed, v w,s (M), whose value is related to the escape velocity from the star, G being the gravitational constant.With respect to v esc , the wind speed includes two additional effects.One is the presence of metals, which increases the effect of photon pressure (line-driven winds), and the second effect depends on the density of the wind itself, which acts by reducing the radiation pressure.These effects are taken into account by introducing a multiplication factor C and the Eddington luminosity of the star L Edd (M) = 4πGm p Mc/σ T , with the Thomson cross section σ T and proton mass m p , such that the stellar wind speed is given by (Kudritzki & Puls 2000) The coefficient C depends on the physical mechanism determining the wind pressure and is a function of the stellar effective temperature T eff , To avoid sharp discontinuities, we adopted a linear interpolation between C = 1 an C = 2.65.The latter value represents linedriven stellar winds, for example, where metal deexcitation contributes to the wind pressure such that the wind speed is expected to be higher than the escape speed at infinity.Finally, it is possible to derive the cluster wind luminosity L w,c through momentum conservation, A118, page 8 of 13 which, once inserted in the expression for the cluster wind luminosity, gives In addition to MS stars, we included the contribution of the observed WR stars to the mass-loss rate and wind luminosities of clusters accounting for the 49 observed WR stars associated with the clusters of our sample.Since the catalog also provides information concerning the spectral class of each star (WN, WC, WO, and their subclasses), we defined the values of the mass-loss rate and speed of their winds.Following Crowther (2007), we assumed an average mass-loss rate of 10 −4.9 M yr −1 , regardless of the spectral type, and wind speeds amounting to 1.6 × 10 8 cm s −1 and 2.3 × 10 8 cm s −1 for WN and WC/WO types, respectively.We assigned these values to the observed WRs associated with our cluster sample, obtaining wind luminosities of 9.3 × 10 36 erg s −1 and 2.3 × 10 37 erg s −1 for WN and WC/WO types, respectively, which we include in the luminosity computation of the clusters.As a result, we found that the most luminous open cluster of our sample is Westerlund 1.
The resulting cluster mass-loss rate and wind luminosity distributions are shown in Figs.7a and b.The median values are 1.3 × 10 −8 M yr −1 and 3 × 10 34 erg s −1 , but the high-energy tail extends up to ∼3 × 10 38 erg s −1 .For clarity, the panel with the wind luminosity shows histograms obtained both with and without the contribution of the observed WRs, indicating the crucial role of this stellar population in achieving higher energetics of the system.We additionally note that this result should be taken with the caveat that the star cluster and the WR samples are not complete at the same level.We may therefore have missed some WR stars associated with clusters in principle, especially for those located at distances 3 kpc.Based on the discussion about the cluster masses, we stress that the result for the wind luminosity should also be regarded as a lower limit.
As a final comment, we note that the above expressions for the wind speed and mass-loss rate assume solar metallicity.Young stellar clusters, however, presumably have a higher metal content than the Sun.It might therefore be wondered how the wind luminosity is affected.A proper description of the metallicity dependence goes beyond the aim of the present work because it is a nontrivial task to determine this.The metallicity has only been measured for a few star clusters, and this measure often only relied on one star or a few stars.We limit ourselves to comment on its possible impact.Both Ṁ and v w depend on metallicity.Vink & Sander (2021) found that massive O stars have a mass-loss rate ∝Z 0.42 , while the wind speed is v w ∝ Z 0.19 , implying that the final wind luminosity is L w ∝ Z 0.8 .For the less massive B-type stars, the dependence is far weaker.Considering that O-type stars dominate the wind power, we can state that if the average metallicity of YMSCs were ∼50% higher than in the Sun, the total luminosity would increase by ∼40%.
The full compilation of cluster masses, mass-loss rates, and wind luminosities is reported in Appendix C for further use.

Conclusions
Very high-energy and ultra high-energy gamma-ray observations of stellar clusters to date only enumerate a few bright systems in our Galaxy, which include Wd 1 and 2, the Cygnus Cocoon, and W43 (Aharonian et al. 2019;Cao et al. 2024).The nature of their energetic emission is still unclear, even though the hadronic scenario is preferred (Aharonian et al. 2019;Mestre et al. 2021).In a few cases, the reported emission requires protons and nuclei to be accelerated up to ∼PeV energies.
From a theoretical perspective, assessing the role of star clusters as extreme particle accelerators at a population level relies on the estimation of their masses and wind luminosities because they affect the energetics of non-thermal particles.Nonetheless, obtaining mass estimates of open clusters is a very difficult task, particularly for the young systems we are interested in, to probe the effectiveness of the WTS scenario.Young clusters are indeed more uncertain than older clusters, both in the use of stellar population synthesis models and observationally because velocity dispersion measurements are frequently lacking.Hence, finding the most reliable technique for estimating their mass remains an open issue in the astronomical community that is very much debated and often dependent on the system parameters.
We have developed a mass estimation method for star clusters based on the properties of their individual member stars.We applied it to the Gaia DR2 open cluster catalog and obtained cluster masses ranging from ∼10 to ∼10 4 solar masses, in agreement with the expectations for the disk population of the Milky Way.The highest mass value obtained for the systems in our sample is for Wd 1, with a lower limit of M Wd1 ≥ 2.2 × 10 4 M .A118, page 9 of 13 Celli, S., et al.: A&A, 686, A118 (2024) Literature studies with a completeness-corrected IMF reported its mass to be 4.91 +1.79  −1.49 × 10 4 M (Gennaro et al. 2011), which is higher by a factor ∼2 than our estimate.It is worth highlighting that strong evidence of mass segregation was found in this cluster, despite its young age of 4.0 ± 0.5 Myr (Gennaro et al. 2011).
We have further explored the usage of WR stars to quantify the uncertainty on the estimated masses.By developing a simplified model for WR occurrence in the star clusters of our Galaxy and comparing it to observations, we conclude that at least a factor 3 toward higher values should be considered for a realistic estimate of the cluster masses.This factor encloses several limitations of our procedure, which neglects binary systems as well as unresolved stars and cluster members that might belong to the cluster, but were not identified by the membership algorithms (Buckner et al. 2024).This possibility is further supported by the comparison with alternative mass estimation methods available in the literature, which are restricted to a smaller cluster sample and lower cluster mass range, however.From both the WR estimation method and the comparison with cluster masses available in the literature, evidence is given that some of the cluster ages claimed in the catalogs are also unsuitable to describe their properties.
The estimated masses were used to also calculate the total wind luminosity, which is the most important parameter for establishing the gamma-ray luminosity for young stellar clusters, and which will be used in two companion papers to evaluate the detection prospects of known clusters in the gamma-ray band.As discussed throughout the text, the method developed here for estimating cluster masses relies on the assumption of a stellar mass distribution within the clusters that is normalized to reproduce the number of stars observed by Gaia in each cluster.To do this, it is necessary to define the mass range of the detected stars as integration limits of Eq. ( 6).We inferred these values from available observations as described in Sec. 3, relying on the measured G-band stellar magnitudes, and converting them into bolometric magnitudes first, to then obtain stellar absolute luminosities after correction for the G-band extinction along the line of sight.We finally adopted the mass-luminosity relation in Eq. ( 5) relatively to the MS evolution to obtain stellar masses.The highest magnitude star per cluster hence determines the minimum stellar mass M * min , while the lowest magnitude star determines the maximum stellar mass observed M * max .In reality, the magnitude distribution of stars within clusters does not show a sharp drop at the highest magnitude, but it rather appears as    Buzzoni (2002).The functional form of the red line, obtained through a minimum residual technique, is provided in Eq. ( 8). a peaked distribution, whose mode is determined by the limiting magnitude of the observation in the given direction.Hence, we adopted the mode of the magnitude distribution rather that its highest value to infer the minimum stellar mass, shown in After fixing M * min and M * max , the normalization constant k of the cluster stellar mass function was obtained.We then determined the expected cluster masses and luminosities by extending the mass range of stars populating these systems, as suggested by theoretical arguments rather than being limited by instrument detection capabilities.To do this, we simply extended the integration mass range between 0.08 M and a maximum value, determined by either the cluster mass or its age, whichever yielded the lower value.
Older clusters are expected to host less massive stars than younger clusters because of the occurrence of SN explosions.The stellar lifetimes depend inversely on their mass.A 10 M star will explode in an SN after about 10 Myr, while a 100 M star will explode after only ∼1 Myr.We parameterized this A118, page 11 of 13 Celli, S., et al.: A&A, 686, A118 (2024) Weidner & Kroupa (2006).The gray markers refer to the selected young cluster sample.
relation through Eq. ( 8), which we fit to stellar evolution data from Buzzoni (2002)  For the effects of the parent cluster mass on the maximum stellar mass that might form within it, we followed the description from Weidner & Kroupa (2006), represented in Fig. A.5 as the solid black line.Because of the non-linearity intrinsic to the application of this constraint, which would require knowing the cluster mass to determine the maximum stellar mass hosted, we implemented an iterative technique that took as input the seed cluster mass, evaluated with a maximum stellar mass limited by either the cluster age (Buzzoni 2002) or by 150 M .From this seed mass, we computed the expected M max following Weidner & Kroupa (2006).For all clusters of our sample for which the M max previously evaluated through the age-limited maximum stellar mass was higher than the masses evaluated based on the seed cluster mass, we first corrected M max with the latter value.As a consequence, we computed the updated cluster mass value M c from Eq. ( 7), which would now correspond to a different maximum stellar mass according to Weidner & Kroupa (2006).We proceeded with these iterations until convergence was achieved on the cluster mass values.The final M max versus M c values are shown as dots in Fig. A.5 for the entire young cluster sample we analyzed.The markers closely following the theoretical line refer to the clusters for which the iterative procedure was effective.As visible, the parent cluster mass strongly constraints the maximum stellar mass that will be formed in lowmass clusters.

Appendix B: Energetic considerations about the contributions of SNe and stellar winds
The methods developed here for estimating the mass and wind luminosity of star clusters further allowed us to infer the temporal evolution of the stellar population of which individual stellar The obtained median values amount to 7 × 10 50 erg and 5 × 10 48 erg for SNe and winds, respectively.This illustrates the relevance of SNe in the energetics considerations of stellar clusters, and therefore also the relevance for particle acceleration and related gamma-ray emission.We note, however, that Eq. (B.2) is expected to be higher when the contribution of WRs is included.Nonetheless, this simplified estimate already suggests that the energy injected into winds throughout the cluster history is higher than 10 51 erg for some systems, namely Wd1, Danks1, Danks2, Patchick94, UBC344, and NGC3603.

Fig. 1 .
Fig. 1.Distributions of physical parameters regarding the selected cluster sample from Gaia DR2 (i.e.younger than 30 Myr).(a) Cluster ages in units of Myr.(b) Cluster distances in units of pc.(c) Number of observed member stars per cluster.

Fig. 3 .
Fig. 3.Estimated stellar and cluster masses from the Gaia DR2 young cluster sample.(a) Comparison of the maximum stellar mass per cluster obtained in our calculation (ordinate) and that inferred from the Gaia observation (abscissa).(b) Distribution of cluster masses in units of solar masses for the selected young clusters.

Fig. 4 .
Fig. 4. Comparison of the cluster masses estimated in this work (abscissa) and by Almeida et al. (2023, ordinate) for the common sample of 102 young clusters from Gaia DR2.The dashed red line shows the identity.

Fig. 5 .
Fig.5.Theoretical prediction for the number of WR stars, N WR , inside a cluster as a function of its age and for different cluster masses as labeled.The duration of the SF process is assumed to be 0 and 1 Myr for the solid and dashed lines, respectively.N WR is divided by M c /(10 5 M ).

Fig. 7 .
Fig. 7. Distribution of the cluster mass-loss rates (panel a) in units of solar masses per year and of the wind luminosity (panel b) in units of erg s −1 for the selected young clusters.The latter shows the cluster wind luminosity computed both with (blue) and without (red) the contribution of observed WRs.The shaded part of the histogram represents overlapping values.

Fig
Fig. A.2. Distribution of the maximum stellar masses as inferred from Gaia observations, described in Sec. 3, for the clusters in the selected sample.
Fig. A.3.Scatter plot between the minimum and maximum stellar mass values inferred from Gaia observations as described in Sec. 3 for the clusters in the selected sample.The dashed red line represents the bisector.

Fig
Fig. A.4. Mass-age relation for main-sequence stars with data fromBuzzoni (2002).The functional form of the red line, obtained through a minimum residual technique, is provided in Eq. (8).
Fig. A.1.The maximum stellar mass per cluster is in turn represented in Fig. A.2, its highest value amounting to ∼18 M .The scatter plot among the two is shown in Fig. A.3.
Fig. A.5. Maximum stellar masses M max as a function of the seed cluster masses, shown on the abscissa.The solid line represents the numerical solution we derived followingWeidner & Kroupa (2006).The gray markers refer to the selected young cluster sample.
in the mass range −0.2 ≤ log 10 [M/M ] ≤ 2.1, as shown in Fig. A.4.The result is shown in Fig. A.4, where the red curve was obtained for the following set of values: A = 4.35 ± 0.07, B = 1.30 ± 0.04, and C = 1.101 ± 0.035, with a χ 2 of 0.5.
Fig. B.1.Comparison of the expected energetics injected by cluster winds (blue histogram) and SN explosions (red histogram) in the selected sample of young stellar clusters.
Almeida et al. (2023)r distances smaller than 2 kpc.In contrast, for farther clustersAlmeida et al. (2023)obtained significantly smaller distances than Gaia/MWSC catalogs.More interestingly, the revised age estimates largely disagree with the values provided by Gaia/MWSC given that the age distribution in theAlmeida et al. (2023)sample is quite narrow and peaked around (7.2 ± 0.30) Myr.These results indicate that the age values provided by the Gaia/MWSC catalogs might be underestimated in the case of younger clusters (<6 Myr), but overestimated for the older ones.This behavior is also supported by the investigation of the observed versus expected WR stars in clusters that we present in Sect. 5.

Table 1 .
Comparison of our estimates of mass (second column) and the values provided by other works (fourth column) for a subset of clusters.