Characterizing Dust Extinction and Spatially Resolved Paschen-$\alpha$ Emission within 97 Galaxies at $1<z<1.6$ with JWST NIRCam Slitless Spectroscopy

For source identification, we start with a catalog generated in Morishita et al. (2024b). The work retrieved the fully processed images and spectroscopic catalogs made available by the JADES team (Rieke et al. 2023). A photometric catalog was then constructed following the method described in Morishita & Stiavelli (2023), using borgpipe (Morishita, 2021). In our following analysis, we adopt the aperture flux from the catalog, which was corrected to represent the total fluxes of the galaxies. For further details regarding the catalog, the readers are referred to Morishita et al. (2024b).

In this analysis, we utilize publicly available NIRCam WFSS data from the First Reionization Epoch Spectroscopic Complete Survey (FRESCO; GO-1895; PI: P. Oesch; Oesch et al. 2023). FRESCO is a JWST cycle-1 program which features comprehensive deep NIRCam imaging and slitless spectroscopic observations using the F444W filter. The obtained spectra achieve a resolution of R $\sim$ 1600 and covers a wavelength range from 3.8 to 5.0 $\mu$m, which enables us to study multiple hydrogen recombination lines across a wide range of cosmic time with a spatially resolved view.

We follow the methodology outlined in Sun et al. (2023)¹¹1https://github.com/fengwusun/nircam_grism/ and reduce the WFSS data using a combination of the official JWST pipeline and several customized steps as detailed below. We retrieve the stage-1 products (rate.fits) from the MAST archive and assign the world coordinate system (WCS) for each frame. Our background subtraction process involves two steps: we first subtract a median background that is separately generated for each module and pupil. To further enhance data quality and eliminate any potential residual background, we employ additional subtraction using SExtractor (Bertin & Arnouts, 1996). To isolate the emission line from the continuum, we follow the method described in Kashino et al. (2023). This involves using a median filter technique that models the continuum with a sliding window, or “kernel”. By subtracting this modeled continuum from the original data, we effectively extract the emission line image.

With the fully calibrated images, we then measure the astrometry offset between the short wavelength (SW) direct images and the reduced F444W mosaic to align each grism exposure with the direct image, this step is essential to allow the spectral tracing models work properly. The spectral tracing, grism dispersion and sensitivity models are produced from multiple JWST/NIRCam Cycle-1 programs (PID: 1076/1536/1537/1538), the details are described in Sun et al. (2023). We extract the 1d spectra with optimal extraction (Horne, 1986), and fit the emission line flux with a single Gaussian profile.

We utilize the archival HST grism observations from the CANDELS Lyman Alpha Emission At Reionization (CLEAR) survey (PI: Casey Papovich), as detailed in Simons et al. (2023), which provided 12-orbit depth observations with the HST/WFC3 using the G102 grism, covering 12 fields across both the GOODS-N and GOODS-S. Simons et al. (2023) expanded the original CLEAR survey dataset by incorporating observations from additional HST projects that utilize the G102 and G141 grisms overlapping with the CLEAR footprint. This approach achieved complete spectral coverage from 0.8 to 1.7 $\mu$m throughout the observed areas. Leveraging this rich dataset, Simons et al. (2023) have compiled two spectroscopic catalogs that present emission line fluxes and redshifts, derived from a combination of photometry and grism spectroscopy.

In our analysis, we use the redshifts and combined ${\rm H\alpha}$ + [N ii] flux from the CLEAR catalog. To correct for [N ii] contamination, we apply the [N ii]/${\rm H\alpha}$ ratios as a function of stellar mass at $z\sim 1.5$, as calibrated by the $\text{KMOS}^{3D}$ program (Wuyts et al., 2016). For an in-depth understanding of the CLEAR survey’s methodology, including survey design, data reduction processes, and the galaxy catalog compilation, we refer readers to Simons et al. (2023).

We start with the NIRCam + HST catalog as described in Sec. 2.1 and select sources for further analysis as follows:

• •

  First, we perform a cross-match between the NIRCam+HST catalog (Sec. 2.1) with the CLEAR spectroscopic catalog (Sec. 2.3) using a matching radius of $0.\!^{\prime\prime}5$, and retain only the galaxies with available spectroscopic redshift from the CLEAR survey.

