The Timescales of Star Cluster Emergence: The Case of NGC 4449

We survey the young star cluster population in the dwarf galaxy NGC 4449 with the goal of investigating how stellar feedback may depend on the clusters’ properties. Using ultraviolet (UV)–optical–near-infrared photometry obtained with the Hubble Space Telescope, we have recovered 99 compact sources exhibiting emission in the Paβ hydrogen recombination line. Our analysis reveals these sources possess masses of 102 < M ⊙ < 105, ages of 1–20 Myr, and a color excess E(B − V) in the range 0–1.4. After selecting clusters with masses above 3000 M ⊙ to mitigate stochastic sampling of the stellar initial mass function, we find that our IR-selected clusters have a median mass ∼ 7 × 103 M ⊙ and remain embedded in their surrounding gas and dust for 5–6 Myr. In contrast, line-emitting sources selected from existing UV/optical catalogs have a median mass ∼ 3.5 × 104 M ⊙ and have cleared their surroundings by 4 Myr. We further find that the environment in NGC 4449 has too low pressure to drive these differences. We interpret these findings as evidence that the clearing timescale from presupernova and supernova feedback is cluster mass dependent. Even in clusters with masses ∼ 7000 M ⊙, stochastic sampling of the upper end of the stellar initial mass function is present, randomly decreasing the number of massive stars available to inject energy and momentum into the surrounding medium. This effect may increase the clearing timescales in these clusters by decreasing the effectiveness of both presupernova and supernova feedback; neither models nor observations have so far explored such dependence explicitly. Future studies and observations with, e.g., the JWST, will fill this gap.


INTRODUCTION
Star formation occurs predominantly in clusters that emerge from the gravitational collapse of massive clouds of molecular gas and dust.When a molecular cloud cools sufficiently, its gravitational potential energy overcomes the kinetic energy of the gas, forming overdensities which then feed the formation of new stars (Krumholz et al. 2019).The formation of massive stars is however complicated by the presence of numerous counteracting mechanisms.While molecular gas, under the force of gravity, accretes onto a forming star, the star simultaneously injects energy into its surroundings (Pellegrini et al. 2011;Dale et al. 2012;Krause et al. 2013).This energy, referred to as stellar feedback, comes in several forms including stellar winds, photoionization, radiation pressure and supernovae explosions.Feedback is able to counteract the force of gravity and expel the gas and dust surrounding the region of star formation (Lucas et al. 2020;Grudić et al. 2022).When enough material has been cleared, star formation effectively ceases within the cluster.Furthermore, the energy and momentum injected into the surrounding interstellar medium drives turbulence and limits the amount of cold, dense gas that can collapse, thus contributing to the low star formation efficiency (Hennebelle & Chabrier 2011;Hopkins et al. 2012;Dobbs 2015;Goldbaum et al. 2016).The impact of feedback on star formation is so significant that the conversion efficiency of turning molecular gas into stars is on the order of a few % (Ostriker et al. 2010;Hopkins et al. 2014Hopkins et al. , 2020;;Grudić et al. 2018).Therefore, stellar feedback is not only critical to regulating star formation on a local level, but on a galactic level as well (Krumholz & McKee 2005;Ceverino & Klypin 2009;Dobbs et al. 2011;Krumholz et al. 2019).
Massive star feedback is conventionally divided into pre-supernova (pre-SN) feedback (photoionization, direct and indirect radiation pressure and stellar winds) and supernova feedback.Pre-SN acts on short timescales, typically ≲4 Myr (Pellegrini et al. 2011;Dale et al. 2012;Krause et al. 2013;Krumholz et al. 2019), and helps clear the medium before supernova explosions occur starting around 4 Myr (Leitherer et al. 2014).Some models indicate that radiative feedback may be key for regulating star formation within galaxies (Hopkins et al. 2020;Bending et al. 2022).
A number of prior observational studies have concluded that pre-SN feedback is the main mechanism for regulating star formation in clusters; these conclusions are drawn from measurements of the timescales necessary for the clusters to emerge from their natal clouds and become optically visible.Several of these studies have examined UV/optical photometry from the Hubble Space Telescope (HST) of young star clusters to derive clearing timescales <4-5 Myr and as short as 2 Myr, in nearby galaxies (Whitmore et al. 2011;Hollyhead et al. 2015;Hannon et al. 2019Hannon et al. , 2022)).However, these surveys, which only utilize optical and UV observations, are inherently limited by the effects of dust attenuation, and may miss highly dust-attenuated clusters, i.e., the clusters most closely associated with their natal clouds.Several other studies have attempted to remedy this by including CO data from ground-based observatories or 24 µm observations from the Spitzer Space Telescope (Matthews et al. 2018;Grasha et al. 2018Grasha et al. , 2019;;Kruijssen et al. 2019;Kim et al. 2021;Chevance et al. 2022).These studies have similarly derived a short clearing time, 3-5 Myr.The use of CO data or infrared images is, however, limited by low resolution.For instance, CO data typically subtends 50 pc, far greater than the 3 pc radius of young clusters (Ryon et al. 2017;Brown & Gnedin 2021).Although Corbelli et al. (2017) leverages the proximity of the Local Group galaxy M 33 (distance ∼840 kpc) to increase the spatial resolution of the Spitzer Space Telescope's 24 µm images and conclude that the embedded phase of star formation is short, ∼2 Myr, the subtended spatial scale is still 26 pc, much larger than a single star cluster.Messa et al. (2021) utilized HST near-UV to near-IR photometry to study the cluster populations of the nearby galaxy NGC 1313.They found between 40-60% of young clusters (<6 Myr) are missed when selection is only based on UV and optical data.Their analysis also extended the clearing time from 2-3 Myr to 3-4 Myr.Messa et al. (2021) is however restricted by the requirement that sources still need to be detectable at optical wavelengths in order to determine their ages and masses.Calzetti et al. (2023), using HST optical and IR images of NGC 4449, attempted to obtain a sample of highly extincted clusters by selecting sources that emit in Paβ but are undetected at optical and shorter wavelengths.The specific selection requirements used in that study resulted in a limited number of sources (34) with a narrow range of properties, which leaves open the question of whether those authors' approach includes potential biases.
The studies above, with their variety of timescales, hint at the possibility that the clearing timescale may be influenced by the galactic environment and/or by the cluster's properties.The dependence of pre-SN feedback on galactic environment has been investigated by McLeod et al. (2021), Della Bruna et al. (2022) and Chevance et al. (2022) for a few local spiral galaxies.These authors find that the galactic environment does have an influence on the effectiveness of pre-SN feedback.While Chevance et al. (2022) conclude that presence of morphological features introduce variations in the pre-SN clearing timescales, McLeod et al. (2021) find a dependency on galactocentric distance.McLeod et al. (2021) conclude that supernova feedback is enhanced by the clearing done by pre-SN feedback in the outer, less dense regions of galaxies.Similarly, Della Bruna et al. (2022) find that the HII regions in the center of the starburst galaxy M83 are confined by the ambient pressure, while they are over-pressurized and expanding in the disk.
In this work we extend the analysis of Messa et al. (2021) to the nearby Magellanic Irregular dwarf galaxy NGC 4449.NGC 4449 is located at a distance of 4.2 Mpc (Tully et al. 2013), implying that our HST nearIR imaging data resolves regions as small as ≲2-3 pc, i.e., comparable to a cluster's size (Ryon et al. 2017;Brown & Gnedin 2021).With a star formation rate of approximately 0.5 M ⊙ yr −1 and a stellar mass of ∼10 9 M ⊙ , NGC 4449 is located about 3σ above the Star Formation Main Sequence for local galaxies (Cook et al. 2014) which qualifies it as a starburst.The metallicity of the galaxy is about 40% solar1 and has a modest galactocentric gradient (Berg et al. 2012;Pilyugin et al. 2015).There are several similarities between NGC 4449 and NGC 1313 (the latter studied by Messa et al. 2021).The two galaxies have comparable distances, oxygen abundances and specific SFRs (SFR/M * ) (see, Messa et al. 2021).Within a factor ∼3, their masses are comparable to that of the Large Magellanic Cloud.However, they differ in their SFR surface density, with NGC 4449 having Σ SF R = 0.03 M ⊙ yr −1 kpc −2 , or three times the SFR surface density of NGC 1313 in the regions targeted by HST (Messa et al. 2021;Calzetti et al. 2023).
We leverage the same HST UV-to-nearIR imaging data analyzed in Calzetti et al. (2023) (Figure 1), but adopt a more general selection process than those authors to recover a larger sample of young star clusters, with a greater diversity of ages, masses and extinctions.We select our sources based on the presence of ionized gas emission identified from the infrared hydrogen recombination line Paβ (λ1.282 µm); this ensures that our sources are typically young, ≲6-7 Myr, since beyond 7 Myr the ionized gas emission is too faint to be detected (Leitherer et al. 1999).However, we do not impose explicit constraints on mass or extinction, thus mitigating selection biases on these parameters.
The paper is structured as follows: section 2 contains a description of the imaging data, section 3 describes the construction of the source catalog, and section 4 details the cross-validation of the sources using a JWST image of NGC 4449.Section 5 describes the multi-wavelength flux measurements using aperture photometry for the sources.The stellar population and dust attenuation models are described in section 6. Results and discussion are in section 7, while the summary and conclusions are in section 8.

