We present a nearinfrared NIR candidate star cluster catalog for the central kiloparsec of M82 based on new JWST NIRCam images. We identify star cluster candidates using the F250M filter finding 1357 star cluster candidates with stellar masses gt10SUP4SUP M SUBSUB. Compared to previous optical catalogs nearly all 87 of the candidates we identify are new. The star cluster candidates have a median intrinsic cluster radius of 1 pc and stellar masses up to 10SUP6SUP M SUBSUB. By comparing the colorcolor diagram to dustfree YGGDRASIL stellar population models we estimate that the star cluster candidates have A SUB V SUB  324 mag corresponding to A SUB2.5mSUB  0.32.1 mag. There is still appreciable dust extinction toward these clusters into the NIR. We measure the stellar masses of the star cluster candidates assuming ages of 0 and 8 Myr. The slope of the resulting cluster mass function is   1.9  0.2 in excellent agreement with studies of star clusters in other galaxies.
['Star clusters', 'Young massive clusters', 'Young star clusters', 'Starburst galaxies', 'Infrared galaxies', '1567', '2049', '1833', '1570', '790', 'Astrophysics - Astrophysics of Galaxies']
JWST Observations of Starbursts Massive Star Clusters in the Central Starburst of M82
{'JWST Observations of Starbursts: Massive Star Clusters in the Central Starburst of M 82': "<!-- image -->  \n1 Steward Observatory, University of Arizona, Tucson, AZ 85721, USA  \n2 Department of Astronomy, University of Maryland, College Park, MD 20742, USA  \n- 3 Joint Space-Science Institute, University of Maryland, College Park, MD 20742, USA  \nInstituto Nacional de Astrof'ısica,  \n'  \nOptica y Electr'onica, Luis Enrique Erro 1, Tonantzintla, 72840 Puebla, Mexico  \n5 Benem'erita Universidad Aut'onoma de Puebla, Av. San Manuel, 72000 Puebla, Mexico  \n- 6 Space Telescope Science Institute, 3700 San Martin Drive, Baltimore, MD 21218, USA  \n7 Max Planck Institut for Astronomy, Konigstuhl 17, 69117 Heidelberg, Germany  \n8 European Space Agency, c/o STScI, 3700 San Martin Drive, Baltimore, MD 21218, USA  \n9 Department of Physics and Astronomy, University of Wyoming, Laramie, WY 82071, USA  \n10 Department of Physics and Astronomy, University of Kansas, 1251 Wescoe Hall Drive, Lawrence, KS 66045, USA  \n- 11 Coll'ege de France, 11 Pl. Marcelin Berthelot, 75231 Paris, France  \n12 Observatoire de Paris, 61 avenue de l'Observatoire, 75014 Paris, France  \nCentre for Astrophysics and Supercomputing, Swinburne University of Technology, Hawthorn, VIC 3122, Australia  \n13  \n14 ARC Centre of Excellence for All Sky Astrophysics in 3 Dimensions (ASTRO 3D)  \n15 Universitat Heidelberg, Zentrum fur Astronomie, Institut fur Theoretische Astrophysik, Albert-Ueberle-Str. 2, D-69120 Heidelberg,  \nGermany  \n16 Departamento de Astronom'ıa, Universidad de Concepci'on, Barrio Universitario, Concepci'on, Chile  \n17 Observatorio Astron'omico Nacional (IGN), C/Alfonso XII, 3, E-28014 Madrid, Spain  \n18 Centro de Desarrollos Tecnol'ogicos, Observatorio de Yebes (IGN), 19141 Yebes, Guadalajara, Spain  \n19 Universitat Heidelberg, Interdisziplinares Zentrum fur Wissenschaftliches Rechnen, Im Neuenheimer Feld 205, D-69120 Heidelberg, Germany  \n20  \nStratospheric Observatory for Infrared Astronomy, NASA Ames Research Center, Mail Stop 204-14, Moffett Field, CA 94035, USA  \n21  \nJet Propulsion Laboratory, California Institute of Technology, 4800 Oak Grove Dr., Pasadena, CA 91109, USA  \n- 22 Department of Astronomy, The Ohio State University, Columbus, OH 43210, USA\n- 23 Sterrenkundig Observatorium, Ghent University, Krijgslaan 281 - S9, B-9000 Gent, Belgium\n- 24 New Mexico Institute of Mining and Technology, 801 Leroy Place, Socorro, NM 87801, USA\n- 25 National Radio Astronomy Observatory, P.O. Box O, 1003 Lopezville Road, Socorro, NM 87801, USA\n- 26 Dept. F'ısica Te'orica y del Cosmos, Universidad de Granada, 18071, Granada, Spain\n- 27 Leiden Observatory, Leiden University, P.O. Box 9513, 2300 RA Leiden, The Netherlands", 'ABSTRACT': 'We present a near infrared (NIR) candidate star cluster catalog for the central kiloparsec of M 82 based on new JWST NIRCam images. We identify star cluster candidates using the F250M filter, finding 1357 star cluster candidates with stellar masses > 10 4 M ⊙ . Compared to previous optical catalogs, nearly all (87%) of the candidates we identify are new. The star cluster candidates have a median intrinsic cluster radius of ≈ 1 pc and have stellar masses up to 10 6 M ⊙ . By comparing the color-color diagram to dust-free yggdrasil stellar population models, we estimate that the star cluster candidates have A V ∼ 3 -24 mag, corresponding to A 2 . 5 µ m ∼ 0 . 3 -2 . 1 mag. There is still appreciable dust extinction towards these clusters into the NIR. We measure the stellar masses of the  \n4  \nstar cluster candidates, assuming ages of 0 and 8 Myr. The slope of the resulting cluster mass function is β = 1 . 9 ± 0 . 2, in excellent agreement with studies of star clusters in other galaxies.  \nKeywords: Star clusters (1567), Young massive clusters (2049), Young star clusters (1833), Starburst galaxies (1570), Infrared galaxies (790)', '1. INTRODUCTION': "The cycle of star formation is a key driver of the evolution of galaxies. At high surface densities of star formation, the majority of stars form in gravitationally bound clusters (e.g., Kruijssen 2012; Krause et al. 2020). Various phases of star cluster formation inject energy and momentum (i.e., feedback) into the surrounding interstellar medium (ISM). Recently, it has become clear that pre-supernovae feedback can be a dominant mechanism to clear the natal gas from a star cluster (e.g., Chevance et al. 2022; Farias et al. 2024). However, the earliest phases of star cluster formation and evolution are enshrouded in dust and gas, making them difficult to observe. Multiwavelength observations of star clusters are therefore necessary to probe the full cycle of star cluster evolution.  \nIn regions of high molecular gas surface densities, the star formation rates (SFRs) are higher than expected for a simple scaling of the Kennicutt-Schmidt relation in more typical conditions (e.g., Kennicutt & Evans 2012). This suggests that the process of star formation in starburst regions is physically different than other regimes. One consequence, seen in both observations and models, is that the fraction of stars formed in clusters increases with the surface density of star formation (e.g., Goddard et al. 2010; Adamo et al. 2011; Silva-Villa & Larsen 2011; Kruijssen 2012; Leroy et al. 2018, though see Cook et al. 2023). Moreover, this high activity mode tends to form 'super' star clusters, massive (M ∗ ≳ 10 5 M ⊙ ) and compact ( r ≲ 1 pc) bound groups of stars (e.g., Portegies Zwart et al. 2010).  \nThe small, irregular galaxy M 82 (NGC 3034; M dyn ≈ 10 10 M ⊙ , R dyn ≈ 4 kpc; Greco et al. 2012; Dale et al. 2023) is the archetypal starburst in the local Universe. Located in the M81 group at a distance of 3.6 Mpc (Freedman et al. 1994), the starburst in M 82 is driven by the tidal interactions in the group (e.g., Yun et al. 1994; Forster Schreiber et al. 2003). The starburst region in M 82 roughly extends over the central kiloparsec. M82 has experienced two recent bursts of star formation (Forster Schreiber et al. 2003). The first occurred between 8 -15 Myr ago, peaked at a SFR of 160 M ⊙ yr -1 ,  \nand was driven primarily by tidal interactions. The second occurred 4 -6 Myr ago, peaked at a SFR of 40 M ⊙ yr -1 , and was driven by a bar inflow (the formation of which likely resulted from the interaction, e.g., Forster Schreiber et al. 2003). M82 has a 'quiescent' total infrared luminosity-derived SFR ≃ 12 M ⊙ yr -1 (Herrera-Camus et al. 2018a,b), which traces the SFR over the last 100 Myr (e.g., Calzetti 2013). A result of the starburst phases in the recent past is the massive multiphase outflow of material from the central starburst region (see e.g., Bolatto et al. 2024, and references therein) - a defining feature of M 82.  \nA number of massive star clusters have been identified and characterized in the central starburst of M 82 using several facilities, including with the Hubble Space Telescope (HST) using optical (O'Connell et al. 1995; Melo et al. 2005; Mayya et al. 2008) and NIR (McCrady et al. 2003) images, optical-infrared spectroscopy (McCrady et al. 2003, 2005; McCrady & Graham 2007; Westmoquette et al. 2014), mid-infrared (MIR) imaging (Gandhi et al. 2011), and Submillimeter Array (SMA) observations (M. J. Jim'enez-Donaire et al. in prep.). The HSTbased cluster catalog of Mayya et al. (2008) is the most complete accounting of the optical star clusters in M 82. Using HST F435W, F555W, and F814W images, they identified 260 star clusters within the central 900 pc and 363 star clusters outside of the central starburst. The pixel scale of these images is 0.05 '' (0.88 pc). In order to reject compact 'starlike' objects from their catalog, Mayya et al. (2008) require that the candidates have a minimum (maximum) diameter of 3 (30) pixels and an area of > 50 pixels. Therefore, the clusters in their catalog have Gaussian radii > 1 . 3 pc and areas of ≳ 40 pc 2 . Their clusters tend to have A V ≲ 6, and their cluster mass functions have a power-law slope of 1.8 for clusters in the starburst, assuming an age of 8 Myr.  \nThe central region has a patchy distribution of dust with extinction in excess of A V = 30 mag for the NIRdetected star clusters (Forster Schreiber et al. 2003). As young star clusters are compact (radii ≲ 1 pc; e.g., Portegies Zwart et al. 2010; Krumholz et al. 2019; Brown &Gnedin 2021; Levy et al. 2022) and tend to resolve into smaller sub-structures (e.g., Levy et al. 2021, 2022), high spatial resolution observations are needed to obtain an accurate census of star cluster populations. JWST offers  \nparsec-scale spatial resolutions for nearby galaxies such as M82 at unprecedented sensitivities, which provides a unique opportunity to identify and study clusters over a range of evolutionary stages from embedded to more exposed. The reduced extinction at NIR wavelengths as compared to traditional optical wavelengths allows us to obtain a more complete catalog of all star clusters. Both NIRCam and MIRI onboard JWST have already proven to be powerful tools to detect and characterize star clusters in nearby galaxies (e.g., Linden et al. 2023; Rodr'ıguez et al. 2023; Schinnerer et al. 2023; Sun et al. 2024).  \nIn this Letter, we construct and present the JWST NIRCam star cluster catalog in the central region of the M82 starburst. We briefly present the observations and data processing in Section 2. In Section 3, we describe the NIRCam star cluster catalog. We discuss the overlaps with star cluster catalogs at other wavelengths, color-color diagrams, and the cluster mass function in Section 4. We summarize our results in Section 5.", '2. OBSERVATIONS AND DATA REDUCTION': "We identify clusters using NIRCam data obtained as part of JWST Cycle 1 GO program #1701 (P.I. Alberto Bolatto). We direct the reader to Bolatto et al. (2024) for full details of the observations and data reduction. Briefly, we use NIRCam (Rieke et al. 2023) SUB640 observations of the central 50 '' ( ≈ 880 pc) of M82 1 . A three color image of the continuum filters from Program #1701 (F140M, F250M, and F360M) is shown in Figure 1, which was produced with multicolorfits (Cigan 2019). The left panel of each inset pair shows the same filters over four cluster-rich 50 pc × 50 pc regions.  \nOur JWST program observed in three 'continuum' filters: F140M, F250M, and F360M to enable a broad range of science goals (see Bolatto et al. 2024 for more details). While the F140M filter has higher angular resolution among these three filters, it suffers the most from dust extinction meaning some clusters will be missed (e.g., Figure 1 and Figure 2 of Bolatto et al. 2024). F360M, on the other hand, has the lowest angular resolution and the continuum is contaminated by the 3.3 µ m polycyclic aromatic hydrocarbon (PAH) feature as well as a water-ice absorption feature, which may be particularly prominent in the young star clusters (see Sandstrom et al. 2023 and Bolatto et al. 2024 for a more in-depth discussion). Based on earlier spectroscopy of M82 (Sturm et al. 2000; Forster Schreiber et al. 2001), we expect the F250M filter to be largely free from line  \nemission and to provide a clean estimate of the stellar continuum. Therefore, F250M is the ideal filter in our suite of NIRCam observations with which to identify NIR-bright star clusters. The F250M filter has a point-spread-function (PSF) full-width-half-maximum (FWHM) of 0.085 '' (1.5 pc) and a pixel scale of 0.042 '' (0.7 pc).", '2.1. Ancillary Data': "For comparison, in Figure 1 we show data from HST in the right panel of each zoom-in pair. The selected filters - F435W, F555W, F814W - were downloaded from MAST 2 and are the same as those used by Mayya et al. (2008) to catalog the star clusters across M 82. These images have pixel scales of 0.05 '' . It is known that some HST images have astrometric offsets ≲ 1 '' with respect to e.g., interferometric images (e.g., Leroy et al. 2018; Levy et al. 2021). Because Bolatto et al. (2024) found that these JWST NIRCam images are well aligned to a radio catalog of supernova remnants, we assume that the absolute astrometry of the NIRCam images is correct. To align the HST images to the JWST NIRCam data, we measured offsets to shift the HST astrometry by ( α , δ ) = (-1.52 '' , -0.67 '' ) so that bright clusters visible in both images were cospatial. The light blue circles in the HST insets of Figure 1 show the cataloged star clusters in the central region from Mayya et al. (2008).", '3.1. Cluster Identification': 'Star cluster candidates were identified from the F250M image using SourceXtractor++ 3 (Bertin et al. 2022; Kummel et al. 2022). We started from the default configuration file, with the following changes. The detection threshold above the background ( detection-threshold ; i.e., the signal-tonoise ratio relative to the measured local background) was increased to 10 due to the presence of numerous bright point sources (see the insets in Figure 1). The minimum number of contiguous pixels above the detection threshold to qualify as a detection ( detection-minimum-area ) was increased to 9 pixels to remove definite point sources 4 from the initial catalog. We also enable source cleaning ( use-cleaning=1 ) which removes false detections near bright objects in the preliminary catalog (though, as discussed below, this  \nDecl. (J2000)  \nFigure 1. The central 870 pc of M 82 seen with NIRCam showing F140M (blue), F250M (green), and F360M (red) on an asinh scale. Insets show zoom-ins to 50 pc × 50 pc regions around some representative massive star cluster candidates. The circles indicate the massive clusters identified in the NIRCam image, where the radius of the circle is the fitted (deconvolved) radius (see Section 3 and Table 1). Blue (pink) circles show clusters which do (do not) overlap with the HST catalog (Mayya et al. 2008, see Section 4.1). Clusters 1 and 2 are labeled (see Sections 3.2 and 4.1). The right insets show HST three-color images at the same locations (blue: F435W, green: F555W, red: F814W). The light blue circles in the HST insets show the star clusters identified by Mayya et al. (2008) in these selected regions using the same HST data, where the radius of each circle is 5 pc. The pink and blue circles are the same as in the left panels.  \n<!-- image -->  \nstep removes very few sources). With these settings, SourceXtractor++ identified 2472 candidates.  \nSourceXtractor++ sometimes defines two clusters that significantly overlap. From visual inspection, we determine that such cases are better described as a single cluster. The fractional overlap of each candidate with every other candidate is measured using their center positions and a radius equal to the PSF FWHM. We flag those candidates with > 75% of their areas overlapping. We replace the coordinates of the overlapping clusters with the median of their coordinates determined by SourceXtractor++ . As described above, SourceXtractor++ tends to find false sources near bright objects. Although the source cleaning step should remove these, we found that, in practice, very few clusters were excluded during this step. To account for this, we remove candidates with slightly less overlap but where the weaker source does not appear to be cluster-like based on visual inspection. We determine the fractional overlap with a radius equal to twice the PSF FWHM, and we flag those candidates with > 50% of their areas overlapping. We extract the flux in a circular aperture of the same radius and keep only the brightest overlapping cluster. SourceXtractor++ also tends to pick up false sources along the edges of the image. We remove these from the preliminary catalog by discarding sources within 10 pixels of the image edge. After these initial cuts to remove spurious sources, we are left with 1767 star cluster and stellar candidates. We will discuss further refinements to this catalog in Section 3.4.', '3.2. Radial Profiles and Radii Measurements': "We measure the size of each star cluster candidate by constructing a radial profile. Our methodology closely follows that of Levy et al. (2022, see their Section 4.3). First we select a square region around each cluster that is 10 × the PSF FWHM on each side (0.85 '' ≈ 15 pc). We then mask out other cluster candidates in this region (based on the positions in the preliminary star cluster candidate catalog) using a circle with a diameter equal to twice the PSF FWHM. This mitigates the effects of other bright objects from affecting the radial profile of the desired source. We extract the radial profile for each candidate in circular annuli centered on the coordinates determined by SourceXtractor++ and with a width equal to 1 pixel. We take the median of the surface brightness in each ring, and the uncertainty is the standard deviation. An example for one cluster is shown in Figure 2.  \nWe then fit each radial profile with a Gaussian function plus a constant offset (e.g., Equation 1 of Levy et al. 2022). This enables us to model the cluster sur-  \nFigure 2. An example radial profile and Gaussian fit, shown here for Cluster 2 (marked in the upper right JWST inset in Figure 1). The black points show the radial profile extracted from the F250M image, where the error bars reflect the standard deviation in each ring. The blue curve and shaded region show the Gaussian fit to the radial profile and the uncertainties. The vertical blue dashed line indicates the HWHM of the fitted Gaussian. The gray vertical shaded region shows the HWHM of the F250M PSF. The red curve and shaded region show the deconvolved Gaussian fit, and the vertical red dashed line shows the HWHM of the deconvolved Gaussian (reported in Table 1).  \n<!-- image -->  \nface brightness profile as a Gaussian and to determine the local background level (with other bright sources removed), which we subtract from the radial profile and subsequent fits. To obtain uncertainties on the fitted Gaussian, we perform a Monte Carlo simulation over the surface brightness uncertainties. Our uncertainties come from the standard deviation of 500 trials (shown as the blue shaded region in Figure 2).  \nWith this Gaussian model, we deconvolve the PSF in a simplified way. For this first analysis, we assume that the JWST PSF is Gaussian 5 . Our estimate of the deconvolved (i.e., intrinsic) cluster radii is achieved by removing in quadrature the PSF half-width-halfmaximum (HWHM) from the cluster Gaussian fitted HWHM. Through flux conservation, we calculate the deconvolved Gaussian profile and uncertainties, shown in red in Figure 2. We report these values in Table 1.  \nBecause of our approximate method to remove the PSF from our radius measurements, the deconvolved radii reported in Table 1 should be treated as a first esti-  \nmate, not a robust measurement, of the intrinsic cluster radii. A more robust analysis is planned by our team (B. Cuevas-Otahola et al. in prep.) accounting for the actual PSF shape (e.g., WebbPSF ; Perrin et al. 2014, see also Rigby et al. 2023).", '3.3. Stellar Mass Estimates': 'We measure the spectral flux density of each cluster from the area under the Gaussian fit to the F250M radial profile before deconvolution: F = 2 πI peak r 2 HWHM , where F is the spectral flux density in mJy, I peak is the peak intensity of the Gaussian profile, and r HWHM is the cluster radius before PSF deconvolution 6 . We report these spetcral flux density values and propagated uncertainties in Table 1.  \nTo convert the measured spectral flux densities into a stellar mass via a mass-to-light ratio (Υ ∗ ), we must assume a cluster age, as the F250M data alone cannot constrain the cluster ages. For this first analysis, we take two age estimates as bounds to the true age of this cluster population. For a lower age limit, we assume a zero age main sequence (ZAMS) stellar population. This limit is applicable to the youngest, most embedded, cluster candidates in the burst. In the yggdrasil SSP models, all of the stars start on the main-sequence at age = 0 Myr 7 .  \nWe assume 8 Myr as our upper age limit, following the analysis of Mayya et al. (2008) which built on primarily optical results from Satyapal et al. (1997), Forster Schreiber et al. (2003), and Melo et al. (2005). Because the NIR is less susceptible to dust extinction, we are able to probe down to more embedded, likely younger, clusters than optical studies. Older, more evolved clusters are visible in both optical and NIR observations. Therefore, the mean age of clusters probed with NIR measurements will be slightly younger than those detected in optical studies. Indeed, Satyapal et al. (1997) find an average age of 6 Myr for their NIR-identified star clusters in M 82 compared to 8 Myr for the optically-i  \nentified clusters 8 . To be conservative, we take 8 Myr as a likely upper limit on the (representative) cluster age.  \nWe determine the Υ ∗ using the yggdrasil single stellar population (SSP) models (Zackrisson et al. 2011). This code (which has also been used to analyze NIRCam observations of star clusters in luminous infrared galaxies; Linden et al. 2023) computes the evolution of an instantaneous burst of star formation in a 10 6 M ⊙ cluster assuming a Kroupa (2001) initial mass function and Padova-AGB stellar evolution tracks (e.g., Bertelli et al. 2008, 2009). We adopt models at solar metallicity, since the metallicity in the central 500 pc of M 82 is roughly solar (Lopez et al. 2020). For simplicity, we assume maximal nebular emission ( f cov = 1). We determine the luminosity of the yggdrasil cluster in the F250M filter by convolving the spectral energy distribution (SED) from the yggdrasil model with the F250M filter throughput (version 5.0, November 2022) 9 , yielding the luminosity of a 10 6 M ⊙ cluster in the F250M filter for each assumed age, which we then convert to the equivalent specrtral flux density. The resulting Υ ∗ in the F250M filter are 0.35 M ⊙ /L ⊙ for 0 Myr and 0.18 M ⊙ /L ⊙ for 8 Myr. At the distance of M 82, these translate to 1 . 7 × 10 5 M ⊙ mJy -1 and 8 . 8 × 10 4 M ⊙ mJy -1 respectively. We multiply the measured spectral flux density (and uncertainties) of each cluster by Υ ∗ to obtain their stellar masses (M ∗ ) for each assumed age (Table 1). As we discuss in Appendix A, the ZAMS stellar masses correspond to the median for a < 10 Myr star cluster population (where the older starburst in M 82 occurred 10 Myr ago; Forster Schreiber et al. 2003). The minimum Υ ∗ occurs at ∼ 8 Myr, and hence those stellar masses reflect the minimum cluster masses.', '3.4. The Final Candidate Massive Star Cluster Catalog': "After the cluster identification, radius measurements, and stellar mass estimates, we implement the following cuts to build the final massive candidate star cluster catalog, presented in Table 1. We note the number of candidates removed in parentheses at the end of each step below.  \nTable 1. NIRCam Catalog of Massive Star Cluster Candidates  \nNote -Table 1 is published in its entirety in the machine-readable format. A portion is shown here for guidance regarding its form and content.  \n- 1. Red supergiants (RSGs), asymptotic giant branch (AGB) stars, and other massive stars emit strongly at 2.5 µ m and can be confused with massive star clusters based on their intensity alone (e.g., Bolatto et al. 2007; Levesque 2018; Boyer et al. 2024). In order to remove these objects from our star cluster catalog, we measure the concentration of the star cluster from the radial profile to determine if the source is extended or consistent with a point source. We remove cluster candidates with ≥ 90% of the flux concentrated in the central 2 pixels. (-99)\n- 2. Next, we remove candidates with bad radius fits, indicating a poor Gaussian fit which tends to produce solutions with unphysically large radii. These correspond to weak clusters which are not well separated from the local background, resulting in a relatively flat profile and hence a poor Gaussian fit with a large radius. We remove clusters with deconvolved diameters > 6 pc. When the images of candidates with diameters > 6 pc are examined, these objects either appear as very small, weak sources in high background regions or highly extended and not cluster-like. (-155)\n- 3. In order to present a catalog of massive star clusters, we implement a stellar mass cut of 10 4 M ⊙ assuming a ZAMS population. In addition, the distinction between star clusters and OB stars be-\n- mes ambiguous below ∼ 10 3 M ⊙ 10 . Therefore, a lower mass limit of 10 4 M ⊙ places us solidly in the star cluster regime. (-154)\n- 4. Finally, we remove one candidate which overlaps with a known background galaxy (2MASS J09555095+6940302; Skrutskie et al. 2006) and one candidate which overlaps with the galactic center and appears highly elongated rather than cluster-like. (-2)  \nAfter these cuts, the catalog (Table 1) consists of 1357 star cluster candidates. We show a mass-radius diagram, along with histograms of the masses and radii, of the final cluster catalog in Figure 3. Our catalog only captures massive (M ∗ > 10 4 M ⊙ ) clusters which are bright and somewhat extended in the F250M band. Deeply embedded, very young, or very compact clusters may be missed by our selection criteria. Future observations and analysis from this JWST program will result in a more complete census of the embedded, young clusters.  \nWe investigate whether there are trends between the measured masses and radii and the deprojected galactocentric radius ( R GC ) of each star cluster candidate. To calculate R GC , we assume that the center of M 82 is 09 h 55 m 51 . 6 s +69 · 40 ' 45 . 6 '' (J2000) (Bolatto et al. 2024) and that M82 has a major axis position angle of 67 · (Martini et al. 2018) and an inclination of 80 · (Lynds & Sandage 1963). We find that there is no corre-  \n10 Based on the slug models, which account for effects of stochastically sampling the IMF at low cluster masses, Krumholz et al. (2015) find that a 10 3 M ⊙ cluster has a bolometric luminosity of 5 ± 4 × 10 5 L ⊙ . Individual O and (massive) B stars have bolometric luminosities which overlap with this range.  \nlation between R GC and cluster radius (Spearman rank correlation coefficient, r s = -0.05). There is, however, a weak negative correlation between R GC and the stellar mass, with r s = -0.48. We note that r s = -0.41 if we assume M82 is perfectly edge-on. There is no significant difference in r s for star clusters which are detected with HST versus only detected with NIRCam.  \nForster Schreiber et al. (2003) used NIR and MIR spectroscopy at 1.5 '' (26 pc) spatial resolution to build a model of the star formation history of M 82. They find evidence for two bursts of star formation occurring 5 Myr and 10 Myr ago. From their models, the younger (older) burst produced 4 . 2 × 10 7 M ⊙ (1 . 6 × 10 8 M ⊙ ) of stars in the inner 500 pc, resulting in a total of 2 × 10 8 M ⊙ of stars from both bursts (i.e., within the last 10 Myr). The total stellar mass of our star cluster candidates (assuming an age of 8 Myr to give a lower limit; see Figure 3) is ≈ 4 × 10 7 M ⊙ . This stellar mass is a factor of 4 lower than the stellar mass expected from the older burst model from Forster Schreiber et al. (2003). On the other hand, if we assume that all of our clusters have an age of 0 Myr, we find a total stellar mass of ≈ 8 × 10 7 M ⊙ . This stellar mass is a factor of 2 larger than the stellar mass expected from the younger burst model from Forster Schreiber et al. (2003). This points to the star clusters we identify in the center of M82 having a spread in age, though we cannot yet say whether they are consistent with finite bursts or a more continuous mode of star formation.", '4.1. Comparison to other Star Cluster Catalogs': "Anumber of star cluster catalogs have been defined for M82. Here we focus on comparisons to two catalogs 11 : 1) the HST- and Keck-based catalog of McCrady et al. (2003) and 2) the HST-based catalog of Mayya et al. (2008). We refer the reader to those papers for details on observations and selection criteria. For both HST-based catalogs, we found it necessary to shift the positions to match the astrometry of the NIRCam data; these adjustments were done by-eye to match several of the brightest cluster candidates. We shifted the McCrady et al. (2003) positions by ( α, δ ) = ( -1 . 93 '' , 0 . 82 '' ) and the Mayya et al. (2008) positions by ( α, δ ) = ( -1 . 52 '' , -0 . 67 '' ). Both of these catalogs cover regions that are larger than the field of view (FoV) of our NIRCam sub640 images. Within this FoV, the McCrady et al. (2003) catalog con-  \nFigure 3. (Bottom left) The mass-radius diagram of the final cluster catalog assuming age = 0 Myr. Typical uncertainties on both quantities are ∼ 10%. Contours show the 10, 25, 50, and 75% density levels for clusters which overlap with the Mayya et al. (2008) catalog and those only detected with NIRCam (pink). The horizontal gray shaded region shows r HWHM less than the F250M HWHM PSF size. The vertical gray shaded region shows M ∗ ≤ 10 4 M ⊙ . The orange arrow shows the shift assuming age = 8 Myr. (Right) A histogram of the deconvolved cluster radii given in Table 1 (gray). The blue and pink histograms show the distribution for clusters which overlap with the Mayya et al. (2008) catalog and those which are only detected in NIRCam, respectively. The horizontal axis shows the number of clusters. The gray shaded region is the same as in the mass-radius diagram. (Top) A histogram of the stellar masses for age = 0 Myr given in Table 1 (gray). The pink and blue histograms are as in the right panel. The orange arrow shows the shift assuming age = 8 Myr. The vertical axis shows the number of clusters. The gray shaded region is the same as in the mass-radius diagram.  \n<!-- image -->  \n∗  \n/circledot  \ntains 19 star clusters and the Mayya et al. (2008) catalog contains 284 star clusters.  \nFor a NIRCam star cluster candidate to be associted with a cluster in these previous catalogs, the coordinates of the cluster should be separated by less than the sum of the cluster radii in each catalog. For the NIRCam catalog, we use the radii listed in Table 1. For the McCrady et al. (2003) catalog, we use 0.6 '' (10.5 pc), which is the quoted precision of their astrometry. For the Mayya et al. (2008) catalog, we use a radius of 0.29 '' (5 pc), as shown in Figure 1. Our cross-matching algorithm allows for multiple NIRCam sources to be matched to the same cluster in the other catalogs, and thus accounts for  \nthe source fragmentation due to our increased angular resolution.  \nOur cluster catalog recovers 100% of of the stars clusters identified by McCrady et al. (2003) and 51% of those identified by Mayya et al. (2008) (within the NIRCam sub640 FoV). Compared to the Mayya et al. (2008) HST catalog, we find 1183 previously unknown star cluster candidates. In other words, 87% of the cluster candidates we identify with these NIRCam data are new. As shown in Figure 3, the clusters that are newly detected with NIRCam tend to have smaller sizes and stellar masses than those that are also detected by Mayya et al. (2008) with HST.  \nDifferences between our NIRCam star cluster candidate catalog and that from HST (Mayya et al. 2008) are due primarily to two effects. First, that we only recover 51% of the sources from the Mayya et al. (2008) catalog is primarily due to the increased angular resolution of JWST compared to HST. As is illustrated in the lower left panels of Figure 1, there are many clusters identified with HST (light blue circles in the right panel) In this region of M 82. In the JWST image, however, many of the large, bright clusters seen with HST break apart into much smaller objects, which do not satisfy our star cluster selection criteria (Section 3.4). As a result, only 8 of the 26 clusters seen with HST in this region are classified as massive star cluster candidates here (blue circles). Second, that we detect 1183 new star cluster candidates than previous HST catalogs is primarily due to dust extinction. As is illustrated in the lower middle and right panels of Figure 1, we detect many star cluster candidates with NIRCam (pink circles) in these regions of M 82 that are not present or classified as star clusters in the HST image and catalog. We note that removing our stellar mass cut only increases the fraction of recovered HST clusters by 1%.  \nMcCrady et al. (2003) obtained Keck NIRSPEC spectroscopy on two massive star clusters in the center of M82, MGG-9 and MGG-11, allowing them to robustly measure their stellar masses and ages. Clusters MGG-9 and MGG-11 clusters correspond to our clusters 1 and 2, respectively, which are highlighted in Figure 1. Using the HST NICMOS NIC2 F160W image, which has a pixel size of 1.3 pc (0.075 '' ), McCrady et al. (2003) found half-light radii of 2.6 pc and 1.2 pc for MGG-9 and MGG-11, respectively. With our 1.8 × smaller pixel scale for NIRCam data, we measure radii of 1.8 pc and 1.6 pc, respectively. By combining velocity dispersions measured from their high spectral resolution Keck NIRSPEC data with these radii, McCrady et al. (2003) report kinematic (virial) masses of (1 . 5 ± 0 . 3) × 10 6 M ⊙ and (3 . 5 ± 0 . 7) × 10 5 M ⊙ . They also estimate the ages of these  \nFigure 4. A color-color plot of the clusters (black dots). The contours show the density of the points; contour levels show 10, 25, 50, and 75% density levels for clusters which overlap with the Mayya et al. (2008) catalog (blue) and those only detected with NIRCam (pink). The colored tracks show yggdrasil SSP models assuming no (f cov = 0; blue) and maximal (f cov = 1; red) contributions from nebular emission. The green arrow shows the slope of the reddening vector with A V = 3 -24 mag marked in increments of 3. The clusters in M82 detected with NIRCam show significant reddening with respect to the dust-free SSP tracks.  \n<!-- image -->  \n-  \nclusters to be 10 Myr and 9 Myr respectively. Assuming these ages, we find stellar masses of (0 . 9 ± 0 . 8) × 10 6 M ⊙ and (4 . 0 ± 1 . 3) × 10 5 M ⊙ for clusters 1 and 2, respectively. We conclude that the masses we derive for these clusters agree with McCrady et al. (2003) within the uncertainties.", '4.2. Color-Color Diagrams': 'To investigate the relative dust content in these candidate massive star clusters, we plot a color-color diagram using the F140M, F250M, and F360M filters. The F140M filter is most affected by dust extinction whereas we expect that dust continuum emission contributes to  \nthe F360M emission in addition to the stellar continuum 12 . Thus, these colors can give us a handle on the relative dust content of these candidate star clusters, especially when compared to dust-free SSP models.  \nWe show the AB magnitudes derived from aperture photometry, described in Appendix B, in Figure 4. We also show tracks from the yggdrasil (Zackrisson et al. 2011) models for pure stellar emission ( f cov = 0) and maximal nebular contribution ( f cov = 1). These two tracks converge for ages ≳ 7 Myr because, under the assumption of an instantaneous burst, most massive stars have died and are not replenished, so that the Lyman continuum flux is too low to sustain appreciable nebular emission.  \nWe plot the reddening vector in Figure 4 following Equation 9 of Salim & Narayanan (2020) for a range of A V . The clusters which are only detected with NIRCam show somewhat larger extinctions compared to those also detected with HST by Mayya et al. (2008), as expected (c.f. the pink and blue contours in Figure 4). In particular, the clusters which are detected with both HST and NIRCam have minimum extinctions ≈ 0 mag, whereas those that are only detected with NIRCam have minimum extinctions of ≈ 3 mag. Overall, the cluster candidates show significant reddening, with A V ∼ 3 -24 mag. This range of extinctions agrees well with the extinctions measured in the NIR by McCrady et al. (2003), who found A V ≈ 11 and 7 for MGG-9 and MGG-11, respectively (clusters 2 and 3 in our sample). For the optically-identified HST clusters, Mayya et al. (2008) report A V ≲ 6; we are able to probe more highly extincted clusters with these NIRCam data. Our range of A V is also in excellent agreement with A V = 10 ± 5 derived by Forster Schreiber et al. (2001) for a foreground screen model using multiple recombination lines in the integrated emission in the central region. In the NIRCam bands themselves, A V ∼ 3 -24 mag corresponds to A 1 . 4 µ m ∼ 0 . 7 -5 . 4 mag, A 2 . 5 µ m ∼ 0 . 3 -2 . 1 mag, and A 3 . 6 µ m ∼ 0 . 1 -1 . 2 mag.  \nThe reddening vector in Figure 4 is anchored to the f cov = 1 SSP at an age of 4 Myr. For this choice, the reddening vector is approximately aligned to the major axis of the density contours. While there is a large  \nspread in the potential cluster ages, this suggests that we are probing a relatively young, heavily reddened cluster population in the center of M 82.', '4.3. Cluster Mass Function': "The cumulative cluster mass function (CMF) measures the number of star clusters above a given stellar mass as a function of the stellar mass and may provide insights into the efficiency of star formation and links between the giant molecular cloud core mass function and the stellar initial mass function. CMFs are generally fit by a single power-law of the form dN/dm ∝ m -β , and observations generally measure an index β ≈ 2 . 0 ± 0 . 2 (e.g., Zhang & Fall 1999; Mayya et al. 2008; Krumholz et al. 2019; Emig et al. 2020; Mok et al. 2020; Whitmore et al. 2021). However, some studies suggest that there may be a cutoff of the power-law at the high-mass end (e.g., Larsen 2009; Konstantopoulos et al. 2013; Adamo et al. 2015; Messa et al. 2018), potentially indicating that there is some maximum stable cluster mass above which proto-clusters fragment, though other studies find that such cutoffs are statistical rather than physical (e.g., Cook et al. 2019; Mok et al. 2019, 2020; Whitmore et al. 2020). While most studies of star clusters in other galaxies are cataloged using optical data (e.g., broadband HST photometry; Zhang & Fall 1999; Mayya et al. 2008; Mok et al. 2020; Whitmore et al. 2020), the CMF slopes are consistent even when high-resolution submillimeter data are used, probing the most embedded, youngest star cluster populations (e.g., Emig et al. 2020). Given the resolution, sensitivity, and range of extinctions probed by these JWST NIRCam data, we are in a prime position to investigate the shape of the CMF in the extreme starburst of M 82.  \nWe construct the cumulative CMF for various assumed ages for the NIRCam-detected clusters. We follow the methodology described in Section 3.3 to calculate the stellar masses for 0, 2, 5, 8, and 50 Myr. We plot these cumulative CMFs (in terms of the fraction of clusters, rather than the absolute number) in Figure 5. We discuss the progression of the CMFs (i.e., Υ ∗ ) with age in Appendix A.  \nThe cumulative CMFs exhibit significant flattening at the low-mass end, which results from our choice to remove clusters with a ZAMS M ∗ < 10 4 M ⊙ . All of the cumulative CMFs show a steep turnover at the highmass end, which is commonly seen in other galaxies, using both optical (e.g., Mok et al. 2020) and submillimeter (e.g., Emig et al. 2020) observations proving a range of cluster ages and embeddedness.  \nWe compare our 8 Myr CMF to the CMF derived for the 'nuclear' (i.e., central, non-disk) clusters identified  \nFigure 5. The cumulative CMFs shown as a fraction of the number of clusters for different assumed ages (colored solid curves). The blue and orange shaded regions on the 0 Myr and 8 Myr curves indicates the approximate range of mass increase due to extinction which scales with M ∗ (dark) or which is constant with M ∗ (light); see Section 4.3.1 for details. The gray shaded region shows masses below 10 4 M ⊙ which are not included in the final catalog. The purple line and shaded region shows the power-law fit for the 0 Myr CMF for M ∗ > 10 4 M ⊙ . The dashed orange curve shows the cumulative CMF for the nuclear star clusters identified by Mayya et al. (2008) which have an assumed age of 8 Myr.  \n<!-- image -->  \n∗  \n/circledot  \nwith HST by Mayya et al. (2008); their cumulative CMF is shown in the orange dashed curve in Figure 5. Overall the cumulative CMFs at 8 Myr show some differences. At the low-mass end, we recover slightly more clusters compared to Mayya et al. (2008), likely due to the increased resolution and sensitivity of JWST. At the highmass end, our cumulative CMF drops off faster than Mayya et al. (2008), likely because large clusters break apart at this increased spatial resolution compared to HST.  \nWe fit our CMFs where M ∗ > 10 4 M ⊙ with a powerlaw of the form log y ∝ -α log x , where β = 1 + α . For a range of ages, we find β = 1 . 9 ± 0 . 2 (where the value is the slope of the ZAMS CMF and the uncertainty reflects the range of slopes measured for different ages), as shown in Figure 5. This is in excellent agreement with CMF slope measurements in other starburst galaxies, which tend to have β ≈ 2 (e.g., Emig et al. 2020; Mok et al. 2020; Whitmore et al. 2021), and with the HSTbased star cluster catalog of Mayya et al. (2008) who found β = 1 . 8 ± 0 . 1. However, from Figure 5 it is clear that a simple power-law is not a particularly accurate description of the CMF shape.", '4.3.1. Effects of Dust Extinction': "Because of dust, the stellar mass estimates calculated in Section 3.3 are likely underestimated since the SSP models used to derive the mass-to-light ratio do  \nnot account for the effects of dust extinction or emission. For the mean F250M AB magnitude (20 ± 1 mag), A 2 . 5 µ m = 0 . 3 -2 . 1 mag (Section 4.2) corresponds to a spectral flux density increase of < 0 . 5 mJy or a stellar mass increase of < 8 × 10 4 M ⊙ ( < 4 . 7 × 10 4 M ⊙ ) for a ZAMS (8 Myr) population. We note, however, that this change will not result in a constant shift for all masses and, therefore, could impact the measured CMF slope. The precise change in mass depends on the extinction of each individual cluster which itself likely depends on the cluster mass (e.g., de Meulenaer et al. 2013) as well as the cluster's exact age. With future spectroscopy, the age and reddening can be constrained individually, leading to much more robust mass measurements.  \nAs the most pessimistic case, we also test the change in the CMF slope assuming a constant A 2 . 5 µ m = 2 . 1 mag across all masses. We show the potential shift in the 0 Myr and 8 Myr cumulative CMFs as the blue and orange light shaded regions in Figure 5. This pessimistic case will result in a steepening of the CMF slope up to β = 3 . 2. Such a steep CMF slope is inconsistent with literature-reported values.  \nAs a more realistic way to quantify the effect of extinction on the shape of the CMF, we assume that the most massive clusters are the most extincted (e.g., de Meulenaer et al. 2013). We implement this a simple linear scaling of A 2 . 5 µ m , ranging from 0.3 mag at the lowest mass to 2.1 mag at the highest mass. We show the po-  \nential shift in the 0 Myr and 8 Myr cumulative CMFs as the blue and orange dark shaded regions in Figure 5. For such a change to the CMF, β = 1 . 8, within the quoted error.", '5. SUMMARY': "JWST has opened a new window into our understanding of star cluster formation and evolution. In this letter, we use the first NIRCam images of the prototypical starburst galaxy M 82 (Bolatto et al. 2024) to catalog and analyze the NIR-emitting star clusters in its center.  \n- 1. We identify 1357 star cluster candidates with M ∗ > 10 4 M ⊙ in the nuclear starburst of M 82 based on the NIRCam F250M images (Table 1; Figure 1). Compared to the HST star cluster catalog of Mayya et al. (2008), we find overlap with ≈ 50% of their sample and identify 1183 previously undetected star cluster candidates.\n- 2. Based on the color-color diagrams and comparing to dust-free yggdrasil SSP models, we find that the star cluster candidates we detected with NIRCam still exhibit heavy dust extinction, with A V ≈ 3 -24 mag, corresponding to A 2 . 5 µ m ∼ 0 . 3 -2 . 1 mag (Figure 4).\n- 3. We estimate stellar masses for the star cluster candidates based on the yggdrasil stellar population synthesis models and construct the CMF (Figure 5). We find a power-law CMF slope of β = 1 . 9 ± 0 . 2, in excellent agreement with studies of star clusters in other starburst galaxies (e.g., Emig et al. 2020; Mok et al. 2020).  \nSpectroscopy of these clusters in the NIR and MIR (e.g., with NIRSpec and MIRI MRS onboard JWST) will enable us to measure robust ages and extinctions for the star clusters, which is currently a key source of uncertainty in our mass measurement. In the future, MIRI imaging and spectroscopy data from this JWST GO program #1701 will enable us to find younger and more embedded star clusters than is possible with the F250M NIRCam filter alone. By cataloging the clusters identified with HST (e.g., McCrady et al. 2003; Mayya et al. 2008), NIRCam, MIRI, and the SMA (e.g., M. J. Jim'enez-Donaire et al. in prep.), we will have a complete census of the star cluster population across evolutionary stages in the prototypical starburst in M 82.  \nR.C.L. thanks Erik Zackrisson for producing the yggdrasil model magnitudes for the NIRCam filters used here. R.C.L. also thanks Sean Linden and the PHANGS-JWST/HST Stellar Populations Working Group for helpful discussions and advice. R.C.L. acknowledges support for this work provided by a National Science Foundation (NSF) Astronomy and Astrophysics Postdoctoral Fellowship under award AST-2102625. A.D.B. and S.A.C. acknowledge support from the NSF under award AST-2108140. L.L. acknowledges that a portion of their research was carried out at the Jet Propulsion Laboratory, California Institute of Technology, under a contract with the National Aeronautics and Space Administration (80NM0018D0004). R.S.K. acknowledges financial support from the European Research Council via the ERC Synergy Grant 'ECOGAL' (project ID 855130), from the German Excellence Strategy via the Heidelberg Cluster of Excellence (EXC 2181 -390900948) 'STRUCTURES', and from the German Ministry for Economic Affairs and Climate Action in project 'MAINN' (funding ID 50OO2206). V.V. acknowledges support from the ALMA-ANID Postdoctoral Fellowship under the award ASTRO21-0062. M.R. acknowledges support from project PID2020-114414GB100, financed by MCIN/AEI/10.13039/501100011033. M.J.J.D. acknowledges support from the Spanish grant PID2022-138560NB-I00, funded by MCIN/AEI/10.13039/501100011033/FEDER, EU. This work is based on observations made with the NASA/ESA/CSA JWST. The data were obtained from the Mikulski Archive for Space Telescopes at the Space Telescope Science Institute, which is operated by the Association of Universities for Research in Astronomy, Inc., under NASA contract NAS 5-03127 for JWST. These observations are associated with program JWSTGO-01701. Support for program JWST-GO-01701 is provided by NASA through a grant from the Space Telescope Science Institute, which is operated by the Association of Universities for Research in Astronomy, Inc., under NASA contract NAS 5-03127. This research made use of SourceXtractor++ 13 , an open source software package developed for the Euclid satellite project. This research has made use of NASA's Astrophysics Data System Bibliographic Services.", 'Facilities: JWST(NIRCam)': 'Software: Astropy (Astropy Collaboration et al. 2022), MatPlotLib (Caswell et al. 2020), multicolorfits (Cigan 2019), NumPy (Harris et al. 2020), pandas (Reback et al. 2020), photutils (Bradley  \net al. 2021), SciPy (Virtanen et al. 2020), seaborn (Waskom et al. 2014), SourceXtractor++ (Bertin et al.  \n2022; Kummel et al. 2022), yggdrasil (Zackrisson et al. 2011)', 'REFERENCES': 'Figure 6. The progression of Υ ∗ in the F250M filter with age from the yggdrasil models at solar metallicity and maximal nebular emission (see Section 3.3). Key phases of stellar evolution, which govern changes in Υ ∗ , are marked. See the text in Appendix A for details.  \n<!-- image -->', 'A. Υ ∗ AS A FUNCTION OF AGE': "As shown in Figure 5, the cumulative CMFs do not shift monotonically with age. This is due to how Υ ∗ varies with age in the yggdrasil models, shown in Figure 6 for the first 100 Myr, which is governed by stellar evolution. Here we present a broad-brush and highly simplified interpretation of the behavior of Υ ∗ with age in the F250M filter, where the numbers corresponds to those in Figure 6.  \n- (0) For the first ∼ 3 Myr of the cluster's evolution, Υ ∗ is roughly constant as stars age on the main-sequence (MS).\n- (1) After ∼ 3 Myr, the most-massive and most-luminous O-stars will rapidly evolve off the MS and explode as supernovae (SN). As these stars die, the luminosity of the cluster decreases resulting in an increase in Υ ∗ . Although the stellar mass of the cluster also decreases, the OB-stars dominate the cluster luminosity, so the decrease in luminosity is the dominant effect.\n- (2) From ∼ 5 -9 Myr, Υ ∗ decreases. During this time, less massive O and B stars are evolving off the MS. While the supergiant phases have roughly constant bolometric luminosity, the stars become redder in color as they cool and expand. This increases the luminosity in the F250M band. The red supergiant (RSG) phase plays a key role here, as the bolometric luminosity, and hence F250M, increases (see also section 4a of Charlot & Bruzual A. 1991). Together, these effects significantly decrease Υ ∗ .\n- (3) Between ∼ 9 -100 Myr, the luminosity in the F250M filter increases primarily due to the presence of increasing number of asymptotic giant branch (AGB) stars (see also section 4a of Charlot & Bruzual A. 1991). This increased luminosity decreases Υ ∗ .  \nIn the first 10 Myr (corresponding to the age of the older starburst in M 82; Forster Schreiber et al. 2003), the median Υ ∗ = 0 . 35 M ⊙ /L ⊙ , equivalent to Υ ∗ for a ZAMS. Therefore, the stellar masses we compute for a ZAMS are likely the best reflection of the stellar mass distribution of the star clusters in the center of M 82 when assuming a single Υ ∗ . The stellar masses derived for an 8 Myr population, on the other hand, are near the minimum possible masses, as they correspond to the minimum Υ ∗ . Therefore, the stellar masses reported in Table 1 reflect the median (0 Myr) and minimum (8 Myr) stellar masses for the star clusters in the center of M 82.", 'B. APERTURE PHOTOMETRY': "To construct the color-color diagram presented in Section 4.2 and Figure 4, we measure the cluster brightness in the F140M, F250M, and F360M images in apertures (as opposed to the flux based on the radial profile fitting described in Section 3.2 and reported in Table 1). We place circular apertures at the cluster positions in Table 1. Our aperture radii  \ncorrespond to the 70% encircled energy of the PSF ( r 70 ) given in the NIRCam APCORR Reference File 14 . For the F140M, F250M, and F360M filters r 70 corresponds to 0.065 '' , 0.092 '' , and 0.123 '' respectively. After extracting the flux in each aperture, we apply the corresponding aperture correction in the NIRCam APCORR Reference File ( apcorr70 ); for all three filters, apcorr70 ≈ 1 . 46. We note that these aperture corrections are measured for point sources. For our clusters, which are slightly extended compared to the PSF, these point-source apertures corrections are not strictly correct. However, for the purposes of this preliminary investigation, we will proceed with these corrections.  \nWe extract the background in an annulus around and with the same area as the circular aperture. We measure the median (50%), 16%, and 84% background levels in the annulus. We use the median for the background subtraction and propagate the 16% and 84% values into the uncertainties.  \nFinally, we convert the flux extracted in each aperture and in each filter to AB magnitudes 15 which we list in Table 1. The color-color diagram for the clusters following this procedure is shown in Figure 4."}