• •

  We further refine our selection to include only those galaxies within the redshift range of $1<z<1.6$, which aligns with the coverage of the NIRCam WFSS F444W filter for capturing the ${\rm Pa\alpha}$ line and the HST G141 grism for capturing the ${\rm H\alpha}$ line. The collection of galaxies that pass this phase of selection is referred to as our parent sample.

• •

  Additionally, selected galaxies must have an AB magnitude brighter than 25 in F444W and both ${\rm H\alpha}$ and ${\rm Pa\alpha}$ detections with a SNR greater than 5.

• •

  Lastly, we conduct spectral energy distribution (SED) fitting for the selected galaxies in the previous step. To ensure reliable results, we further restrict our sample by including only those with the SED fitting results with a reduced $\chi^{2}_{\nu}$ value of less than 3.

Our final sample consists of 42 galaxies from GOODS-S and 55 from GOODS-N. The distribution of F444W magnitudes for both our parent and final samples is displayed in Fig. 1.

We conduct SED fitting using the photometric data to estimate stellar masses. We use Code Investigating GALaxy Emission (CIGALE, Boquien et al. 2019). CIGALE is particularly well-suited for analyzing data across a broad wavelength range, extending from X-ray to radio frequencies. Our SED fitting employs a delayed exponential model for star-formation history (sfhdelayed), mathematically represented as ${\rm SFR}(t)\propto t\exp(-t/\tau)$. For the stellar population synthesis, we utilize the model from Bruzual & Charlot (2003), assuming solar metallicity. The model incorporates dust attenuation based on the Calzetti et al. (2000) extinction law. A comprehensive description of the parameters used in our SED fitting can be found in Liu et al. (2023).

In this work, we adopt Case B recombination assumption with ${\rm T_{e}=10^{4}}$ K and ${\rm n_{e}=10^{2}}$ ${\rm cm^{-3}}$ as determined by photoionization models using CLOUDY version 17.02 (Ferland et al., 2017). Under these conditions, the model predicts intrinsic intensity ratios of ${\rm H\alpha}$/${\rm H\beta}$ = 2.79 and ${\rm Pa\alpha}$/${\rm H\beta}$ = 0.305. By applying the Calzetti et al. (2000) reddening curve ($R_{\rm{v}}$ = 4.05) and comparing the observed ${\rm Pa\alpha}$ to ${\rm H\alpha}$ flux ratios with Case B predictions, we calculate $E(B-V)_{\text{neb}}$ for our sample:

Here, $k(\lambda)$ represents the extinction coefficient at wavelength $\lambda$, as prescribed by Calzetti et al. (2000). From the calculated color excess, we derive the extinction values A(${\rm H\alpha}$) and A(${\rm Pa\alpha}$) using:

Subsequently, with the estimated dust extinction from the ${\rm Pa\alpha}$ and ${\rm H\alpha}$ line flux ratios, we compute the intrinsic SFRs. The selection of the parent sample in this study relies on the ${\rm Pa\alpha}$ detection, which is sensitive to HII regions that signal active star formation. Thus, we expect the majority of our selected sources to be star-forming galaxies. To assess their star-forming activities more accurately, we derive dust-corrected SFRs using ${\rm Pa\alpha}$ line flux and color excess $E(B-V)_{\text{neb}}$. The ${\rm Pa\alpha}$ flux is determined by fitting a single Gaussian profile to continuum-subtracted 1D spectra (Sec. 2.2). Upon determining the intrinsic ${\rm Pa\alpha}$ flux using the Calzetti et al. (2000) approach, we convert it to the ${\rm H\alpha}$ flux based on the line ratio predicted under Case B recombination. SFRs are then estimated from the ${\rm H\alpha}$ luminosity following the SFR-$L_{\mathrm{H\alpha}}$ relation established by Kennicutt (1998):

To align with the Chabrier (2003) IMF, we scale the SFRs by dividing them by a factor of 1.53 (Driver et al., 2013). As illustrated in Fig. 2, most of our sample galaxies are located on or near the star formation main sequence predicted by Whitaker et al. (2014). Our analysis indicates that the ${\rm Pa\alpha}$ line is a robust indicator of SFRs, relatively unaffected by dust extinction and offers a more direct measure of the ionizing radiation from young O and B stars. This is particularly true for less massive galaxies with a stellar mass below $\rm{10^{10}~{}M_{\odot}}$, which are less affected by dust extinction, as demonstrated in the comparison between the left and right panels of Fig. 2.