IMAGING DATA
The data used in this survey is comprised of HST imaging covering from nearUV (NUV) to nearIR (NIR) in 10 bands.With the WFC3/UVIS and ACS instruments, we sample the NUV to optical range with broad-band filters (F275W, F336W, F435W, F555W and F814W), a Hα+[NII] narrow-band filter (F657N), and a medium-band filter sampling the line-free continuum at V (F550M) (Table 1).Further, with WFC3/IR we sample the NIR with two broad-band (F110W and F160W) and one narrow-band (F128N) filter to cover the Paβ hydrogen recombination line emission.In summary, the 10 filters spanning the NUV-NIR include two hydrogen recombination lines (Table 1).For each IR filter, the standard calibration pipeline CALWFC3 v. 3.5.2was used to process individual frames into the final images, after correction for bias, dark, and flat fields (Messa et al. 2021).The NUV and optical images were retrieved from the MAST 2 archive already processed through their respective instrumental pipelines.All images were aligned and mosaiced to the Gaia DR2 reference frame and then resampled to 0.04"/pix (NUV/optical) or 0.08"/pix (NIR), using DrizzlePac.The NIR images are resampled to twice the pixel size of the NUV/optical images because the Point Spread Function (PSF) of the WFC3/IR images is twice that of the optical and UV images.At the distance of NGC 4449, the pixel scale of 0.04" corresponds to 0.81 pc.Using the header keyword PHOTFLAM, images were then converted from instrumental units (e − /s) to physical flux units (erg/s/cm 2 / Å).Table 2 contains each filter's respective pivot wavelength and FWHM.Since the WFC3/IR images have the narrowest field of view of the sample, the source search area of NGC 4449 is limited by their coverage (Figure 1).This area covers the central 2.8×2.4 kpc 2 of NGC 4449 and captures 67% of the total SFR of the galaxy.
We produce emission line images by subtracting stellar continuum from the narrow-band images, F658N and F128N.For the F658N image (Hα+[NII]), we use the F550M and F814W images to perform an interpolation of the continuum flux.For the F128N (Paβ) image, we use the F110W and F160W images for the interpolation and the resulting continuum maps are then scaled using isolated foreground stars within the images which do not have a contribution from any nebular emission.The presence of Paβ in the F110W filter requires the subtraction to be performed iteratively to adequately remove the continuum from the final line image.After producing the first F128N line emission map, we subtract this map from the original F110W image to create a line-free 1.1um continuum image.We then repeat the interpolation, scaling, and subtraction process until the final image differs from the previous version by less than 1%; this is achieved with two iterations.The continuum subtracted images are then multiplied by their filter widths, given by their full width at half maximums, and corrected for the filter's transmission at the galaxy's redshift to obtain line fluxes.The line emission in the F658N image is corrected for the [NII] contribution using a [NII]/Hα of 0.11 for NGC 4449 (Berg et al. 2012).