In our redshift range of interest, the F444W grism captures the ${\rm Pa\alpha}$ emission line, a robust tracer of SFRs, while the F444W image captures the rest-frame K-band image ($\sim$ 2 $\mu$m), providing an accurate measure of the galaxy’s stellar mass unaffected by star formation history (e.g., Kauffmann & Charlot, 1998; Bell & de Jong, 2001; Drory et al., 2004). Therefore, comparing the flux distributions from the F444W image and the F444W grism enables us to study the mechanisms governing star formation and mass assembly within galaxies, with minimal impact from dust extinction.

However, the NIRCam grism achieves a high resolution of R $\sim$ 1600 at 3.95 $\mu$m, which corresponds to a velocity dispersion of $\sigma$ $\sim$ 80 km/s (Nelson et al., 2023). This resolution will lead to morphological distortion of emission-line maps along the spectral axis of galaxies having velocity dispersion greater than 80 km/s. Therefore, the 2D spectra we extract are not purely emission line maps; they also incorporate kinematic properties along the spectral axis. Overcoming such degeneracy is possible with an additional medium-band filter that serves to extract emission line map (e.g., Nelson et al., 2023). Yet, this approach significantly narrows the applicable redshift range due to the coverage of medium filters. Recognizing these constraints, our analysis in this study is thus focused on examining the light profile distribution in the cross-dispersion direction, which remains unaffected by velocity dispersion. Our method involves the following steps:

• •

  Model fitting: We employ the technique outlined in Szomoru et al. (2010), which involves using a Sérsic $R^{1/n}$ profile (Sérsic, 1963) to model the observed galaxy images and 2D spectra. We first select an unsaturated star within the same field of view to represent the PSF. Subsequently, we begin the Sérsic $R^{1/n}$ fitting process using the initial parameters derived from SExtractor. To refine our model, we employ optimization techniques to adjust the parameters. The best-fit PSF-convolved Sérsic model is determined by the minimal residual.

• •

  Intrinsic Profile Reconstruction: Using the previously identified best-fit parameters, we construct a Sérsic profile and incorporate the residuals obtained from comparing the observed data with the PSF-convolved model. This profile, enhanced with the residuals, is considered the intrinsic flux profile for further size measurement analysis. Although the residuals are still affected by the PSF, this technique has proven to accurately reconstruct the true flux distribution, even for galaxies poorly described by a Sérsic profile (Szomoru et al., 2010). Fig. 3 schematically illustrates the reconstruction process.

• •

  Flux Distribution Analysis: With the reconstructed intrinsic profiles for both 2D spectra and images, we compute the one-dimensional cross-dispersion light intensity profile by summing pixel values across each row. Specifically, we begin by creating cutouts of the 2D spectra and images, each with dimensions of $3.\!^{\prime\prime}0\times 3.\!^{\prime\prime}0$; we then rotate F444W images to align with the 2D spectra, ensuring the horizontal direction is set as the dispersion direction. Additionally, we mask out pixels identified in the segmentation maps as belonging to other galaxies and those with negative values from the continuum removal process, thereby effectively minimizing the effects of nearby galaxies and artifacts from bad pixels.

Based on the method above, we obtain flux profiles for both the stellar and star formation components in the cross-dispersion direction. To quantitatively compare the distributions of stellar mass and star formation in individual galaxies, we measure the half-light radius along the cross-dispersion direction. Given the different depths of the grisms and direct images, the SNR of the stellar components is higher than that of the SFR components. For a fair comparison, we define the total flux as the amount captured within a specified radius in the flux profile, with this boundary determined by the point at which the SNR of the 2D spectra drops to 1, thereby accounting for the spectra’s shallower depth. The half-light radius is then calculated by locating the point where the cumulative flux reaches half of the total flux. This point, measured from the center of the cutout along the y-axis, is defined as the half-light radius. Although this diverges from the traditional definition of half-light radius, this method facilitates a direct comparison of the spatial distributions of star formation and stellar mass within galaxies, based on their light profiles in the cross-dispersion direction.

To enhance the SNR and measure the average properties of ${\rm Pa\alpha}$ and F444W profiles, our sample is divided into two sets of bins: one based on stellar mass and another based on specific star formation rate (sSFR). For each categorization, galaxies are classified into low, medium, and high bins, with 32, 32, and 33 galaxies in each respective category. Following the methodology described in Nelson et al. (2016a), we stack the ${\rm Pa\alpha}$ spectra and F444W images in each bin to create the mean profiles. Before proceeding with the stacking, we mask nearby sources using segmentation maps, deconvolve the PSF, and rotate the F444W images to ensure that the horizontal direction corresponds to the dispersion direction, as detailed in 3.3. To prevent the final stacked profile from being dominated by brighter sources, we apply weights based on their F444W flux. According to Nelson et al. (2016a), deprojecting, rotating, or scaling the images does not significantly affect our main conclusions; therefore, we do not apply additional rotation to the dispersion axis-aligned profiles. Furthermore, we calculate the half-light radius in the cross-dispersion direction of the stacked profiles using the same method applied to individual galaxies.

We examine the nebular reddening within our sample galaxies and investigate the relationship between dust extinction and stellar mass. In Sec. 3.2, we derive the $E(B-V)_{\text{neb}}$ using the ratio of the ${\rm Pa\alpha}$ and ${\rm H\alpha}$ emission lines. We then calculate the dust extinction A(${\rm H\alpha}$) and A(${\rm Pa\alpha}$) from these ratios. We present our findings in Fig. 4, with error bars indicating the flux measurement uncertainties.

Our results reveal a positive correlation between stellar mass and dust extinction, as inferred from ${\rm Pa\alpha}$-to-${\rm H\alpha}$ line flux ratios. This is consistent with previous studies that employed Balmer-line ratios to assess dust extinction, as well as those using far-IR and UV continuum measurements (e.g., Garn & Best, 2010; Zahid et al., 2013; Reddy et al., 2015; Fudamoto et al., 2020; Shapley et al., 2023). Based on the model predictions by Calzetti et al. (2000), A(${\rm Pa\alpha}$) is significantly smaller than A(${\rm H\alpha}$) due to the much smaller k($\lambda$). Our results indicate that at our target redshift, A(${\rm Pa\alpha}$) is typically around $\sim$ 0.1 mag, with the median value increasing to only $\sim$ 0.2 mag in the most massive bin. While ${\rm Pa\alpha}$ is intrinsically fainter than ${\rm H\alpha}$ (I_Paα/I_Hα = 0.11 under Case B), this results demonstrate the robustness of ${\rm Pa\alpha}$ as a star formation indicator, less attenuated by dust extinction and more sensitive to the intrinsic SFR. Such conclusions are broadly in agreement with Reddy et al. (2023), who investigated the Paschen lines of 63 galaxies at redshifts $z=1.0$–$3.1$ using NIRSpec. We note that a small fraction of galaxies in our sample exhibit a negative extinction. This happens when the observed ratio of ${\rm Pa\alpha}$ to ${\rm H\alpha}$ is less than the dust-free minimum value of 0.109, which is not physical. Such observations could be attributed to the relative faintness of the ${\rm Pa\alpha}$ line in comparison to the ${\rm H\alpha}$ line, possibly leading to incomplete detection of the galaxy’s fainter regions in the 2D spectra. Furthermore, our reliance on empirical calibration for ${\rm H\alpha}$ flux, based on the [N ii]/${\rm H\alpha}$ ratio, introduces additional uncertainties into the ${\rm H\alpha}$ flux determination.

Thanks to the HST grism, previous studies (e.g., Nelson et al., 2016a; Vulcani et al., 2016; Nelson et al., 2019; Matharu et al., 2022) have significantly advanced our understanding of star formation regions within galaxies through spatially resolved ${\rm H\alpha}$ emission line maps. As one of the most prominent hydrogen recombination lines, ${\rm H\alpha}$ is closely linked to ionized regions and serves as a reliable tracer of star formation. However, even ${\rm H\alpha}$ is susceptible to significant dust extinction. This presents a challenge for studying star formation locations as the dust amount may differ significantly across various regions of a galaxy (Nelson et al., 2016b; Tacchella et al., 2018; Matharu et al., 2023). Such variability introduces notable uncertainties in quantifying galaxy sizes and analyzing morphological features. The utilization of spatially resolved ${\rm Pa\alpha}$ emission line map offers a solution to overcome these challenges. In this section, we present the galaxy sizes derived from ${\rm Pa\alpha}$ 2D spectra and F444W continuum data, offering insights into the spatial distributions of star formation activities and existing stellar mass within galaxies. In order to better understand the observed size difference in our galaxies, we investigate the relationship between the SFR-to-stellar size ratio, sSFR, and stellar mass.