SOURCE CATALOG
Within the area of NGC 4449 determined by the field of view of the NIR pointing (Figure 1), we identify sources that have compact Paβ recombination line emission, therefore, are potentially young star clusters.This approach is similar to the one employed by Messa et al. (2021).Those authors use a second method for specifically isolating dusty sources by searching for high-extinction compact sources in extinction maps.Here we only concentrate on the compact Paβ method, because we are interested in isolating dusty sources that are also young.There is considerable overlap between the catalogs selected with the two methods: 65% of the sources from the extinction map method are in common with the sources from the compact Paβ method; the latter method is also very efficient, since it retrieves 84% of the sources in common between the two catalogs (Messa et al. 2021).
Star cluster candidates are extracted from both the ACS F814W image and the continuum-subtracted Paβ line emission image at their respective pixel scale (0.04"/pix and 0.08"/pix, respectively).We use SExtractor (Bertin & Arnouts 1996) to perform the source extraction.Input parameters are optimized to identify sources with a ≥5σ detection in 5 and 3 contiguous pixels in the F814W and F128N, respectively.This identifies 224,342 and 1240 sources in the F814W and Paβ images, respectively.We then match the two catalogs within a 2 pixel (0.08" or 1 pixels in the IR channels) tolerance, retaining 1119 common sources.Visual checks were performed to ensure that the low tolerance was not discarding potential matches.This initial catalog is then inspected by eye, using the full suite of images.Apart from image artifacts and edge detections, the main source of contamination in this catalog are poorlysubtracted sources which appear as detections in our Paβ line map with Signal-to-Noise ratio S/N≥ 5 relative to the locally-measured background.In regions where the stellar continuum is over(under)-subtracted, the sources that are not perfectly subtracted will be preferentially removed (retained) in our catalog.These sources are straightforward to visually identify and remove from the catalog as they include negative pixel values and live in regions where the overall background is over-subtracted.Since these false detections arise due to poor continuum subtraction rather than genuine Paβ emission, omitting them does not bias our catalog against additional young clusters, nor does the issue of over-subtraction mask any real sources.We also remove sources whose peak Paβ emission does not appear coincident with the peak stellar continuum emission, since these may be the result of several sources in superposition in a compact association or may be slightly evolved sources that have ejected their surrounding gas.After this visual inspection we retain 464 sources.

CONFIRMATION OF SOURCES WITH JWST DATA
We cross reference the locations of the remaining 464 sources to a stellar-continuum subtracted Paα (λ 1.8756 µm) line emission map acquired with JWST.The JWST Paα image has been obtained using the NIRCam F187N filter, continuum-subtracted using adjacent broad-band filters (private communication, A. Adamo 2023).This image has double the angular resolution and is deeper than our HST/WFC3/IR images, thus providing a sanity check for our nebular emission detections.The Paα hydrogen recombination line serves as an additional marker for young star clusters while being less susceptible, by virtue of being at longer wavelength, to the effects of dust attenuation than Paβ.We perform a visual inspection of the Paα image at the location of our sources, discarding all sources that may be poorly subtracted stars.In addition, we impose a detection limit cut-off and only deem acceptable those sources that are detected in Paα with S/N≥3.This last criterion removes 167 faint emission line sources, which may be either old clusters or young low-mass ones.We do not expect that discarding low-mass young clusters will affect our conclusions, which are based on sources with mass≥3,000 M ⊙ (see section 7).We note that we do not recover any of the sources identified by Calzetti et al. (2023); none of the 34 sources selected by those authors show evidence for emission in Paα, although they are sufficiently compact to be indistinguishable from stellar sources and therefore part of the stellar population in this galaxy.After applying both of our selection criteria, we are left with 99 sources (21% of the original sample) that display Paα emission, with >90% of them detected with S/N>7 in the line.This is the sample our analysis centers on.

PHOTOMETRY
Aperture photometry is performed on each of the 99 cluster candidates.We adopt a fixed 3/1.5 (NUV+optical/NIR) pixel radius (0 ′′ .12=2.44 pc) to ensure that photometry is not affected by nearby contaminants.A local sky annulus in the range 5-7/2.5-3.5 (NUV+optical/NIR) pixels is used for the background subtraction of each source.Prior to computing the photometry, each source is centered on the F814W using a centroiding algorithm and fixed in all filters to this point.The F814W is used as reference since it is our reddest optical band, with a point spread function (PSF) roughly half the size of our NIR bands.This makes the F814W ideal as a reference since it can capture the flux of highly extincted sources compared to shorter bands while still having sufficient angular resolution to differentiate between the source and nearby contaminants.
We account for flux outside of the fixed aperture by applying aperture corrections which were determined using the methods described in Messa et al. (2021).Briefly, to estimate the corrections, we convolve isolated stellar sources in each band with a Moffat profile, considered to be the most accurate light profile function for describing clusters (Bastian et al. 2013).We use a R eff value of 2.5 pc, which is the peak of the cluster radius distribution (Ryon et al. 2015(Ryon et al. , 2017)), to obtain an average correction in each band for all clusters.Although the choice of R eff affects the normalization of the SED, it does not affect its overall shape.Thus, when fitting the SEDs to models, the only physical property that is affected is the resulting mass of the cluster.In addition, Messa et al. (2021) found that the resulting change in mass distribution between R eff values of 1.3 pc (the average radius of clusters in their sample) and 2.5 pc is minimal.We apply the following corrections (all expressed in AB mag): -1.54±0.09for F160W, -1.52±0.09for F128N, -1.39±0.09for F110W, -1.19±0.08 for F814W, -1.16±0.08 for F658N, -1.13±0.08 for F555W, -1.08±0.08 for F550M, -1.13±0.08 for F435W, -1.03±0.08 for F336W, -1.04±0.08 for F275W.Finally, we correct the photometric measurements for foreground galactic extinction with the MW foreground extinction, E(B-V)=0.017.Total uncertainty is derived from a combination of photon noise, uncertainty in aperture corrections, and the standard deviation of the background.We check for relative (band-to-band) photometric accuracy between broad and narrow-band filters by measuring several stellar sources in the field of NGC 4449.The stellar sources do not show, or in a few cases only have non-significant (<1 σ), emission lines, either in Hα or Paβ, indicating that the narrow-band filters are accurately calibrated within a few percent, and we should not expect to detect 'spurious' emission lines due to mis-calibrations.Photometry of the 99 sources is listed in Table 3 and their location is shown in Figures 1 and 4.

MODELS AND FITTING
We fit our photometry to synthetic models to derive physical properties of the cluster candidates including age, mass, and extinction.We utilize Yggdrasil SED models (Zackrisson et al. 2011) in combination with dust attenuation/extinction models to produce synthetic photometry similar to the procedure followed in Calzetti et al. (2015a), Adamo et al. (2017), Messa et al. (2021), andCalzetti et al. (2023).
Yggdrasil generates single stellar population (SSP) models by utilizing the Starburst99 (Leitherer et al. 1999) spectral synthesis models and then combining these models with Cloudy (Ferland et al. 2013), to generate the nebular emission lines (Zackrisson et al. 2011).The Starburst99 models already include both stellar and nebular continuum.
The SSP models we use are generated assuming instantaneous star formation, a Kroupa Initial Mass Function (IMF) (Kroupa 2001) in the range 0.1-120 M ⊙ , and metallicity of Z=0.008, consistent with the metallicity of NGC 4449, which $JH0\U is about 40% of the solar value (Berg et al. 2012).Ages are sampled in 1 Myr increments between 1 Myr and 10 Myr, with increasingly sparser age sampling out to 10 Gyr.We discuss later the impact of stochastic sampling of the stellar IMF on our source parameters determinations (Cerviño et al. 2002).
The SSP SEDs are attenuated with both a starburst attenuation curve (Calzetti et al. 2000) and an LMC extinction curve (Fitzpatrick 1999).We use a foreground dust geometry (Calzetti 2001) of the form: where E(B-V) is the color excess and κ(λ) is the attenuation/extinction curve.The models range in E(B-V) from 0-5 mag with increments of 0.02 mag.The dust-attenuated SED models are then convolved with the HST filters' transmission curves to produce synthetic photometry.
We then fit the SEDs derived from aperture photometry to these synthetic photometry models using χ 2 -minimization, taking into account the uncertainties in each band as well as upper limits.A source must possess at least 3 bands with σ<0.3 for the source to be fit.To account for highly extincted sources with weak, highly uncertain detections at shorter wavelengths, we utilize the uncertainty in the F555W as a reference point.If σ F 555W exceeds 0.3 mag, we treat the F555W and all shorter wavelength bands as upper limits during the fitting process.This approach resolves the issue of poor fits caused by a lack of detection at shorter bands.With this criteria we were able to successfully fit all of the 99 sources.î î î :DYHOHQJWKÅ $%0DJ ) : ) : ) : ) 0 ) : ) 1 ) : ) : ) 1 ) : Fitting returns a distribution of reduced-χ 2 (χ 2 red ) values for each model and its corresponding physical properties.We take the model with the minimum χ 2 red to be the best fit.With this, we find that the models with the LMC extinction curve produce better fits than the models with the starburst attenuation curve, in agreement with the findings of Calzetti et al. (2023).We thus adopt the models with the LMC extinction as our reference models.We determine an uncertainty interval for the best-fit parameters of each source to be all fits within 1σ of the minimum χ 2 red value, a corresponding χ 2 red interval of 3.53 based on our 3 degrees of freedom.The list of cluster properties, ages, masses and extinctions, with 1 σ uncertainties and the value of the minimum χ 2 red from the fits is given in Table 4 for all 99 sources.Figure 2 illustrates the resulting distribution of best-fit ages and minimum χ 2 red of our sample.We define a minimum χ 2 red cutoff of 10 to ensure a well-fit sample, following the approach of Adamo et al. (2017).Of the 99 sources with sufficiently accurate photometry to be fit with our method (see above), 80 sources possess χ 2 red < 10. Figure 3 illustrates four SEDs and their best-fits with χ 2 red < 10.These four sources span the whole parameter space of age and color excess in our sample.The location of the 80 sources with χ 2 red < 10 is shown in Figure 4 (upper-right panel).