In Sec. 4.2, we investigated the half-light radius of ${\rm Pa\alpha}$ and F444W, and their potential correlation with stellar mass. Our main findings are illustrated in Fig. 5 and 8, which show that galaxies, especially those with high stellar masses, exhibit ${\rm Pa\alpha}$ extent notably larger than their stellar continuum. These observations align with the findings of Nelson et al. (2016a), who reported a larger effective radius for ${\rm H\alpha}$ compared to F140W in a study of 3200 galaxies at $0.7<z<1.5$. Notably, the size difference was more pronounced with increasing stellar mass. While Nelson et al. (2016a) observed a steeper size evolution with stellar mass than our results suggest, this difference could be attributed to the different tracers used, i.e., ${\rm H\alpha}$ and F140W in their study versus ${\rm Pa\alpha}$ and F444W in ours. Given that our tracers are expected to be more sensitive to attenuated emissions, the difference indicates that galaxy central regions of our galaxies systematically host more attenuated emissions.

It should be noted that the half-light radius, determined from the integrated profile along the cross-dispersion axis, may not directly reflect the overall size of the galaxy due to the significant impacts of orientation and ellipticity. This factor might also contribute to the observed scatter in the relationship between size and stellar mass as shown in Fig. 5. Nevertheless, despite these limitations, the analysis of the ratio of galaxy sizes remains a robust metric.

The observed increase in galaxy size ratio with stellar mass necessitates further investigation into the underlying physical processes that govern galaxy structural formation. It remains a challenging question to understand how galaxies form their central mass concentrations, particularly massive bulges, and evolve into mature, passive systems, given the presence of active star formation regions surrounded by dust. Our study provides a direct probe of star-forming regions, minimizing reliance on assumptions required for dust extinction corrections, such as UV slope and SED fitting techniques. Consequently, we anticipate that the regions traced by ${\rm Pa\alpha}$ will predominantly represent ionized hydrogen areas within the galaxies. Based on this foundation, and the positive correlation between ${\rm Pa\alpha}$-to-stellar size ratio and increasing stellar mass (Fig. 5), our finding that most galaxies exhibit larger ${\rm Pa\alpha}$ sizes supports the inside-out growth model of galaxies (e.g., van Dokkum et al., 2010). This model suggests that galaxies gradually build up their outer regions while the inner regions have completed their evolution more rapidly.

Our study distinguishes itself from previous research that utilized ${\rm H\alpha}$, which faces significant uncertainties due to dust-obscured central regions. Our research confirms the scenario of extended star-forming regions, employing an independent tracer made feasible for the first time by NIRCam WFSS. Although ${\rm Pa\alpha}$ may still be significantly affected by dust extinction in extremely dusty starburst galaxies, as inferred from a recent study by Bik et al. (2023) that studied a galaxy with A_V = 17.2, our sample primarily consists of normal star-forming galaxies, most of which have A(${\rm Pa\alpha}$) $<$ 0.2, with the dustiest having A(${\rm Pa\alpha}$) $\sim$ 0.6.

As discussed in Sec. 5.1, we observe considerable scatter in the measurements of galaxy size ratios, particularly among medium-mass to massive galaxies (i.e. galaxies with stellar mass $\rm{M_{*}>10^{9.5}\,M_{\odot}}$). Within our sample, we identify three distinct flux distributions in the emission line maps of these galaxies. Even among galaxies with similar stellar masses and SFRs, their ${\rm Pa\alpha}$ flux profile can vary significantly. For example, as illustrated in Fig. 9, the galaxy on the left displays extended and clumpy star-forming regions, whereas the middle one shows centrally concentrated star-forming regions; the right galaxy, however, exhibits expanded star-forming regions, suggesting inside-out growth. The diverse star formation regions observed in our massive galaxies are interesting. The overall trend revealed by stacked profiles in Fig. 8 shows that massive galaxies exhibit larger star formation regions, in good agreement with findings from previous studies. However, the larger scatter in the size ratio observed in individual massive galaxies (right panel of Fig. 5), which originates from galaxies with different star formation patterns, warrants further investigation. Swinbank et al. (2012) studied nine ${\rm H\alpha}$-selected galaxies at $0.8<z<2.2$ using SINFONI/VLT and reported a variety of ${\rm H\alpha}$ morphologies; at least six of these galaxies exhibited extended or clumpy features, while two displayed compact features. In a related study, Förster Schreiber et al. (2018) analyzed 35 star-forming galaxies and found that the majority were late-type, disk-dominated systems, with eight displaying prominent rings or off-center clumps. The diversity observed in ${\rm H\alpha}$ is consistent with the variety of morphologies we have seen in ${\rm Pa\alpha}$. In this section, we delve into the underlying physics behind the formation of clumpy, extended, and compact star formation regions in these galaxies.

Both the clumps distributed at the outskirts of galaxies and the extended ${\rm Pa\alpha}$ can lead to a higher effective radius measured with ${\rm Pa\alpha}$ compared to that of the stellar continuum. Gaseous clumps, often contributing to the formation of bulges and disks, are regarded as key drivers of galaxy structure formation and evolution and are commonly observed in star-forming galaxies at $z>1$ (e.g., Dekel et al., 2009a, b). The existence of massive clumpy galaxies can be attributed to dynamically unstable gas, which leads to the formation of giant star-forming clumps within the disk (Krumholz et al., 2018), or it could be indicative of gas-rich minor merger events (Agertz et al., 2009).

To further understand the formation mechanisms of these complex structures, observations of both ionized and molecular gas are essential. There are only a handful of studies focusing on clumpy galaxies; recently, Bik et al. (2023) analyzed the dusty, star-forming galaxy GN20 at $z\sim 4$ using MIRI/MRS, detecting spatially resolved ${\rm Pa\alpha}$. Their emission line map revealed that ${\rm Pa\alpha}$ emissions are concentrated in four distinct clumps. Additionally, GN20 was observed with ALMA (Hodge et al., 2012, 2016), where the CO(2-1) emission also revealed a similar clumpy structure. These findings suggest that the galaxy is in the late stages of a major merger, implying that the clumps in the gas-rich disk are a result of this merger, while the central starburst is directly driven by the event. Similarly, Cochrane et al. (2021) reported bright, extended structures in a $z=2.2$ galaxy, with clear clumps in ${\rm H\alpha}$ and extended dust continuum emissions. This disordered morphology and extreme star formation are linked to a merger event. Although our sample primarily consists of main sequence galaxies, it is possible that some of the clumpy structures we observed are triggered by merger activities. Such sub-structures remain an active area of exploration, and our study indicates that clumpy star-forming populations can be efficiently studied with WFSS. To fully comprehend the impact of dust and molecular gas in shaping the characteristics of HII regions, further investigations using sub-mm observations are crucial. The extended star formation region observed in the middle panel of Fig. 9 has been well-studied in previous ${\rm H\alpha}$ observations utilizing HST grism and adaptive optics (e.g., Nelson et al., 2016a, b; Vulcani et al., 2016; Matharu et al., 2022). Unlike rest-frame UV-optical light, ${\rm H\alpha}$ emission suffers less attenuation by dust and is closely linked to the star formation activities of massive stars. Given the relatively short lifespans of massive stars, their observed population offers a direct measure of the current SFR (e.g., Kennicutt, 1998; Pflamm-Altenburg et al., 2007). It was reported that when galaxies assemble their massive structures, star formation typically starts in the central region and then spreads outward. This process helps build up the galaxy’s outskirts, leading to the larger size of star formation regions observed in massive galaxies.