RESULTS AND DISCUSSION
We analyze the correlation between physical properties obtained from the SED fitting, presenting samples selected with different cuts, and discuss their implications.
Figure 5 illustrates the age and color-excess relation of the sources in our sample, using both the entire sample of 99 sources (left) and the sub-sample of 80 sources with χ 2 red < 10 (right).We observe a plume of sources ranging in E(B-V) from 0 mag to ∼1.4 mag centered at 5-7 Myr.We also observe a limited number (16 total, 14 with χ 2 red <10) of sources with ages ≤3 Myr, none of which possesses E(B-V) exceeding 0.56 mag.The overall shape of both the full sample and the one with χ 2 red < 10 are qualitatively similar.The age-mass plots (Figure 6) show that the 5-7 Myr plume has best-fit masses in the range ∼10 2 -10 4.2 M ⊙ with higher masses, up to 10 5.5 M ⊙ , for the few >20 Myr sources.The low mass end is consistent with the mass limit we can reach given the depth of our images.We derive the mass limit as a function of age from the I-band detection limit, and its track is shown in Figure 6.
While young (≲6 Myr), moderately dusty star clusters tend to have strong Paβ line emission, which is well modeled by our fitting routine (top panels of Figure 3), the fitting routine has difficulty discriminating between the SED of a dusty, ∼6-8 Myr old cluster (Figure 7) and the SED of a much older, but low dust, cluster, leading to some degeneracy in the fits; this is clearly illustrated in the two bottom panels of Figure 3, where a dusty 6 Myr SED (left) and an almost-dust-free 9000 Myr SED (right) are shown.Four sources in our sample (three with χ 2 red < 10) are older than 20 Myr and as old as ∼10 Gyr, according to their best fits.The requirement that sources be detected in Paβ line emission with S/N≥5 in SExtractor is intended to break the degeneracy between "young and dusty" and "old and low-dust" sources.Further, a total of 10 sources with χ 2 red <10 have large age uncertainties (>1 Gyr); these sources are all 6 Myr and older in age, including the three sources with best fit ages >20 Myr, again demonstrating the limitations of the fitting routine when emission lines are weak.An additional consideration is that clusters around the 6-10 Myr age range can contain red supergiants whose brightness can over-shine their host cluster, further diluting evidence for a young stellar population (see discussion in Whitmore et al. 2020;Calzetti et al. 2023).Thus, the four sources with best fit age >20 Myr may be younger than the best fits indicate.For sake of simplicity, in what follows we neglect further consideration of the four sources with best fit ages older than 20 Myr, reducing our sample to 95 sources.
We test for biases in the SED fit results due to potential chance alignment of, e.g., a red supergiant star in front of an otherwise low-dust star cluster (e.g., Whitmore et al. 2020).As red supergiants appear around ∼6 Myr of age, this configuration would produce an unusually red SED which may be mistaken for a dusty one.We test this possibility by performing spot-checks on our sample: we measure the extinction values from the line emission ratio Hα/Paα for three clusters with E(B−V)>1 mag (# 13, 45, and 82), finding that the nebular emission yields E(B−V) values comparable to those obtained from the SED fits.We conclude that the results from the SED fits are robust in terms of the physical parameters they produce.
To mitigate the effect of stochastic sampling of the IMF on the sources' derived physical parameters, we also define a 3000 M ⊙ mass cutoff which retains 22 sources with any χ 2 red and 15 sources with χ 2 red < 10 (Figure 4).The age and color-excess of the resulting sample are presented in Figure 8.The 5-7 Myr plume is still visible, but the sample only contains three (one) sources with age of 4 Myr or less.

Where are the clusters younger than 5 Myr?
Our source extraction parameters (detection in I-band and Paβ line emission) were chosen to ensure an unbiased sample of young clusters.Despite this, we recovered relatively few very young (≤4 Myr) sources and just one with mass exceeding 3000 M ⊙ for the χ 2 red <10 sample.In order to understand whether our selection method is biased against these very young and very dusty (E(B-V)> 1 mag) clusters, we perform a by-eye search of line-emitting sources that are not already in our catalog.The requirement to have line emission ensures that the sources are young.In order to maximize the likelihood that the sources are also dusty, we perform the search around detected radio sources that are located within infrared-emitting regions.The radio sources are from the collection of Reines et al. (2008), who have identified compact radio-emitting regions in NGC 4449 with an accuracy of ∼1 ′′ .Radio emission is less subject to dust extinction than optical and near-infrared wavelengths, so it can ensure an unbiased sample.The only shortcoming is that radio observations are limited to bright regions, i.e., massive star clusters.However, we search in the surroundings of these bright clusters for fainter (lower mass and/or more extincted) clusters.The infrared images we use to locate dusty regions are from the IRAC instrument on the Spitzer Space Telescope, centered at 8 µm, a wavelength that targets dust emission; we use here the image presented in Calzetti et al. (2018).Our by-eye search yields 32 sources, only 18 of which are new, i.e., not already in our sample.Of these, 13 were successfully fit, and only a single source has a best-fit χ 2 red < 10 and mass >3,000 M ⊙ , but it is relatively unextincted with an E(B-V)=0.28mag.The location of this source on the E(B-V)-age plot is shown in Figure 9.Of the remaining 12 fitted radio sources, 10 are low mass (<3,000 M ⊙ ) and 2 have poor χ 2 red >15 and low E(B-V).These findings illustrate that our ability to recover both young and dusty sources in NGC 4449 is not limited by our HST WFC3 UVIS/IR observations, whether through visual inspection or automated searches.
We investigate whether our choice of models may drive some of these results.Whitmore et al. (2020) compares the Yggdrasil models (Zackrisson et al. 2011) with the Bruzual & Charlot models (Bruzual & Charlot 2003), finding that the two sets of models generally yield comparable ages and masses.When discrepancies occur, they are in the direction of Yggdrasil assigning older ages (>40 Myr) to clusters that the Bruzual & Charlot models yield young ages (<10 Myr) for.We believe this should not affect our sample, which is selected to have nebular line emission and, thus, young ages by default.To verify this, we fit a segment of our sources using the fitting code CIGALE (Boquien et al. 2019), which includes the Bruzual & Charlot models and nebular emission, both lines and continuum.We focus on the 11 sources with χ 2 red < 10, masses >3000 M ⊙ , and ages between 5 and 20 Myr.The best fits we obtain from CIGALE indicate shifts towards marginally younger ages, down to 5 Myr, increased E(B-V) values up to 3 mag, and up to a ∼2.5 factor increase in mass.This test demonstrates that the fundamental picture remains consistent, marked by an absence of very young (<5 Myr), extincted sources.
In addition, similar tests performed in Calzetti et al. (2023), using the stellar population inference code Prospector (Johnson et al. 2021), yielded no major changes or systematics in the ages/masses/extinctions of the sources relative to the results from Yggdrasil.In light of this discussion, we conclude that the choice of models/codes does not have a major impact on our ages and extinction values.
The considerations above leave us with two possibilities: that compact and dusty star clusters with age ≲4 Myr are either absent in NGC 4449 or have been somehow missed by our searches.To test either of these possibilities we investigate existing results on both the star formation history (SFH) of NGC 4449 and the optically-selected young star clusters from the LEGUS project (Calzetti et al. 2015b).
The SFH of NGC 4449 has been derived by McQuinn et al. (2010), McQuinn et al. (2012) and Sacchi et al. (2018), using resolved stellar populations detected in the ACS/WFC optical images from the Hubble Space Telescope, and by Cignoni et al. (2018), using both UV (WFC3) and optical (ACS/WFC) images.All authors recover similar shapes for the SFH, although there are differences in the details, especially in the youngest age bins.We concentrate on the results of Sacchi et al. (2018), as these authors include some of the most recent stellar evolution models in their analysis (similarly to Cignoni et al. 2018) and divide the galaxy in sub-regions, helping a better match with our WFC3/IR footprint.Our footprint coincides with the 'Center' region and part of the 'Clump1' and 'Clump2' regions identified by Sacchi et al. (2018).Sacchi et al. (2018) finds that those central regions experienced a mild burst of star formation, that begun about 20 Myr ago, ended about 5-6 Myr ago and peaked around 14 Myr ago.Since the end of this burst, about 5-6 Myr ago, star formation in the central region of NGC 4449 has been negligible, as inferred from resolved star populations.This result is consistent with the one obtained by Cignoni et al. (2018) who recovered the SFH of the whole galaxy over the past ∼200 Myr.However, as the earliest stages of star formation produce clustered systems (unresolved star clusters and associations, see Lada & Lada 2003;Elmegreen 2010;Elmegreen & Hunter 2010) analyses that concentrate on resolved stars may miss such sources, and underestimate the SFH in the youngest age bins and specifically over the most recent few Myr.
To further test this explanation, we investigate the young star clusters identified in NGC 4449 by the LEGUS project (Calzetti et al. 2015b) and discussed in Calzetti et al. (2023).The LEGUS star cluster candidates were selected from UV,U,B,V,I images obtained with the WFC3/UVIS and ACS/WFC instruments on the HST.Whitmore et al. (2020) added candidates to the NGC 4449 cluster catalog based on inspection of the HST Hα images.The requirement that a cluster needs to be detected in at least four bands in order to be attributed an age, mass and extinction from the SED fitting of its photometry constrains the selection to clusters that are detected in at least the UV (WFC3/F275W) or U (WFC3/F336W) (Adamo et al. 2017).This sets also an upper limit to the maximum extinction of the clusters, determined by Whitmore et al. (2020) to be E(B-V)≲1 mag.The approach used by Adamo et al. (2017) to the SED fitting is similar to the one adopted in this work, enabling a direct comparison.Calzetti et al. (2023) isolated 53 star cluster from the LEGUS catalog, with the following characteristics: mass≥3,000 M ⊙ , age≤10 Myr, and χ 2 red ≤10 from the SED fits.
We inspected the location of these 53 star clusters on the narrow-band images Hα and Paβ to identify all LEGUS young clusters that are coincident with hydrogen recombination line emission, either compact or filamentary.We found only 21 such star clusters.The physical parameters (age, mass, extinction) of these 21 LEGUS cluster are plotted in Figures 6 and 9, where we include the small offset in mass due to the different choices in the distance of NGC 4449 (3.9 Mpc for LEGUS versus 4.2 Mpc in this work).An interesting property of the LEGUS line-emitting clusters is that they are all ≤4 Myr old and they extend to high masses, with median mass ∼3.5×104 M ⊙ .For comparison, the high-fidelity (χ 2 red ≤10) Paβ-emitting clusters we identify in our sample with mass ≥3,000 M ⊙ and age ≤10 Myr (8 clusters) have ages ≳5 Myr (with one at 4 Myr age) and median mass ∼7×103 M ⊙ , albeit with large uncertainties.Thus, our Paβ-emitting compact clusters are both systematically older and are a factor ∼5 less massive than the line-emitting clusters in the LEGUS sample.
Figure 10 shows that the ≤10 Myr old, Paβ-emitting sources are also on average dustier than the line-emitting LEGUS sources, as expected from the selection criteria of each sample.While both sets of sources can have low extinction values, the LEGUS clusters have a maximum E(B−V)=0.85mag, with a median value Ē(B − V )=0.4 mag, while the Paβ-emitting sources in our sample have maximum E(B−V)=1.44mag with median Ē(B − V )=1.04 mag or a little over twice that of the LEGUS clusters.The age-dust degeneracy may slightly affect the mass assigned to a source if the age is incorrect.Much of this uncertainty is captured in the 1σ uncertainties reported for our sources.For reference, in the age range 1-6.5 Myr, the mass inferred from the I-band varies by 0.33 dex and from the H-band by 0.19 dex, too small of a variation to affect the results presented in Figure 10.
In addition to differences in median age, mass, and extinction, the two samples of star clusters also differ in numbers.We find 21 line-emitting LEGUS clusters and 8 Paβ-emitting clusters that are ≤10 Myr old, ≥3,000 M ⊙ and with SED fits yielding χ 2 red < 10.However, when limiting the comparisons to the same mass range (3,000-15,000 M ⊙ ) for both samples, the numbers become comparable, since there are only 10 LEGUS clusters in this mass interval.Thus, there are no major discrepancies in the numbers recovered in the two samples, when the mass range is controlled for.One implication is that about half of the LEGUS clusters that are still retaining their gas are more massive than 15,000 M ⊙ , but they require larger dust column densities ( Ē(B − V )=0.54 mag) than their lower mass counterparts ( Ē(B − V )=0.28 mag, Figure10) to do so.
Why have we missed the LEGUS young clusters in our catalog?The line-emitting LEGUS star clusters are typically more extended and/or too faint in Paβ than allowed by our selection criteria, leading to these clusters being missed by our detection algorithm.In addition, many clusters are coincident with filamentary and/or diffuse line emission and are excluded by our matching criteria.
We now consider the sample of young, UV/optically-selected sources derived by Whitmore et al. (2020), who use an approach to the SED fits that differs from the default LEGUS one (Adamo et al. 2017).In Whitmore et al. (2020), the five-band LEGUS photometry is fit with Bruzual & Charlot models (Bruzual & Charlot 2003), but without inclusion of nebular continuum or line emission; the authors add a prediction for the Hα emission from the ionizing photons flux in the models, that they use to fit the photometry in the F658N filter.With these assumptions, Whitmore et al. (2020) derive a slightly different set of young sources for NGC 4449 from the one discussed above.We find 18 sources from the catalog of Whitmore et al. (2020) that are good quality (flag Qual=1), younger than 10 Myr, more massive than 3,000 M ⊙ , 3 and are within the WFC3/IR footprint.All 18 sources are already present in the LEGUS catalog, but with several differences in the derived physical properties.Seven of the >3,000 M ⊙ clusters from the catalog of Whitmore et al. (2020) have masses <3,000 M ⊙ in the LEGUS catalog.This is a direct consequence of the approach implemented by Whitmore et al. (2020): neglecting the nebular continuum at young ages yields overestimated stellar masses (Reines et al. 2010).Of the remaining 11 clusters, seven have derived ages ∼5 Myr in Whitmore et al. (2020) and 10≤age≤30 Myr in LEGUS.Visual inspection of the Hα image reveals that only one of the seven is truly associated with Hα emission; the remaining six clusters are contaminated by neighboring strongly emitting line sources and do not appear to be associated with in-situ Hα emission.The bona-fide Hα-emitting source is ∼5.5 Myr old, with mass ∼2.5×10 5 M ⊙ and E(B-V)∼0.44mag.The remaining 4 sources (of 18) have derived ages ∼4-5 Myr in the catalog of Whitmore et al. (2020) and ∼2-3 Myr in LEGUS, suggesting, again, that differences in the handling of the nebular continuum and lines may be at the root of the observed differences.The overall picture remains, however, unchanged, as already discussed: the young, UV/optically-selected sources have ages≲5 Myr and low extinctions when associated with line emission.We retain below the direct comparison with the LEGUS results, as our approach to SED fitting is similar to the one adopted in that project.
The picture that emerges from this analysis is that, in the central region of NGC 4449, star formation has continued for the past few Myr, producing mainly star clusters.Massive clusters retain their gas (and dust) for the first ∼4 Myr and then disperse it.Lower mass clusters can retain their gas and dust for a longer period of time, and appear compact (in line emission) and dusty up to ∼6 Myr and possibly longer.The age limit for our analysis, ∼7 Myr, is discussed earlier and is due to two concomitant effects: (1) beyond 7 Myr line emission becomes too faint to be recognized as such, and ( 2) there are age degeneracies in the SED fits for ages >6-7 Myr.
We observe that the massive clusters in the LEGUS sample retain their gas only for the first ∼4 Myr, implying that timescales may be too short to observe their dustier progenitors.Conversely, we speculate that the progenitors of the lower-mass, Paβ emitting, compact clusters we identify are buried in dust and are undetectable in our near-IR bands, similar to what occurs in the Milky Way (Gutermuth et al. 2009;Megeath et al. 2022).This may account for why we do not find young and dusty star clusters in our sample.Alternatively, stochastic sampling of the IMF may affect our age determinations (see next sections), although we attempt to mitigate its impact by selecting clusters with mass ≥3,000 M ⊙ .

Comparison with NGC 1313
The study of Messa et al. (2021) of the late-type spiral NGC 1313 is the closest one in methodology to the present study.These authors also find that the Paβ-selected young star clusters in NGC 1313 are low-mass systems and show a plume of high extinctions concentrated around 4-10 Myr, and that the LEGUS star clusters for this galaxy identify typically lower-extinction sources.However, Messa et al. (2021) do not isolate the LEGUS clusters that are coincident with ionized gas emission and do not compare the LEGUS stellar masses with the stellar masses of their Paβ selected sources.Despite these differences, those authors conclude that when including infrared-selected star clusters the clearing timescales extend from ∼2 Myr to ∼4 Myr.Our timescales appear slightly longer, ≈6 Myr, for the Paβdetected clusters.This difference, if confirmed by further studies, may be due to the higher pressure environment of NGC 4449 relative to that of NGC 1313.Based on the SFR surface density presented in the Introduction (section 1), the environmental pressure in NGC 4449 is about 2-3 times higher than in NGC 1313.

Consequences for the Timescale of Cluster Emergence
Our relatively unbiased sample of Paβ-emitting compact sources in the galaxy NGC 4449 shows a dearth of sources younger than 4 Myr, while, at the same time, we find that sources in the age range 5-6 Myr can display high values of dust extinction, up to E(B-V)∼ 1.4 mag.As discussed in section 4, none of the sources identified in Calzetti et al. (2023) with E(B−V)≳2 mag are confirmed by our more detailed inspection using the JWST Paα image, although we still find, like those authors, that we do not recover very young (≤4 Myr) sources.Conversely, the UV(U)-bright clusters identified in the LEGUS sample are associated with line emission only when ≲4 Myr old.
As presented in Section 3, the original sample from SExtractor was visually inspected to remove stars, asterisms, image artifacts, and also sources where the Paβ line emission is not coincident with the stellar continuum emission.This last criterion de facto removes from the sample sources where the ionized gas has been displaced from the location of the ionizing star cluster.This can happen, for instance, if supernova or pre-supernova feedback has occurred in the cluster, indicating that our Paβ-emitting sample will be missing sources where feedback has been active.
Pre-SN feedback acts on short timescales (<4 Myr) before the occurrence of supernovae around ∼4 Myr (Pellegrini et al. 2011;Dale et al. 2012;Krause et al. 2013;Krumholz et al. 2019;Leitherer et al. 2014).In star clusters where pre-SN feedback is efficient, the surrounding ionized gas is ejected, creating a 'shell-like' morphology for the gas emission (Whitmore et al. 2011;Hollyhead et al. 2015;Hannon et al. 2022).These clusters are likely to be young (≲4 Myr), but also have low extinction, which explains why they have been prominently identified in UV-optical samples (Messa et al. 2021;Hannon et al. 2022).
When selecting for compact-emission sources in the near infrared, like our sample, we succeed in locating star clusters with a large range of dust extinction properties, including very dusty ones with E(B-V)>1 mag.However, this sample of compact, line-emitting sources lacks very young (≲4 Myr) clusters, raising the question of how effective pre-SN and supernova feedback is in the ∼3,000-15,000 M ⊙ mass regime of our clusters.All of our clusters contain massive stars (presence of ionized gas emission), yet they show compact morphology and include many dusty sources, suggesting that neither pre-SN feedback nor supernova feedback are effective clearing mechanisms in these sources.
Even considering a range of supernova timescales, ∼4-8 Myr, due to incomplete sampling of the IMF (Chevance et al. 2022;Stanway & Eldridge 2023), pre-SN feedback should still have performed a significant amount of clearing in these clusters.
Using the calculations presented in Calzetti et al. (2023) adapted to our range of extinctions (A V ∼0.8-4.5 mag), we derive external pressures in the range P ext /k B =(0.1-2.3)×10 6K cm −3 , which are significantly smaller, by about an order of magnitude, than the pressure supplied by photoionization, direct radiation and stellar winds from the star clusters (Calzetti et al. 2023).Thus, at face value, the host environment of our clusters does not represent an impediment to the clusters' ability to remove gas and dust from their immediate surroundings.
However, the Paβ-emitting clusters, albeit selected to be >3,000 M ⊙ , are a factor ∼7 lower mass than the LEGUS line-emitting clusters, with a median mass of 7,000 M ⊙ .Therefore, we can expect stochastic sampling of the stellar IMF to have an impact on the cluster's properties, due to the effect this has on the number and mass range of the massive stars contained in the clusters (Cerviño et al. 2002;Stanway & Eldridge 2023).Stanway & Eldridge (2023), using the binary population models BPASS, find that a 7,000 M ⊙ cluster typically produces stars with maximum mass around 18 M ⊙ , with a 16th-84th percentile range of ∼9-25 M ⊙ .These are mid-O/early-B type stars, which, although still able to produce ionizing photons and stellar winds, generate rates that are an order of magnitude or more below those of a fully sampled IMF (Leitherer et al. 1999;Vink & Sander 2021).Chevance et al. (2022) shows that, in models, stochastic sampling plays a role in the spread of expected delay times for supernovae to occur (see their Figure 4), covering the range ∼4-10 Myr for a 3,000 M ⊙ cluster, but settling at ∼4 Myr for clusters heavier than ∼2×10 4 M ⊙ .Stochastic sampling can also affect our age determinations, since we use deterministic stellar population models to derive our ages, masses and extinctions.However, deterministic models generally produce younger ages than stochastic models (Krumholz et al. 2015), thus our clusters may well span a wider range of ages than those we derive, but will also tend to be older.
Since external pressure in NGC 4449 is about an order a magnitude too low to be a factor in determining the clearing timescale, the picture that emerges from these arguments is that pre-SN and supernova feedback are not effective on the same timescale in all star clusters, but will depend on the cluster's mass.The clearing timescale of 5-6 Myr we find is aligned with expectations from stochastic sampling of the IMF in the mass range of our line-emitting star clusters.

SUMMARY AND CONCLUSIONS
We have utilized multi-wavelength HST photometry, including narrowband filters corresponding to Hα and Paβ hydrogen recombination lines, to build a sample of 99 young cluster candidates in the nearby dwarf-starburst NGC 4449.Leveraging a JWST/NIRCam Paα image as well as our own images, we visually confirmed that these sources possess compact morphologies with coincident line emission.We also found that none of the sources identified by Calzetti et al. (2023) shows line emission in the JWST Paα image.After deriving the SEDs of the candidates using aperture photometry and fitting them to synthetic SED models, we were able to derive the ages, masses, and extinctions for all 99 candidates.We isolate a final catalog of 15 sources with χ 2 red <10 and mass >3000 M ⊙ to minimize the effects of stochastic sampling of the stellar IMF.Our analysis of the resulting sources reveals a notable plume of sources at 5-6 Myr.We also observe relatively few very young sources, younger than ∼4 Myr, and only one with mass≥3,000 M ⊙ .A critical evaluation of our selection method indicates that we are unlikely to have missed from our sample very young (≲4 Myr) sources with compact line emission.
A comparison with the UV-optically selected young star clusters from the LEGUS sample shows that, while our clusters retain their gas up to an age of at least 5-6 Myr, possibly longer, the LEGUS clusters lose their gas by 4 Myr age.The main difference between the two samples is that the young LEGUS clusters extend to masses much higher than those in our sample.Considering only clusters with mass ≥3,000 M ⊙ , the median mass for the LEGUS clusters is ∼3.5×10 4 M ⊙ , while for our sample it is ∼7×10 3 M ⊙ .These results suggest that pre-SN and SN feedback are effective in massive star clusters, which, as a result, do not retain their gas (and dust) beyond ∼4 Myr.Conversely, they appear to be less effective and require longer timescales in the lower-mass star clusters (i.e., those in our sample), and these clusters still retain gas and dust up to the maximum age our SED fits can reliably determine (≲7 Myr).We suggest that the progenitors of these lower-mass clusters are still sufficiently buried in their natal cloud to have escaped detection in our sample.Only after about 5-6 Myr do these dusty clusters begin to emerge from their natal clouds.
The scenario implied by our data is that the clearing timescales due to pre-SN and supernova feedback are cluster mass-dependent; the timescales are not universally short but increase for decreasing star cluster mass.This is likely due to both pre-SN and supernova feedback being impacted by stochastic sampling of the stellar IMF, whose effect increases for decreasing mass.If correct, this scenario would have important repercussions on a galaxy's ability to leak ionizing photons, since the steep cluster mass function implies that galaxies contain many low-mass clusters and relatively few massive star clusters, at least at low redshift (e.g.Adamo et al. 2017;Messa et al. 2018;Cook et al. 2019).Observations of nearby galaxies with the JWST possess improved sensitivity and angular resolution in the NIR and will likely provide an improved sample of star clusters covering the full range of ages, masses, and extinctions in the local universe, enabling accurate determinations of the clusters' clearing timescales.

$JH0\U
ORJ0>0 @ $JH0\U ORJ0>0 @ Figure 6.Mass limit (black solid line) as a function of age, as determined from the I-band detection limit of our images.We set a 3000 M⊙ lower limit (black dashed line) to avoid impact on our physical parameters from stochastic sampling of the IMF (see Section 6).(Left) Mass versus age for all 99 fitted sources.(Right) The same plot as to the left, only for the sources with χ 2 red <10.In the right panel, we also show the location of the LEGUS star clusters with χ 2 red <10, mass≥3,000 M⊙ and coincident with ionized gas emission (magenta triangles).See Adamo et al. (2017) for a discussion of the LEGUS star clusters.red < 10 and mass >3000 M⊙, and with the inclusion of the single visually identified source with comparable characteristics (black triangle).This is a relatively unextincted source with E(B-V)=0.28,suggesting that the approach of visual identification is insufficient at returning a sample of young, highly extincted clusters.(Right) The same clusters as in the left panel, with the addition of the LEGUS clusters with χ 2 red < 10, mass >3000 M⊙ and coincident with ionized gas emission (magenta triangles with error bars).See Adamo et al. (2017) for a discussion of the LEGUS star clusters.

Figure 1 .
Figure 1.A three color composite of NGC 4449 (blue=F336W, green=F814W, red=Paβ emission).The locations of the 99 cluster candidates are marked by yellow circles.Due to the requirement that sources be detected in Paβ, all sources lie within the field of view of the F128N filter, traced by the white line.North is up, East is left.

Figure 4 .
Figure 4. Three color composites of NGC 4449 showing each stage of the sample selection, with the sources marked with yellow circles.(Upper Left) Full set of successfully fitted sources (99 sources).(Upper Right) Successfully fitted sources with χ 2 red < 10 (80 sources).(Lower Left) Successfully fitted sources with M>3000 M⊙ (22 sources).(Lower Right) Successfully fitted sources with both χ 2 red < 10 and M>3000 M⊙ (15 sources).

Table 1 .
Data Sources

Table 2 .
Filter Table HST Filter Pivot λ FWHM Standard Filter Figure 2. Minimum χ 2red versus age of 99 fitted sources.The black dashed line denotes a χ 2 red < 10 cutoff below which there are 80 sources.The 1σ uncertainties in age are calculated from a 1σ interval from the minimum χ 2 red .

Table 3 .
Source Location and Photometry