Compact star formation presents another intriguing scenario for further investigation. One hypothesis for the most massive galaxy populations is that they undergo a process known as “compaction”, which leads to the evolution of massive bulges over cosmic time (Dekel & Burkert, 2014; Zolotov et al., 2015; Lacerda et al., 2020; Marques-Chaves et al., 2022). Subsequent starbursts and/or AGN feedback may then drive these galaxies toward quiescence (Genzel et al., 2014; Yesuf et al., 2014), resulting in quenched, dense cores (Damjanov et al., 2009; van Dokkum et al., 2009; Szomoru et al., 2012; Morishita & Ichikawa, 2016). Recent studies have elucidated the formation mechanisms of compact galaxies. Franco et al. (2020) highlighted the significant role of compact galaxies at z $\sim$ 2.5 - 3, associating them with short depletion timescales. Simulations such as Illustris suggest that these galaxies often form early in extreme environments characterized by an abundance of gas and/or in (post-)merging systems where large inflow of gas trigger intense starbursts (Wellons et al., 2015). The acquisition of external gas is believe to be responsible for initiating secondary star formation within the core regions of these galaxies that result in their compact morphology (Chen et al., 2016). Importantly, such populations might be underrepresented or even missed in previous ${\rm H\alpha}$ analyses, at redshifts $z=1-2$, our target epoch corresponding to the peak of cosmic star formation activity, dust extinction affecting ${\rm H\alpha}$ can vary significantly, with magnitudes ranging from $\sim 1$, to $>3$ (e.g., Förster Schreiber et al., 2009; Sobral et al., 2012; Kashino et al., 2013; Koyama et al., 2019). Moreover, dust extinction can vary widely across a galaxy. As discussed in Matharu et al. (2023), dust attenuation is not a constant value within galaxies; instead, the spatial distribution of dust attenuation can vary across different stellar masses. A recent study by Cochrane et al. (2021) investigated a massive galaxy at $z=2.2$ using multi-wavelength data and found A(${\rm H\alpha}$) values ranging from $\sim 2-3$ in the outskirts to $>$ 5 in the central region. Interestingly, their ${\rm H\alpha}$ emission line map, obtained using SINFONI/VLT, showed a clear offset from the peak of the dust continuum emission traced by ALMA. Tadaki et al. (2014) also reported color gradients in galactic clumps, with central clumps being redder, indicating that dusty star-formation activity is concentrated in these nuclear regions. Such dust gradients introduce significant complexity to the analysis of ${\rm H\alpha}$-traced HII regions, potentially biasing detections toward galaxies with extended star formation regions. Galaxies featuring compact, yet dusty, star formation regions might remain undetected in ${\rm H\alpha}$ observations. Therefore, compact star-forming galaxies might still be understudied in previous ${\rm H\alpha}$-based work. Future surveys using better tracers less affected by dust, such as ${\rm Pa\alpha}$, will be crucial to provide a comprehensive view of the different paths of galaxy structure formation.

One might question our size measurements due to the differing depths of the F444W imaging and WFSS spectroscopy, with the former being significantly deeper. Such a depth discrepancy could result in an underestimation of the size of star formation regions, as the fainter outskirts of galaxies may remain undetected. Ideally, a comparison between the continuum captured by WFSS and the observed emission lines would help to mitigate these discrepancies. However, practical implementation of this method faces significant challenges as the continuum captured in dispersed spectra is frequently contaminated by nearby continua, thereby complicating the accuracy of measurements. This issue is particularly pronounced for larger galaxies, where the continuum is more likely to overlap with other sources. Therefore, our analysis relies on direct images from F444W. To further assess the potential effects of this depth difference, we carefully inspect each source and select those without significant contamination from nearby galaxies. We measure their sizes using the modeled continuum (see Sec. 2.2 for the continuum modeling process), applying the same measurement method as used for F444W images and 2D spectra. The results, illustrated in Fig. 10, show that measurements based on the modeled continuum and F444W are generally consistent. Therefore, we assert that the depth difference does not significantly impact our conclusions. Even if the sizes of the 2D spectra are underestimated, our findings that the star formation regions traced by 2D spectra are more extended remain unaffected.

The inherent complexities in galaxy morphology and kinematics, particularly the elongation of morphologies along the spectral axis due to kinematic distortions, prevent us from obtaining accurate and separate information on spatial and spectral information. As discussed in Nelson et al. (2023), for high EW emission lines, such as ${\rm H\alpha}$, a combination of medium-band filters and WFSS can provide insights into both the velocity and morphology of the emission line maps. However, this approach may be challenging for weak lines or for fields without medium-band filters on the targeted line. To apply this method to our research objectives, specifically using ${\rm Pa\alpha}$ to trace dust-free star formation, a larger sample size would be required. A future joint analysis involving multiple programs with available medium-band and WFSS observations will offer more opportunities to study galaxies with strong ${\rm Pa\alpha}$ emission lines.