Finite Resolution Deconvolution of Multi-Wavelength Imaging of 20,000 Galaxies in the COSMOS Field: The Evolution of Clumpy Galaxies Over Cosmic Time

Compact star-forming clumps observed in distant galaxies are often suggested to play a crucial role in galaxy assembly. In this paper, we use a novel approach of applying finite resolution deconvolution on ground-based images of the COSMOS field to resolve 20,185 star-forming galaxies (SFG) at 0.5<z<2 to an angular resolution of 0.3", and study their clumpy fractions. A comparison between the deconvolved and HST images across four different filters shows good agreement and validates the deconvolution. We model spectral energy distributions using the deconvolved 14-band images to provide resolved surface brightness and stellar mass density maps for these galaxies. We find that the fraction of clumpy galaxies decreases with increasing stellar masses, and with increasing redshift: from ~30% at z ~ 0.7 to ~50% at z ~ 1.7. Using abundance matching, we also trace the progenitors for galaxies at z ~ 0.7 and measure the fractional mass contribution of clumps toward their total mass budget. Clumps are observed to have a higher fractional mass contribution toward galaxies at higher redshift: increasing from ~1% at z ~ 0.7 to ~5% at z ~ 1.7. Finally, the majority of clumpy SFGs have higher specific star formation rates (sSFR) compared to the average SFGs at fixed stellar mass. We discuss the implication of this result to in-situ clump formation due to disk instability.


INTRODUCTION
Over the last 3 decades, we have developed a broad understanding of the stellar mass build-up within galaxies over much of cosmic time. The star formation rate density is known to decrease over the past ten billion years (e.g., Lilly et al. 1996;Madau et al. 1996;Madau & Dickinson 2014), while the bulk of the global stellar mass density was formed during the epoch of peak cosmic star formation at 1 < z < 3 (e.g., Dickinson et al. 2003;McLeod et al. 2021). However, a complete picture of where and how stars are formed within galaxies remains elusive as it is still Corresponding author: Visal Sok sokvisal@yorku.ca unclear which mechanisms are driving the stellar mass assembly.
When looking at galaxy populations across cosmic time, we find that there is a transition in galaxy morphologies toward higher redshifts. In particular, peculiar galaxies that host kiloparsec-scale and clump-like structures outnumber Hubble-Sequence galaxies at z > 1 (Elmegreen et al. 2007). These clumpy structures are star-forming regions as indicated by their enhanced specific star formation rates (sSFR) (e.g., Wuyts et al. 2012;Hemmati et al. 2014), and their detection in high-resolution observations of Hα emission line maps (e.g., Wuyts et al. 2013;Mieda et al. 2016). While there is no formal definition, compact star-forming structures in distant galaxies are generally dubbed as "clumps", and their host galaxies as "clumpy galaxies". Clumpy star-forming galaxies (SFGs) have been studied for both field and lensed galaxies, with a variety of multi-wavelength data, including the restframe UV and optical imaging, emission line maps, and CO observations (e.g., Elmegreen & Elmegreen 2005;Elmegreen et al. 2007Elmegreen et al. , 2009aFörster Schreiber et al. 2011;Wuyts et al. 2012;Murata et al. 2014;Livermore et al. 2015;Guo et al. 2015;Shibuya et al. 2016;Cava et al. 2018;Zanella et al. 2019). Despite their ubiquity at high-z, the origin and evolution of clumps are not well constrained by observational studies due to the limiting resolution of our observing facilities. Size estimates of star-forming clumps can vary from tens of parsecs to a few kiloparsecs based on observations at z ∼ 1−4 (e.g., Livermore et al. 2015;Soto et al. 2017;Cava et al. 2018). At these scales, even the Hubble Space Telescope (HST) can only marginally resolve these structures at its limiting angular resolution of ∼0.1", or a physical resolution of ∼1 kpc, at z ∼ 2.
As a result, the formation and evolution of clumps can be explained within a wide range of theoretical models. The mechanisms for clump formation are associated with both in-situ and ex-situ origins. Such clumpy and irregular morphologies can be the result of ex-situ processes such as mergers, which is supported by some observational studies (Guo et al. 2015;Calabrò et al. 2019). It is also not clear whether clumps are just remnants of accreted satellites that have not been tidally disrupted (e.g., Mandelker et al. 2014;Zanella et al. 2019). The other line of evidence for their formation is that these clumps are not only observed in interacting systems, but also galaxies that exhibit a dynamic structure consistent with that of a rotating disk as identified on the basis of their kinematics (e.g., Förster Schreiber et al. 2009;Genzel et al. 2011;Förster Schreiber et al. 2018). A majority of these high-redshift (z ∼ 1−3) disks are often observed to be gaseous and turbulent, with gas mass fraction of 20%-80% (e.g., Daddi et al. 2010;Tacconi et al. 2010) and intrinsic gas velocity dispersion of ∼20−90 km s −1 (e.g., Förster Schreiber et al. 2009;Wisnioski et al. 2015;Förster Schreiber et al. 2018). These two properties can lead to insitu processes where gravitational instabilities create large star-forming clumps: large gas dispersion results in large Jeans lengths, with the high gas fraction sustaining on-going star formation (e.g., Toomre 1964;Dekel et al. 2009).
Due to their unique mode of star formation during the main epoch of massive galaxy building, clumps are often suggested to have a central role in galaxy assembly. However, the details of how clumps affect the growth of galaxies is poorly understood at the moment. For example, the fate of clumps is inconclusive; clumps can be long-lived stellar structures that migrate inward and coalesce onto the progenitor of present-day galactic bulges (e.g., Bournaud 2016, Mandelker et al. 2014, or short-lived phenomena, where strong stellar feedback disperse the gas and unbind the stellar system on time-scale of <100 Myr (Hopkins et al. 2012).
In order to decipher the nature of star-forming clumps at high-z, high resolution and multi-wavelength data are needed to resolve these structures and to infer their physical properties, respectively. At present, HST is the only observatory capable of providing such data. While HST has now been in service for approximately three decades, there are only a few existing fields with truly deep and multi-wavelength imaging in more than a handful of filters (e.g., GOODS; Giavalisco et al. 2004). Data from the GOODS-South and ERS (Windhorst et al. 2011) have already been exploited for resolved stellar populations studies (Wuyts et al. 2012). While recent data such as the Hubble Frontier Field (HFF; Lotz et al. 2017) are available, the amount of future data that will be comparable to these remains small. On the other hand, larger cameras on ground-based observatories can take advantage of the large field of view, which allow for a larger volume of the sky to be surveyed. However, ground-based observatories are subjected to atmospheric turbulence, and the seeing of these images are effectively degraded, rendering them of limited use for resolved studies of galaxies.
Image deconvolution has often been explored as a solution to reverse the atmospheric and instrumental blurring effect. One challenge with deconvolution is that a sampled image cannot be fully deconvolved without violating the sampling theorem (Magain et al. 1998). A solution is to deconvolve images to a finite resolution, such that the resolution of the deconvolved images is compatible with the new sampling interval. Further pioneering work by Cantale et al. (2016a) onto this deconvolution technique has shown that it is possible to achieve close to space-based resolutions in deconvolved ground-based images (see Cantale et al. 2016b). The new algorithm, finite resolution deconvolution (FIREDEC), has two main features that are advantageous compared to previous deconvolution algorithms. (1) FIREDEC deconvolves to a finite resolution, ensuring that solutions are well-sampled and constrained, and (2) FIREDEC introduces a new regularization scheme to characterize noise, allowing the algorithm to better distinguish signal from noise.
In this paper, we present FIREDEC as an alternative method to obtain high resolution images of distant galaxies. We use image deconvolution to provide resolved stellar properties of high-z (0.5 < z < 2) galaxies within the Cosmic Evolution Survey (COSMOS; Scoville et al. 2007) field. The original ground-based images have a typically seeing of 0.6−0.8 arcsec, and are deconvolved to a target resolution of 0.3". This resolution corresponds to a physical size of ∼2.4 kpc at our redshift range. COSMOS is the ideal field to perform this study as it covers an area of ∼2 deg 2 , and has arguably the most impressive array of multi-wavelength coverage (30+ photometric bands).
The paper is structured as follows. Section 2 describes the observational datasets and our selection of high-z starforming galaxies. In Section 3, we give a brief introduction to image deconvolution, and provide a comparison between deconvolved and HST images. In Section 4, we discuss our methodology for constructing spatially-resolved stellar population maps, while Section 5 discusses how the surface mass density and surface brightness maps are used to identify and classify clumpy galaxies. We present the measured fraction of clumpy galaxies in relation to redshift and global properties of galaxies such as stellar masses and star formation rates in Section 6. In Section 7, we trace the progenitors of clumpy SFGs using abundance matching, and present the fractional mass contribution of star-forming clumps as galaxies evolve. Finally, we discuss the correlation between star formation rates and clumpy morphologies, and its implications on the origin of star-forming clumps. We adopt the following cosmological parameters Ω M = 0.3, Ω Λ = 0.7, and H o = 70 km s −1 Mpc −1 . All magnitudes are quoted in the AB magnitude system.

DATA DESCRIPTION
We make use of the COSMOS/UltraVISTA catalog from . The catalog was obtained using PSFmatched photometry in 30 photometric bands that span from GALEX NUV to Spitzer MIPS 24 µm. A total of 262,615 sources were detected within the UltraVISTA K s -band imaging, which reaches a depth of K s,tot < 23.7. The field covers a total area of 1.6 deg 2 . A detailed description of the photometric catalog construction, along with measurements of the photometric redshift, stellar mass, and UV+IR-determined star formation rate, is presented in .

Selection of the Photometric Data
The photometric data from the COSMOS/UltraVISTA catalog consisted of 30 filters with varying seeing. For example, the UV and NIR data have seeing that typical ranges between 0.5−1.2" . Ideally, we would want to deconvolve all the available imaging of COSMOS. However, deconvolving to a higher angular resolution is the equivalent of recovering high frequency components. The deconvolution of data with both poor and good image qualities (IQ) to a similar (higher) resolution introduces signal-to-noise ratio (S/N) variations, as deconvolution introduces high frequency components within its solution. It can therefore become difficult to reconcile structures within deconvolved images when using a non-homogeneous data that range in image qualities, and can bias our measurements of stellar properties from SED modeling.
A solution is to deconvolve photometric data that has good and relatively homogeneous image qualities. Within the large number of filters available in COSMOS, we choose not to use the g + filter of the Subaru/SuprimeCam as the data have the worst seeing compared to the others (1.01" to 1.20"). We also omit the i + -band image as its longer integration time caused most stars to be saturated and unusable for point spread function modeling. Finally, most of the medium-band photometric data have poor image qualities (e.g., median IQ 0.9), and are therefore omitted.
The total number of filters that are used within this study is 14. These consist of four optical broadband images (B j V j r + z + ) and the six IA and IB optical intermediateband images from the Subaru/SuprimeCam. The four nearinfrared images, Y JHK s , from the UltraVISTA survey are also included. Our selected data have a total seeing range of 0.5−0.9", as summarized in Table 1.

Selection of Galaxy Sample
The COSMOS/UltraVISTA catalog contained 262,615 sources. The parent galaxy sample is obtained using the flag, "use=1". Essentially, the "use=1" flag excludes objects that are flagged as stars based on a color-color cut, objects that are badly contaminated by nearby bright stars, and objects that contain bad pixel regions. The total number of galaxies in the parent sample is 160,070.
Using the derived redshift, stellar mass, and star formation rate from the COSMOS/UltraVISTA catalog, we further filter suitable galaxies for the study. Specifically, we select star-forming galaxies at 0.5 < z < 2, as this paper is focused on the evolution of clumpy star-forming galaxies. Also, only galaxies that satisfy a target S/N cut are included in the study. The choice for the target S/N is discussed in Section 2.2.1 as well as the criteria needed to ensure that deconvolution can recover new spatially-resolved information within the galaxies.

Limiting S/N for Image Deconvolution
Image deconvolution is an iterative process where the original image is resampled and pixel intensities are redistributed until the solution minimizes a defined cost function (see Section 3.1). As a result, it is of limited value to deconvolve objects that are already faint and have low S/N, as deconvolution cannot recover fine scale structures that are not already at high S/N in the original images.
The recovered S/N per pixel by deconvolution depends on the surface brightness of the galaxy. In the ideal situation, we would know the shape and size of galaxies a priori and can select galaxies based on its surface brightness to ensure that we achieve some minimum S/N per pixel. In practice, such morphological parameters are not known from ground-based images. Therefore, we use the integrated S/N of the galaxies in the K s -band as a proxy for what their surface brightness may be. The integrated S/N is calculated as the ratio between the total flux and error as determined by SExtractor's flux auto. For galaxies between 0.5 < z < 2, the K s -band is the best single-filter proxy for stellar mass as it probes the rest-frame NIR for galaxies at z ∼ 0.5 and the rest-frame optical at z ∼ 2. Our tests shows that the amount of information recovered from deconvolved galaxies with S/N < 20 is limited, so deconvolution is only applied to galaxies with S/N > 20 in the K s -band. Figure 1 shows the distribution of galaxies based on the integrated K s -band S/N and redshift for different mass bins. The horizontal dotted line shows the S/N cut we use, and galaxies below the dotted line are excluded in the analysis in the remainder of the paper. We also omit galaxies with log(M * /M ) < 9.8 as the majority of them have an integrated S/N below the limiting S/N. Similarly, we exclude all galaxies above z > 2 as only massive galaxies with log(M * /M ) > 11 satisfy the required S/N. The total number of galaxies from the parent sample that satisfies the S/N, mass and redshift cuts is 32,903.
Note that our sample is a S/N limited sample, not a masslimited sample. Using the integrated S/N as a proxy assumes that all galaxies at a given magnitude have similar light profiles, which is not necessarily true. This is however unavoidable as the recovered S/N across the deconvolved images will depend on the surface brightness of the galaxy. In order to obtain a uniform, mass-complete sample with redshift, one would need to model how the surface brightness depends on mass and redshift. However, based on our selection crite-0.5 1.0 1.5 2.0 2.5 z 0 50 100 150 200 250 300 K s Signal-to-Noise Ratio 9.4 < log(M * /M ) < 9.8 9.8 < log(M * /M ) < 10.2 10.2 < log(M * /M ) < 10.6 10.6 < log(M * /M ) < 11.0 11.0 < log(M * /M ) Figure 1. The Ks signal-to-noise ratio as a function of redshift. Galaxies are binned into different redshift and stellar mass bins. Only galaxies between the redshift of 0.5 < z < 2 and galaxies with S/N (Ks) greater than 20 are selected, as denoted by the dotted lines. Top: The distribution of all the galaxies in COSMOS at 0.5 < z < 2 within the UVJ diagrams excluding those with log(M * /M ) < 9.8 and Ks-band S/N < 20. Bottom: The distribution of the selected star-forming galaxies in the SFR-M plane. ria, only low-mass galaxies with low surface brightness are excluded at higher redshifts, and this should cause only minimal bias in the results.

Defining Star-Forming Galaxies
The main focus of the study is on the stellar populations of galaxies during their lifetime prior to quenching, therefore the sample is further separated into star-forming and quiescent galaxies. The morphological bi-modality of galaxy populations has been well studied in the local universe and is observed to persist to higher redshifts (e.g., . Several methods have already been developed to distinguish the two populations. In particular, the classification based on the U − V vs. V − J (UVJ) diagram has been used in many previous studies (e.g., Williams et al. 2009;Whitaker et al. 2011;. The separation of star-forming and quiescent galaxies based on the UVJ diagram has been shown to correlate with the separation based on UV+IR-determined specific star formation rates (hereafter sSFR ≡ SFR/M * , e.g. Williams et al. 2009) and based on SED-determined sSFRs (e.g. Williams et al. 2010). The inverse of the sSFR defines a timescale for the formation of the stellar population within the galaxy, which is why it is often used as a diagnostic of quiescence. Galaxies are labeled as star-forming if their rest-frame (U − V ) and (V − J) colors satisfy, (1) In the top panel of Figure 2, we plot the distribution of galaxies satisfying the limiting S/N and redshift criteria onto the UVJ diagram. The bi-modality of galaxy populations can be clearly seen in the UVJ diagram, where quiescent galaxies lie along the upper envelope (above the solid black line) and star-forming galaxies along the lower envelope. In total, the number of deconvolved star-forming galaxies that satisfy the color, redshift, and S/N criteria is 22,156. In the bottom panel, we show the distribution of these galaxies in the SFR-M * plane, where the SFRs are determined from the UV and IR fluxes (details of how these are derived are presented in . The total number of SFGs in this study is approximately 1-2 orders of magnitude larger in comparison to previous HST studies of the clumpy fraction of SFGs at similar redshift ranges (e.g., ∼700 in Wuyts et al. 2012 and∼3300 in Guo et al. 2015).

Finite Resolution Deconvolution
One of the main features of our deconvolution algorithm is that images are deconvolved to a finite resolution. Deconvolved images therefore have a narrower point spread function (PSF) of their own. Such finite resolution deconvolution ensures that solutions are well-constrained and well-sampled.
The following equations are generalized from Cantale et al. (2016a). We encourage the reader to review Cantale et al. (2016a) for the technical discussion on finite resolution deconvolution.
Here, we use an asterisk to denote the convolution operator. An image (D) taken from a ground-based optical telescope can be mathematically expressed as, where P SF is the point spread function and Z is the noise. M is a model of the image that is unaffected by atmospheric and instrumental blurring. The PSF of an image can further be expressed as, where g is the target PSF, taken to be a 2D Gaussian function, and P is the kernel that transforms g to PSF. Within FIREDEC, the construction of P consists of two fits: an analytic fit and then a numerical fit. The analytic fit models a Moffat-like profile to obtain the main characteristics of the PSF (i.e., the central core and the "wings" of PSF). Instead of fitting an exact Moffat profile, the analytic function is estimated as a sum of three elliptical Gaussians. The numerical fit is then needed to further model structures of the PSF that are not properly captured by the analytic fit. The numerical fit is obtained by deconvolving the residual of the analytic fit by the target PSF. When images are sampled onto a finite grid of pixels and affected by noise, deconvolution can become an ill-posed problem where many solutions are compatible with the data. However, a solution to deconvolution can be obtained by minimizing a cost function that include a regularization term, where the summation goes over each pixel element and σ represents the noise map for the observed image. Since FIREDEC deconvolves to a finite resolution, the cost function is written as, where B = M * g. For both equations, the first term is a chisquared term and the second term is a regularization term, H, weighted by a constant Lagrange parameter, λ. The purpose of the second term is to penalize high frequency components that arise during deconvolution. Within FIREDEC, the regularization term is defined to be, whereB is the denoised version of B, orB = φ(B). The denoising function implemented in FIREDEC, φ, is based on wavelet transform and is reviewed by Cantale et al. (2016a). Their implementation of the regularization term takes into accounts the dynamical range of astronomical images. As the S/N within an image can vary strongly, such implementation ensures that high S/N regions are treated in a similar way as low S/N regions during the deconvolution process, without having to define a spatially varying λ.
The determination of the Lagrange parameter is also crucial, as it ultimately dictates the importance of the regularization term within Equation 5. Traditionally, finding the optimal value for λ requires testing different values until a relative good deconvolution is achieved. This is checked by visually inspecting the noise within the deconvolved images. A quantitative analysis is used within this study to determine λ, and is outlined in Appendix A. Note that two different Lagrange parameters are determined in Appendix A. A Lagrange parameter is needed for the numerical fit of the PSF, where the fit is obtained by deconvolving the residuals of the analytic fit with the target PSF using a cost function similar to Equation 5. A separate Lagrange parameter is required for the deconvolution of the ground-based image (D) to obtain a model of the image with the target PSF of g (B).

Deconvolution in Practice
Instead of deconvolving the entire image of the COSMOS field in a given filter (which can become computationally intensive), we instead deconvolve smaller image cutouts that are centered on each galaxy. Doing so also ensures that that we account for any variations in the PSF across the full COS-MOS field. Indeed, observations of the COSMOS field consisted of multiple observing runs that extended over several years, leading to varying image qualities. Figure 3 shows the seeing variations for each filter. Each map covers an approximate area of 1.4 deg×1.2 deg, i.e., the full COSMOS field. The maps are binned into 17 × 15 different regions, where each region is 25 arcmin 2 . The colormap shows the difference between the median FWHM within a region and the median FWHM of the field (i.e., ∆FWHM). Note that in the optical Subaru filters, the full mosaic consists of a 3×3 grid of SuprimeCam pointings. The variation in image qualities in each of those pointings relative to each other is quite visible in some filters (e.g., IB505 and IA624). There is also PSF variations within the pointings (which is common in large mosaic cameras like SuprimeCam), demonstrating the importance of determining a local PSF when performing deconvolution.
Different deconvolution kernels are created based on the local PSF in order to deconvolve images to the same angular resolution. SExtractor is used to find all potential PSF stars within the field. The brightest stars that are flagged Figure 3. The seeing variations of COSMOS for eaphotometric filter. We bin each mosaic into pixel regions of 25 arcmin 2 , and use SExtractor to calculate the FWHM for all the stars. The colormap shows the differences between the median FWHM within each pixel and the median FWHM over the mosaic. as unblended and unsaturated are selected as suitable stars for PSF modeling. We also verify that each star does not have any close companions by inspecting the residual of its PSF fit for structures that are not associated with the star. It should be noted that the deconvolution kernel for each galaxy is modeled using the closest suitable PSF star to the galaxy. The same star may not be necessarily used as the PSF star in different filters as its magnitude can vary.
At 0.5 < z < 2, the effective radius of galaxies are observed to be ∼3−5 kpc (e.g., van der Wel et al. 2014), corresponding to an angular size of ∼1−2". We use an image cutout size of 7.8" by 7.8", which is roughly 4 effective radii bigger than the estimated size of galaxies. The target resolution is set as 0.3", which corresponds to a gain of a factor of 2−3 in resolution per band. A resolution of 0.3" also enables us to probe a physical size of ∼2.4 kpc at z ∼ 1. While this is approximately several times bigger than the estimated size of clumps (∼1 kpc; e.g., Soto et al. 2017), Section 6 shows that star-forming clumps can still be probed, even at this coarser resolution.
The target resolution of 0.3" is chosen to avoid oversampling pixels and mitigate noise that arise from deconvolution. Our ground-based images typically have a FWHM of ∼0.7" and a pixel scale of 0.15", corresponding to a sampling factor of around ∼4 for the PSF. Deconvolving to a higher resolution would require for smaller sampling interval by each pixel to allow for higher frequency components, which include both signal and noise. At our target resolution, each pixel is only resampled into 2 2 = 4 sub-pixels. The 4× resampling factor allows for the angular resolution to be improved to a level that enable us to partially resolve galaxies and detect clumpy structures, but also does not oversample the data to a point where it is difficult to distinguish clumps and noise. As we demonstrate in the remainder of the paper, our deconvolved images look like HST in similar filters and our measurements are also consistent with HST-based studies. This shows that this choice in resolution is optimal for maintaining sufficient S/N per sub-pixel, yet still recovering key features of galaxy properties that are required to classify them as clumpy or not.

Examples of Deconvolved Images
We compare the deconvolved images to HST images in order to see how well the deconvolution parameters are determined for the COSMOS data. Note that a direct comparison to validate whether the deconvolved images match the ground-based images (given the modeled PSF) was also presented in Cantale et al. (2016b) at lower z (≤ 0.8).
The COSMOS field was imaged in the F814W filter, with partial imaging in the F475W filter. There are also available F606W images from the CANDELS/COSMOS survey (Grogin et al. 2011;Koekemoer et al. 2011), and F160W images from COSMOS-DASH survey (Mowla et al. 2019). In total, there is an area of ∼0.03 deg 2 that has all 4 filters. We emphasize that HST images are only used to verify structures that are resolved by FIREDEC. These images are not used as additional constraints for deconvolution.

Quantitative Comparison with F814W
A quantitative comparison between the Subaru z + and F814W images is complicated by the fact that the two filters are not truly identical. The wavelength coverage of F814W extends from ∼7000−9500Å and its transmission curve peaks around 7500Å, while the Subaru broadband z + filter spans from ∼8000−10000Å. The PSF of both images also differ: z + -band images are affected by atmospheric blurring and their PSF can be characterized as a Moffat profile, while the PSF of F814W images is diffraction-limited, with large and extended PSF wings. Some differences between show the ground-based and HST images, while the deconvolved images are shown in column (c). We PSF-match the HST images to the deconvolved images at an angular resolution of 0.3" (column d). The differences between the deconvolved and PSF-matched HST images are shown as R(F814W, 0.3") in column (e). In general, the residual maps illustrate that deconvolved structures are consistent with HST images, with only minor structures within the central regions.
the deconvolved z + images and their HST counterparts can therefore be attributed to the effects of the morphological Kcorrection, where galaxy morphologies change as a function of wavelength, and differing PSFs. Keeping these points in mind, we present the quantitative comparison between the deconvolved and HST images.
The following notations are used to simplify the descriptions: P (x, a) represents PSF-matching an image x to a resolution of a, while R(x, a) represents taking the differences between an image x and the deconvolved images at that resolution. The F814W images are first degraded to match the resolution of the deconvolved images. This is done by convolving them by a Gaussian function so that their resolution becomes 0.3". The PSF-matched HST images are shown as P (F814W, 0.3") in column (d) of Figure 4, and the differ- Figure 5. Examples of deconvolution using multi-wavelength inages. Each row corresponds to a different galaxy. Panels (a) and (b) show examples of deconvolved images at a resolution of 0.3" in the rest-frame UV and optical, respectively, for galaxies between 0.5 < z < 1. Panel (c) shows the composite images. Within each panel, we show the ground-based, deconvolved and HST image for that galaxies from the left to right column. These galaxies are chosen mainly for their clumpy morphologies as observed in the deconvolved images. The HST/ACS F606W and F814W images are PSF-matched to the F160W image, with an angular resolution of 0.18". Similarly, the HST composite images are shown at a resolution of 0.18". ences between them and the deconvolved images are shown in column (e) as R(F814W, 0.3"). We also account for any misaligned astrometry by cross-correlating the HST images with the deconvolved images. In general, we find that the residuals of R(F814W, 0.3") are fairly uniform near the outskirts of galaxies, with minor structures near the central regions. These residuals are most likely due to the imperfect characterization of the PSF since it is impossible to know the exact profile of the PSF at the location of each galaxy due to the PSF variations across the field. However, the outer regions of the galaxies within the residual maps show negligible structures and indicates that key aspects such as clumps are clearly captured by deconvolution. In particular, examples such as ID116369 and ID118818 further illustrate that clumps are barely visible in the ground-based, but show nicely within the deconvolved images.

Multi-Wavelength Comparison
Figures 5 and 6 compares the deconvolved multiwavelength images of several galaxies with their corresponding HST images. The information across the different filters are not used as additional constraints for deconvolution, but instead the images are deconvolved separately. Each row corresponds to a different galaxy, while the columns show the rest-frame UV and optical imaging, along with the composite image of each galaxy. We show galaxies that have been imaged in multiple filters by HST. The galaxies are also selected mainly for their clumpy morphologies as observed in the deconvolved images. A visual inspection of the deconvolved and HST images reveals that deconvolved images (FWHM = 0.3") are strikingly similar to their HST counterparts (FWHM = 0.18"), even at the coarser resolution. Structures resolved in the deconvolved images are also observed within the HST images. Examples such as ID119634, ID124963 and ID115641 demonstrate the potential of FIREDEC where deconvolution recovers distinct clumps that are barely discernible in the ground-based V or z + -band images. The deconvolved structures are also observed across different filters, which further indicates that these structures are real, and not artifacts of deconvolution. A closer inspection of the deconvolved images from Figures 5 and 6 also shows that deconvolution is limited by the S/N of the ground-based images. In cases where the ground-based images have low signal, the low contrast between signal and noise causes deconvolved noise to appear more prominently as clumpy structures. This is noticeable particularly within the background of low S/N images such as the deconvolved V -band of ID117052 and ID117310. We stress that this does not have an effect on our analyses as noise is not correlated across different filters. Combining multi-wavelength observations to model the spectral energy distributions mitigates the effects of noise in our study of starforming clumps. Indeed, since clumps are identified using the rest-frame U and V surface brightness maps (see Section 5), noisy structures that are presented in single images are not detected within these maps. Finally, while individual clumps may not be resolved at the target resolution of 0.3", we iterate that the main goal of this paper is to broadly study the morphology of galaxies and to showcase the capability of finite resolution deconvolution. As Figures 5 and 6 show, deconvolution, while imperfect, offers a gain in angular resolution and presents a possible solution to obtain observational data that are similar in angular resolution to HST imaging.

METHODOLOGY FOR CREATING SPATIALLY-RESOLVED STELLAR POPULATION MAPS
The following sections give an overview of how multiwavelength images are combined so that resolved stellar properties can be inferred from SED models. We explain how star-forming clumps are detected within the resolved stellar mass density and surface brightness maps, and give a formal definition for what constitutes as a clumpy star-forming galaxy.

Segmentation Maps
Segmentation maps are needed to distinguish different objects within an image cutout and assign specific pixels to them. To create these maps, we use a watershed segmentation algorithm from scikit-image (van der Walt et al. 2014). The watershed algorithm treats the intensity of a grey-scale image as an overturned topological surface. The algorithm then uses defined markers to flood the valley until different markers merge to create the boundary of the segmentation maps. In practice, the markers are defined to be the coordinate of objects within the image cutout. For each image cutout, we mask the background, and the masked image is then normalized by the negative of its sum to create the overturned topological map. The background of image is determined by 3-σ clipping the image.

Misaligned Astrometry
In general, COSMOS has good astrometry. The error associated with the astrometric solutions is small (e.g. < 0.2" for the optical data and 0.1" for the NIR data; Capak et al. 2007;McCracken et al. 2012). While this is sufficient for photometry, it is not ideal for resolved studies of galaxies. Notably, such offsets are roughly the size of clumps, and therefore can have significant impact on our clump analyses. This is particularly true for compact galaxies, where a misalignment at the shorter rest-frame UV wavelengths can result in the misidentification of off-center clumps.
While the morphology of galaxies can change considerably between the optical and NIR, it is comparable between filters that are adjacent in wavelength. As a correction for potential offsets, we calculate the relative offsets in pixel coordinates by iteratively cross-correlating each photometric cutout to the adjacent photometric cutout at longer wavelength. This is applied on the ground-based images. To this end, we fix the astrometry for the K s -band images, and align the other image cutouts with respect to it. We use the integrated S/N of the galaxy (taken from the UltraVISTA catalog) as an indicator of whether a correction should be made. If the integrated S/N for a given photometric cutout is less than 5, we do not apply any corrections as there is already very little signal within the image cutout.

Pixel Binning
Modeling the spectral energy distributions on the pixel-bypixel basis can lead to unreliable measurements of the underlying stellar populations as the photometry of individual pixels can have low S/N. A solution is to bin pixels, so that each bin has some minimum S/N and ensures that the correct stellar properties can be obtained. However, galaxies spatially exhibit a wide range of S/N, so binning methods such as the Voronoi binning technique (Cappellari & Copin 2003) are used to homogenize the S/N across the galaxies, given a constraint on the minimum S/N. An example is shown in Wuyts et al. (2012), where they binned pixels based on the H-band in order to study the underlying stellar populations of galaxies at z ∼ 2.
However, binning based on the NIR alone can cause localized signals such as clumps to spread over a large area as the morphology of galaxies can vary between the restframe UV to NIR. Recently, Fetherolf et al. (2020) improved the Voronoi binning technique by incorporating the S/N information across multiple filters in the binning process. Such binning technique takes into account the morphological changes, but requires high S/N data as additional constraints from multiple filters generally lead to bigger bin size (and subsequently lower angular resolution) for low S/N data. While the global S/N of an image does not change as flux is recentered during deconvolution, regions of low surface brightness that are not initially observed in ground-based images are uncovered by deconvolution. Since our galaxies are selected to have an integrated S/N as low as 20, we opt to bin pixels based on the photometric information from two filters, rather than multiple filters in order to optimize the S/N per bin and the bin size, while accounting for the morphological changes across wavelengths.
We modify Cappellari & Copin (2003)'s Voronoi binning technique so that there are two binning channels: one based on the filter that probes the rest-frame UV of the galaxy (i.e., dependent on the redshift) and the other based on the observed K s -band. In the original Voronoi binning technique, a pixel is accreted to a bin if the addition of that pixel improves the S/N of the bin. The accretion of pixels stop once the bin reaches the desired S/N threshold. In our modified version, the binning process is similar, except a bin is decided if it reaches a minimum S/N in either one of the two bands. We define the minimum S/N to be 5. The pixel noise of each deconvolved filter is estimated by taking the standard deviation of the flux within empty apertures. These empty apertures are randomly placed across the deconvolved images, and have a diameter of 0.15", which correspond to the original pixel scale of the ground-based images. As alluded to earlier, binning also degrades the angular resolution of the data, so the inferred stellar population maps are rescaled using the deconvolved images in order to recover the deconvolved resolution of 0.3". This is discussed in Section 4.2.1.

Resolved Spectral Energy Distribution Modeling
We use EAZY (Brammer et al. 2008) to compute the restframe U and V luminosities for each spatial bin for all galaxies. EAZY calculates the colors by integrating the best-fit SED through the redshifted filter curves over the appropriate wavelength range. We use the response curve defined in Maíz Apellániz (2006) for the U and V filters. We fix the redshift of each spatial bin to the photometric redshift of its respective galaxy.
The stellar populations for each bin are determined by individually fitting Bruzual & Charlot (2003)'s stellar population synthesis models to the 14-band B-to-K s SED with FAST/C++ 1 . This is essentially the C++ version of FAST (Kriek et al. 2009). A Chabrier (2003)'s initial mass function and Calzetti et al. (2000)'s dust extinction law are adopted, and we assume a uniform solar metallicity and exponentially declining star formation history. We allow log(τ ) to vary between 7.0 and 10.0, in increments of 0.2, and log(t) to vary between 7.0 and 10.0, in increments of 0.1. We also restrict t to be less than the age of the universe at the observed redshift of each galaxy. The visual attenuation (A V ) of each fit is allowed to vary between 0 and 4.
Note that the SED models make no distinction between clumps and non-clumpy regions, and therefore assumes that clumps have the same properties as the surrounding region (e.g., the same star formation history; SFH hereafter). While this may not be necessarily true, little is known about the SFHs of high-z clumps. We choose to use the τ -model as it is commonly assumed for distant galaxies and it has been shown to recover stellar masses and SFRs with typical systematic offsets of ∼0.3 dex (Lee et al. 2018 Modeling the SEDs of each spatial bin allows us to construct the stellar mass density and surface brightness in the rest-frame U and V. However, Voronoi binning effectively degrades the angular resolution of these maps. In order to recover the resolution of 0.3", light variations within the deconvolved images are used to spatially calibrate the stellar masses and surface brightness. In essence, the deconvolved images in the K s band are used to scale the stellar mass density maps as this is our best filter for probing stellar masses at our redshift range. Similarly, the surface brightness (in the rest-frame U and V ) are scaled with the respective deconvolved images that probe the rest-frame UV and optical.
To this end, the deconvolved images are first scaled so that the minimum and maximum values within each deconvolved 1 https://github.com/cschreib/fastpp/releases/tag/v1.3.1 image are between 0 to 1. An image, x, is therefore scaled using the following equation, Using the minimum and maximum values of the whole image (as opposed to a Voronoi bin) ensures that the relative brightness within the galaxy remains preserved. Given a measured quantity (Q) of a Voronoi bin, the resolution of 0.3" is recovered by scaling Q using the following relation, In the notation above, i represents pixel elements within the given Voronoi bin. Qx i is normalized by the sum of x i within the bin to ensure that the inferred stellar mass or luminosities from SED fitting is preserved after rescaling.

IDENTIFYING AND CLASSIFYING CLUMPY STRUCTURES
In the following sections, we address the spatial variations within the stellar mass density and surface brightness maps. This is done by creating normalized radial profiles (Wuyts et al. 2012). A normalized profile is a two-dimensional profile that simultaneously quantifies the colors and radial dependencies on stellar properties, and enables us to detect any spatial variations within the stellar population maps. We briefly discuss how a 2-dimensional profile is made, and the classification scheme used to identify clumpy regions. We refer the reader to Wuyts et al. (2012) for further details.

Normalized Light and Stellar Mass Profiles
The stellar mass-weighted center of the galaxy is adopted as a reference point for all measurements. This is motivated by the fact that the mass-weighted center should correspond to a location within a galaxy where the gravitational potential is strongest. The normalized brightness profiles are then constructed as follows.
We first associate the morphology of each galaxy with two elliptical parameters: the position angle of the galaxy and its axial ratio. The former is obtained by measuring the position angle of each pixel element within the segmentation map (as obtained from Section 4.1.1), and constructing a distribution of the position angles. The mode of that distribution is taken to be the position angle of the galaxy. Similarly, a distribution of the pixel elements' distance to the mass-weighted center is made, and the axial ratio of the galaxy is defined to be the ratio between the mode of the distribution and the extent of the distribution. For the rest-frame U and V surface brightness map, we then define the half-light radius (R e ) of a galaxy to be the semi-major axis length at which an elliptical aperture contains 50% of the respective light. This radius is derived by constructing a curve of growth using elliptical apertures (defined by the position angle and axial ratio), which are placed at the mass-weighted center. Having established a half-light radius, we measure the average surface brightness within the half-light radius (Σ e ). The surface brightness of each pixel element is then normalized by the average surface brightness, and its galactocentric distance is normalized by the half-light radius. In essence, the whole process normalizes the surface brightness of each pixel element and remaps them into a new parameter space. The constructed normalized light profile will typically decrease as a function of radius, with the exceptions being for when pixel elements have an enhanced surface brightness. These enhanced surface brightness are most likely due to star-forming clumps in high-z galaxies (e.g., Wuyts et al. 2012).
The same procedure is used to create the normalized stellar mass profile for each galaxy. The exception is that we compute the half-mass radius as the semi-major axis length at which 50% of the total stellar mass of the galaxy is contained. Note that the elliptical parameters are derived based on the segmentation map (as opposed to using the surface brightness) of the galaxy, so that the shape of the elliptical aperture is driven by the outer isophotes of a galaxy, rather than by individual bright clumps. Figure 7 shows four examples of the normalized radial profiles alongside the resolved stellar mass distribution (panel ii), U rest (panel iii) and V rest (panel iv) surface brightness maps. The composite image of the galaxy is also shown in panel (i), while the rest-frame (U − V ) color map is shown in panel (v). Since clumps are regions of enhanced surface brightness compared to the underlying disk, the radial profiles can be used to identify clumpy features within the galaxies. Panels (vi), (vii) and (viii) indicate whether pixel elements are identified as clumps, or part of the inner/outer "smooth" regions based on the stellar mass distribution and surface brightness maps. The classification scheme for whether pixel elements are clumps is later discussed in details in Section 5.2.

Examples of Normalized Mass/Light Profiles
A number of observations can be noted from the normalized profiles. All three radial profiles (panels vi, vii, viii) are flat within the half-mass/light radius. Since galaxies are typically described as Sersic profiles, the flattening of the radial profiles within the central regions suggests that our profiles are limited by angular resolution (and indeed they are, as we only deconvolve to a resolution 0.3", which is 2.4 kpc at z ∼ 1), and implies that it may be difficult to identify centrally-located clumps. However, the concerns for the non-detection of clumps near the central regions are not severe, as previous studies of clumpy galaxies based on HST images demonstrate that clumps are typically found near the outskirts of galaxies (e.g., Zanella et al. 2019). Even at higher resolution, it may still be difficult to identify clumps within the half-mass/light radius as the central brightness of galaxies are much higher. Similarly, Wuyts et al. (2012) also observed the flattening of the radial profiles with HST images, presumably from the finite resolution of 0.18" in the F160W filter.
At larger radii within the normalized light profiles, we find that the surface brightness drops and shows a considerable scatter of ∼1 dex at R R e . The large scatter in the U rest (and in some cases, V rest ) profile is driven by localized regions of enhanced brightness (e.g., ID114072 and ID119634). For the remainder of the paper, we refer to any off-center regions with an enhanced level of surface brightness as "clumps". This nomenclature is inspired by the morphological appearance of the galaxies as shown in Figure 7. However, we stress that this definition encompasses all regions that have an elevated surface brightness, regardless of their geometric shapes (e.g., structures such as spiral arms). It is not intuitive whether a distinction can be made between actual star-forming clumps and other structures based on the excess of surface brightness alone. The normalized profiles also indicate that there is a color gradient; we typically find redder colors within the inner region of the objects (R < 0.32R e ) as shown by the colormap of the radial profiles. The comparison between the mass and light profiles also shows that the stellar mass distribution of galaxies are generally smooth, even though their light profiles exhibit clumpy morphologies. Clumps are also easily detected within the U rest map compared to V rest . This is expected if we consider the M/L effect where younger stellar populations are brighter at the same stellar mass. Such effect is also observed through the fraction of clumpy galaxies (see Section 6.1).

Classification Scheme for Clumpy Structures
The methodologies for detecting clumps vary from studies to studies, and no universal way have been defined. A number of methods have been used to detect clumpy structures, but the majority are based on identifying pixels with enhanced intensity through a computer algorithm (e.g., Wuyts et al. 2012;Guo et al. 2015;Zanella et al. 2019;Huertas-Company et al. 2020). Our normalized profiles naturally allow us to perform a similar analysis. Following the classification scheme of Wuyts et al. (2012), we split the normalized parameter space into three regions; [inner] x < −0.5 where x ≡ log(R/R e ) and y ≡ log(Σ/Σ e ). Pixel elements are considered to lie in the inner, outer, or clump regime based on their location in the (Σ/Σ e , R/R e ) parameter space. Note that the dashed curve of Equation 9 does not truly define the light curve for all galaxies. This is why there are gaps between the normalized profiles and the dashed curve in Figure 7, and similar observations is also observed in the normalized profiles shown by Wuyts et al. (2012). These dividing lines were derived by Wuyts et al. (2012) based on their stacked profiles, and are fixed for all galaxies. We argue that the same lines are also sufficient for classifying clumpy structures in our galaxies as our mass/light profiles are obtained in a similar way. Each mass/light profile is obtained by normalizing the mass density/surface brightness of each Voronoi bin by the mean mass density/surface brightness of the galaxy. While the dividing lines do not vary from galaxies to galaxies, the normalization ensures that all clumpy pixels from different galaxies are mapped to a similar region within the normalized parameter space. As shown by Wuyts et al. (2012) and this study, it is visible that these divisions are sufficient for identifying. This is illustrated in the right column of the panels in Figure 7, where pixels are colorcoded based on their pixel type (i.e., inner, outer and clump as red, grey and yellow, respectively). In general, pixels that are identified as clumpy coincide with regions of blue color as observed in the composite images. The normalized radial profiles also illustrate this, as clumpy pixels generally have blue rest-frame (U − V ) color. Furthermore, for cases where the galaxies' light profile (e.g., bottom panels ofFigure 7) are smooth no clumps are detected.
However, the distinction between a clumpy and nonclumpy galaxy is still subjective: are all galaxies with clumpy pixels "clumpy"? In general, studies of star-forming clumps typically use the ratio of the clumps' luminosity to the total luminosity of the host galaxy (i.e., the fractional luminosity) to make the distinction. Wuyts et al. (2012) defined clumpy galaxies to be galaxies with clumpy pixels that contribute more than 5% of the total luminosity. On the other hand, Guo et al. (2015) and Shibuya et al. (2016) defined clumpy galaxies as those that have off-center clumps with a fractional luminosity of at least 8%. This threshold was derived by comparing the fractional luminosity of high-z starforming clumps and local star-forming regions that were artificially redshifted. The fractional luminosity of 8% represents a point where the number counts of high-z star-forming clumps deviates from local star-forming regions (see Guo et al. 2015 for more details). For consistency of comparison, we use the same definition of 8% as the dividing line between clumpy and non-clumpy galaxies. Given that deconvolution is a new approach to studying clumps, having such consistency while measuring the clumpy fractions from the HST and deconvolved data allows us to show that de-convolution can be a powerful tool for analyzing small-scale features within galaxies from ground-based images.
In essence, the resolved stellar mass and surface brightness maps are used to identify clumpy pixels. A galaxy is defined as clumpy in a given band if pixels originating from the clump regime contribute to at least 8% of the galaxy's total luminosity in the respective band. Similarly, galaxies are identified as clumpy in stellar mass if at least 8% of the total stellar mass of the galaxy originates from pixels in the clumpy regime of the normalized mass profile.
Again, it is emphasized that the nomenclature for "clump" includes all off-center structures that have an excess of surface brightness or stellar mass, regardless of their geometric shape. In particular, galaxies that appear to be edge-on (or galaxies that are highly eccentric with e > 0.8) can have pixels that lie further along the R/R e axis in the clumpy region within the normalized profile plots. We remove these galaxies from the analysis to avoid biasing our sample of clumpy galaxies. The eccentricity of the galaxies are calculated using the segmentation maps from Section 4.1.1, where we obtain the semi-major and semi-minor axes by fitting an ellipse to the edge of the segmentation maps. The total number of galaxies that are not eccentric, and therefore is included in the analysis of the paper is 20,185. For the rest of the paper, we refer to non-clumpy galaxies as "regular" galaxies.

THE CLUMPY FRACTIONS
In the following sections, we measure the clumpy fractions based on the surface brightness and stellar mass density of the galaxies (f U rest clumpy , f V rest clumpy , f mass clumpy ), and show the evolution of the clumpy fractions with redshift. Our measurements are also compared to other studies that are based on HST data. Finally, we also investigate the relation between the abundance of clumps and physical properties of the host galaxy such as stellar masses and star formation rates.

Evolution of the Clumpy Fractions
The evolution of the fraction of clumpy galaxies is shown in Figure 8. Our galaxy sample is binned into 3 different mass bins: 9.8 < log(M * /M ) < 10.2, 10.2 < log(M * /M ) < 10.6 and 10.6 < log(M * /M ). Within each mass bin, we also separate galaxies into redshift bins. The fraction of clumpy galaxies within each bin is calculated. We find that the clumpy fraction changes considerably depending on how clumps are detected, with f clumpy dropping by ∼20% when clumps are selected based on the V rest surface brightness maps compared to U rest . When clumps are identified based on the stellar mass distribution, we find that the general trend of f mass clumpy is similar to f V rest clumpy , but at lower amplitudes.
The clumpy fractions are found to increase with redshifts, irrespective of the stellar mass considered. However, mas- . The clumpy fractions for SFGs at 0.5 < z < 2. From left to right, clumpy galaxies are identified based on the rest-frame U luminosity, V luminosity, and the stellar mass distribution map. The green, red and blue markers are the clumpy fraction of the low, intermediate, and highmass bin respectively. The shaded region represents that errors derived from the Poisson error of the galaxy counts in those bins. The best-fit slopes for the clumpy fraction evolution are given in Table 2 and are plotted as the dashed lines.
sive galaxies are less likely to host clumpy structures compared to low-mass systems. The slope for the evolution of the clumpy fractions also varies, where f clumpy can decrease between 10 − 30% from the highest to lowest redshift bin. The evolution of the clumpy fractions can be quantified by fitting f clumpy to power-laws of the form C(1 + z) α . The values for each fit are listed in Table 2. We find that the slope of f U rest clumpy is shallower compared to f V rest clumpy and f mass clumpy . The clumpy fractions for the intermediate and high-mass bins flatten out at z < 1. However, the opposite is observed in the low-mass bin, where the clumpy fractions flatten out at z > 1, with the exception only for clumps detected in the stellar mass maps.
The flattening of the clumpy fractions as observed in U rest and V rest for low-mass systems at z > 1 could be a S/N effect. Since galaxies are defined as clumpy based on the fractional luminosity (or mass) contribution of clumps, low S/N clumps are generally not as easily detected against the low surface brightness of low-mass galaxies. We note that while our observed clumpy fractions could be biased toward lower values for such galaxies, previous studies that used a similar definition are also affected by a similar bias. To conclude, we list the notable points of our findings; • Low-mass galaxies are generally clumpier than highmass galaxies at all redshift, irrespective of the selection method for clumpiness (i.e., in U rest , V rest , or mass) • Regardless of stellar mass, the fraction of clumpy galaxies increases toward higher redshifts. Given that S/N decreases with redshift, this trend may be even stronger than what we measure.
• Galaxies have higher clumpy fractions when selected in U rest compared to V rest or stellar mass. While the  Table 2. Listed are values obtained from fitting the clumpy fractions to C(1+z) α . We find that the clumpy fraction of low-mass galaxies evolves slowly with redshift compared to higher mass systems. For a given mass bin, we find that the clumpy fraction evolves strongly with redshift when using the surface mass density map, followed by the Vrest and Urest surface brightness maps.
rest-frame UV trace younger stars, V rest probe older stellar populations and can trace stellar mass, which can explain why f V rest clumpy and f mass clumpy are similar.

Comparison to HST Studies on Clumpy Fractions
In Section 3.3, we tested our deconvolution algorithm by comparing deconvolved images to HST images across a wide range of wavelengths and showed that deconvolved multiwavelength images are generally consistent with their HST counterparts. As an indirect test of how well FIREDEC can resolve structures, we also compare our measured clumpy fractions to other measurements based on HST imaging (i.e. Wuyts et al. 2012, Murata et al. 2014, Guo et al. 2015  . Comparison of f U rest clumpy with other HST-based studies. The top panel shows the clumpy fractions of the whole galaxy sample, without binning galaxies by stellar masses. The middle and bottom panels shows the clumpy fractions for the low and highmass bins, respectively. Our clumpy fractions are derived from the Urest surface brightness maps. The error bars are based on Poisson counting statistic. The black and white markers are values based on HST data, where clumps are detected within single-filter images that probe the rest-frame UV. In general, we find our values are broadly consistent with these studies, and are only slightly lower at higher redshift. Shibuya et al. 2016). While our galaxy selection differs slightly from previous HST studies (e.g., in term of stellar masses and redshift range), such a comparison still serves as a fiducial test for whether deconvolution can correctly recover spatial information within ground-based images and produce results with similar accuracy as HST-based studies.
The top panel of Figure 9 shows the clumpy fractions for all galaxies from different studies. The clumpy fractions are derived using the rest-frame U -band, either by modeling the SEDs (as in the case for this study and Wuyts et al. 2012's), or by using a selection of filters to probe the rest-frame UV of galaxies at different redshifts (i.e., Guo et al. 2015;Shibuya et al. 2016). The only exception is Murata et al. (2014), who used the F814W filter to identify clumps in all galaxies at 0.5 < z < 1. In general, the overall distribution of our f clumpy is in agreement with other measured values.
We find that f clumpy drops from ∼50% to ∼30% between 0.5 < z < 2. Wuyts et al. (2012) found a slightly evolving clumpy fraction that decreases from 60% to 57% between the redshift bins of 0.5 < z < 1.5 and 1.5 < z < 2.5. Similarly, both Guo et al. (2015) and Shibuya et al. (2016) reported values that are slightly higher compared to ours, but show a similar trend in the redshift evolution. Our lower fractions could be attributed to our coarser resolution of 0.3", which can blend and smear clumps. We also note that there are inconsistencies between our clumpy fractions and those reported from Murata et al. (2014), where a steeper slope was observed in their study. While they found that f clumpy drop from 35% to 5% between 0.2 < z < 1, their measurements are based on a single HST filter. The discrepancy could therefore be affected by the fact that galaxy morphologies change as a function of wavelength. Such effects have been observed in this study and Wuyts et al. (2012)'s, where the clumpy fractions are lower when clumps are selected at longer rest wavelengths. If we take into account of the morphological K-correction, the slope of the clumpy fractions from Murata et al. (2014) should dampen and be better aligned to ours. Considering the variety of methodologies and datasets used to distinguish star-forming clumps, the systematic differences between all studies are small.
The middle and bottom panels of Figure 9 show the clumpy fraction at two different mass bins. As an approximate comparison of the literature, we split our sample into a low-mass bin at 9.8 < log(M * /M ) < 10.4, and a high-mass bin at 10.4 < log(M * /M ). In general, we find that our measurements are similar to the literature, with minor discrepancies. The evolution of f clumpy is slightly different, with our values showing a slower evolution with redshift compared to both Guo et al. (2015) and Shibuya et al. (2016). In the low-mass bin, our fraction of clumpy galaxies drops by ∼15% between 0.5 < z < 2, compared to the ∼20 − 30% decrease observed in the HST-based studies. At z ∼ 2, we also find that our values are lower (by ∼ 15%) compared to those reported in Shibuya et al. (2016). For the high-mass bin, our clumpy  Figure 10. Galaxies are binned within the M * and SFR parameter space, where only bins with more than 10 galaxies are shown. The colormap shows the clumpy fraction of each bin. We further separate the sample into two redshift bins, with a low-z bin (z < 1) on the left and the high-z bin (z ≥ 1) on the right. The red line is the star-forming main sequence from Whitaker et al. (2012) using the fiducial redshift of 0.75 (1.5) for the low-z (high-z) bin.
fraction decreases by ∼30% from z ∼ 2 to z ∼ 1, whereas studies such as Guo et al. (2015); Shibuya et al. (2016) reported a similar decrease, but on a shorter timescale. However, in general, our fractions are found to be mostly within the errors of other measurements. Obtaining resolved measurements of the physical properties of clumps and their host galaxies are crucial if we truly want to understand the role that clumps play in galaxy evolution. We remind the reader that the majority of the HSTbased studies presented within this paper are limited to one filter per galaxy, and are therefore not adequate for SEDs modeling. The only exemption is the resolved studies of ∼700 galaxies at 0.5 < z < 2.5 by Wuyts et al. (2012). Within this comparison of f clumpy , we show that image deconvolution can be used to further extract information from existing ground-based images. Indeed, in the binary view of clumpiness, we find that the deconvolved data are largely consistent with what was observed by previous HST-based studies. However, in addition to tracking the fraction of clumpy galaxies, the multi-wavelength deconvolved data allow us to spatially model SEDs and obtain resolved stellar mass density and surface brightness of ∼20,000 galaxies. Our larger sample also resulted in smaller statistical uncertainties in the clumpy fractions compared to the HST-based studies.

Dependency of Clumps on sSFRs
Since SFRs and stellar masses are fundamental properties of galaxies, it is intriguing to question whether clumpy morphologies are associated with galaxies whose global properties deviate from those on the star-forming main sequence. In the following analyses, we further exclude galaxies with a specific star formation rate (sSFR ≡ SFR/M * ) less than 10 −10 yr −1 in order to compare our UVJ-selected SFGs to those on the main sequence. However, our result does not change even when we include these galaxies. Figure 10 shows the clumpy fraction (based on U rest ) along the M * -SFR parameter space. We bin galaxies based on their mass and SFR, and measure the clumpy fraction of each bin. Our sample is also separated into two different redshift bins; low-z (z < 1) and high-z (z ≥ 1). The red dashed line denotes the star-forming main sequence relation from Whitaker et al. (2012), using the fiducial redshift of 0.75 and 1.5 for the low and high-z bin respectively. Our SFRs are taken from the COSMOS/UltraVISTA catalog, which is derived from the UV+IR luminosities.  showed that the star-forming main sequence derived from the data are consistent with other measurements such as Whitaker et al. (2012).
Galaxies have higher SFRs toward higher redshifts, which is observed as the upward shift of the star-forming main sequence relation. Interestingly, the population of galaxies that lies above the main sequence generally have a higher clumpy fraction compared to those that are situated below it. The trend persists in both redshift bins and appears to be independent of stellar mass. This is not inconsistent with Figure 8 if enhanced SFRs are correlated with higher clumpy fractions, such that the increase in the clumpy fractions at high-z is explained as an increase in SFR within galaxies. At fixed stellar mass bins, clumpy galaxies will generally have higher SFRs compared to non-clumpy galaxies. Equivalently, clumpy galaxies have higher sSFRs compared to the average galaxies at fixed stellar mass bin. Although biases could be introduced as a result of our fine SFR and stellar mass bins in Figures 10, our findings are in good agreement with the results from Murata et al. (2014). They found that the distribution of clumpy galaxies within the M * -SFR parameter space typically reside in the upper envelope of the main sequence.  Figure 11 shows the dependence of the clumpy fraction on the sSFR, binned by stellar masses. We again separate the sample into two different redshift bins, and compare our results to similar works that are based on HST imaging. For both redshift bins, we find that f clumpy increases by ∼20% from our lowest sSFR bin to the highest bin. Galaxies in the high-mass bin (log(M * /M ) > 10.6) typically have clumpy fractions that are slightly lower compared to those with lower stellar masses, although both measurements are within errors of one another. While studies such as Murata et al. (2014) and Shibuya et al. (2016) used slightly different stellar mass bins, they reported similar trends. At z < 1, Murata et al. (2014) reported slightly lower values for the clumpy fraction as it increases from 0% to ∼20% with sSFRs. At z ≥ 1, Shibuya et al. (2016) reported f clumpy that are higher compared to ours, and their strong f clumpy evolution with redshift was also observed in the bottom panels of Figure 9. While there are differences in the normalization factor of f clumpy , they also found a similar trend where f clumpy is increasing with sSFR in the z ≥ 1 bin, but only at the low-mass end of galaxies (log(M * /M ) < 10.6). These differences could Galaxies are separated into two redshift bins, with the blue and red color denoting the low-z bin (z < 1) and the high-z bin (z ≥ 1), respectively. We find an increase in the clumpy fraction toward higher ∆ log(sSFR). At fixed ∆ log(sSFR) bin, the measured clumpy fraction is larger at higher redshift (z ≥ 1). The error bar represents the Poisson count statistic. The white markers are measurements made by Huertas-Company et al. (2020). They classified clumpy galaxies as those having clumps with stellar mass greater than 10 10 M , and measured ∆ log(sSFR) in relation to the M * -sSFR relation from Fang et al. (2018) be systematic errors as various methodologies are used for detecting clumps. For example, Shibuya et al. (2016) relied on the F606W filter to detect clumps for all galaxies at z < 2, as opposed to this work where a selection of filters is used. Their measurements are also associated with higher uncertainties, with error bars spanning between 20−40% in f clumpy , compared to our errors, which benefited from having a larger sample size. On the other hand, it is possible that lower clumpy fractions are observed in deconvolved images for clumps at z > 1.5 due to having lower S/N. In order to investigate how the relation between f clumpy and sSFR changes compared to galaxies on the main sequence, we plot the clumpy fraction as a function of ∆ log(sSFR) in Figure 12. Here, ∆ log(sSFR) measures how far the specific star formation rate of our galaxies deviate from the M * -SFR relation of Whitaker et al. (2012) at fixed stellar mass. Our sample is again binned into two redshift bins, and ∆ log(sSFR) is binned between −1 to 1. We find that the clumpy fraction is slightly dependent on ∆ log(sSFR) for both redshift bins. Lower clumpy fractions are found in galaxies that have low sSFRs compared to the main sequence, with the clumpy fractions increasing by ∼10% at the highest ∆ log(sSFR) bin. This is different what was reported by Huertas-Company et al. (2020), where they observed lower clumpy fractions that show no dependencies on ∆ log(sSFR).
One would expect that the majority of star-forming clumps to be found in galaxies above the main sequence if in-situ clump formation is mainly dependent on gas mass density. Since the sSFR is an indicator of the on-going star formation and could be correlated to the gas density, our observations provide a self-consistent picture that links high SFR and high gas fraction to clumpy morphologies. Indeed, the observed increase of the clumpy fractions toward high redshift in Figure 8 can be attributed to the increase in star formation rate within galaxies and higher gas mass fractions. We also note that our definition of clumpiness is binary. While the threshold for determining whether a galaxy is clumpy or not is physically motivated, such binary classification can affect the relation between clumps and sSFRs. In Appendix B, we show that similar results to Figures 10 and 12 are obtained when using the fractional U rest contribution of clumps to their host galaxy (hereafter, C UV , as opposed to the clumpy fraction).

EVOLUTION OF CLUMPY SFGS LINKED AS PROGENITORS-DESCENDANTS USING ABUNDANCE MATCHING
The stellar mass assembled in galaxies evolves with cosmic time as stellar masses can grow from mergers and ongoing star formation. Inferring the mass assembly history of galaxies can be challenging, as it requires us to accurately link progenitors to descendants. This process is complicated as galaxies can only be observed at one snapshot in time. It is possible to connect descendant galaxies to their most likely progenitors by observing how galaxy populations change in different parameter spaces (e.g., number density). The simplest method to derive the masses of progenitors is by matching their cumulative number density. This method begins with the assumption that as galaxies grow in masses, their cumulative number density would remain constant if galaxies retain their rank order in stellar mass. However, the assumption that the cumulative number density is fixed through redshift as galaxies evolve is only true if there are no mergers and galaxies have similar SFHs. The effect of mergers can be predicted using abundance matching in simulations (e.g., Behroozi et al. 2013;Wellons & Torrey 2017). Hill et al. (2017) showed the progenitor mass evolution for four galaxy populations with masses of log(M * /M ) ∼ 11.5, 11.0, 10.5, and 10.0 at z = 0.1. This was done by connecting the observed mass functions to the number densities from Behroozi et al. (2013). We use the mass evolution calculated by those authors to trace the progenitors of two populations with with log(M * /M ) ∼ 10.9 and log(M * /M ) ∼ 10.2 at z ∼ 0.7. At each progenitor mass bin, we measure the fraction of clumpy galaxies. The evolution of the clumpy fractions for these progenitors is shown in Figure 13. We find that the clumpy fractions decrease as galaxies grow in stellar mass. This dependency suggests an interplay between clump formation and global properties such as gas availability and stellar masses. Cosmological simulations suggest that gas accretion onto massive systems is mainly through hot-mode accretion and low-mass systems through cold-mode gas accretion (Kereš et al. 2005), which implies that cold gas accretion is dominant at higher redshifts when galaxies are less massive. As galaxies grow in mass, there is a suppression of gas supply, and star formation activity decreases, leading massive galaxies to passively evolve into the "red and dead" populations. However, note that cold gas accretion may also be present in massive systems at z ∼ 2 (Dekel & Birnboim 2006), which explains how some massive galaxies at highz can exhibit clumpy morphologies. Our results therefore support in-situ clump formation, where the unstable, gas-rich disk of high-z galaxies is supplied by cold gas accretion.

Evolution of the Clumpy Fractions using Abundance Matching
Furthermore, when comparing the evolution of the clumpy fraction for the two galaxy populations, we find that the progenitors of the low-mass population generally have higher clumpy fractions compared to the high-mass population at fixed redshift. This result illustrates the so-called "downsizing effect" in the clumpy morphology, which have been noted by several studies (e.g. Elmegreen et al. 2009b, Murata et al. 2014): while massive clumpy galaxies exist at z > 1, such clumpy morphologies are mainly observed in low-mass systems in the local universe.

Contribution of Clumps to the Mass of SFGs
Star-forming clumps are often suggested to be an important part of the stellar mass growth of galaxies due to their ubiquity in star-forming galaxies during cosmic noon. However, measuring the contribution of clumps toward the stellar mass of the galaxy is difficult as it requires accurate measurements of the photometry of the clumps. Even with current state-of-the-art observatories, resolution effects can cause an overestimation in the size and mass of clumps (e.g., Tamburello et al. 2017;Cava et al. 2018;Meng & Gnedin 2020). Noting this caveat, the following analysis looks at the fractional contribution of the stellar mass that lies within clumpy structures to the integrated mass as galaxy populations grow over time.
Following the progenitors of the two populations in Figure 13, we list the fractional contribution of the integrated clump mass (as identified in U rest ) to the overall stellar mass budget of SFGs within the respective redshift bin in Table 3. For galaxies with log(M * /M ) ∼ 10.9 at z ∼ 0.7, we find that their clumps have a higher fractional mass contribution at earlier times: clumps account for ∼5% of the total mass at z ∼ 1.7, while only ∼1% at z ∼ 0.7. At face value, this . The evolution of clumpy fraction for the progenitors of galaxies with masses of 10 10.2 M (squares points) and 10 10.9 M (circle points) at z ∼ 0.7. The mass evolution for progenitors is obtained from Hill et al. (2017), by connecting the observed mass functions to the number density from Behroozi et al. (2013). At low redshifts, clumpy morphologies are found mainly in low-mass systems, while a similar clumpy fraction is only observed in more massive galaxies at higher redshift.
result suggests that the stellar mass growth within galaxies relies in part on clumps. However, the mass of the progenitors are also evolving. The mean galaxy mass decreases by 0.6 dex from z ∼ 0.7 to z ∼ 1.7. This corresponds to a decrease by a factor of ∼4 in the mean galaxy mass, which is only slightly lower compared to the increase of a factor of ∼5 in the fractional mass contribution from clumps. Therefore, the total mass within clumps only increases slightly toward higher redshifts, if at all. The clumpy fraction also decreases with cosmic time (i.e., from ∼50% to ∼25%), suggesting that galaxies will host fewer, but more massive clumps as they evolve in mass. A similar result is obtained when we compare two populations of different masses at the lowest redshift bin (z ∼ 0.7): the more massive population have a lower clumpy fraction, but clumps from both populations have similar fractional mass contribution. When comparing two populations of similar mass, but at different redshift bins (e.g. log(M * /M ) ∼ 10.2 at z ∼ 0.7 and z ∼ 1.7), we find that the population at high-z have a higher fractional mass contribution from clumps. The clumpy fraction is also higher, supporting the scenario that clumps become more important to stellar mass growth at earlier cosmic times. We again remind the reader that the nomenclature for clumps encompasses all structures that have an enhanced surface brightness (or mass) with respect to the underlying surface brightness. The measured clump mass in this Section is based on the identification of clumps using the normalized light profile, and effects such as clumps blending and background flux are not taken into account. These effects may cause an overestimation in the measured mass values. In future studies, it will be beneficial to isolate individual clumps and measure their properties (e.g., Guo et al. 2018 Table 4. Summary of clumps in the progenitors of galaxies with log(M * /M ) ∼ 10.9 at z ∼ 0.7. Similar to Table 3, we list the clumpy fraction and the fractional contribution for the progenitors.
At fixed stellar masses, clumpy SFGs are more likely to have higher SFRs compared to regular SFGs. The distribution of clumpy and regular SFGs within the M * -SFR parameter space (Fig. 10) is not surprising and can be explained within the context of disk instability. The stability of a gas disk is described by the Toomre's Q parameter (Toomre 1964), where σ g and v c represent the gas velocity dispersion and circular velocity, and Σ g is the gas density at radius r. The gas disk is considered to be stable if Q > 1. It is clear that in the simple case of high cold gas density, the gas disk is more susceptible to gravitational instabilities, and can fragment into clumpy star-forming regions. The relation between gas density and the star formation rate density was established by Kennicutt (1998), and is observed to persist for high-z galaxies ). If we assume that this relation holds for our galaxies, then at fixed stellar masses, galaxies with high sSFRs (i.e., galaxies above the star-formation main sequence) are expected to have high gas density and subsequently low Toomre parameter values. This could explain the relation between ∆ log(sSFR) and the clumpy fractions in Figure 12. Indeed, the formation of massive clumps within a gravitationally unstable gas-rich disk have been often noted in the literature (e.g., Escala & Larson 2008;Dekel et al. 2009). It has also been observed in numerical simulations, where massive clumps are more common in galaxies with a high gas fraction compared to those with a lower gas fraction (Fensch & Bournaud 2020). Similar relations between the morphologies, star formation rate and gas fraction are also observed within the local universe. Locally, galaxy morphologies are based on the Hubble Sequence, with late-type galaxies referring to spiral and irregular galaxies, and early-type galaxies referring to ellipticals. Eales et al. (2017) found that morphologies change gradually with respect to sSFRs and stellar masses, where late-type galaxies typically have higher sSFRs at lower masses and early-type galaxies are more massive and have lower sSFRs. As late-type galaxies tend to have a higher cold gas reservoir compared to early-type galaxies (Calette et al. 2018), the difference in the morphology of nearby galaxies can also be explained through gas availability and star formation.
Finally, the Toomre parameter suggests a crucial relation between the gas kinematics and the gas fraction. If the disk is expected to be marginally stable, we can set Q g = 1 and rearrange Equation 10 to express v c /σ g in term of the gas fraction (see Genzel et al. 2011 andGlazebrook 2013), At large scales, the kinematics of high-z galaxies are found to be mostly rotating disks, but with rather large gas velocity dispersion (σ ≈ 30−90 km s −1 ; see Förster Schreiber et al. 2009). Theoretically, lower v c /σ g values and larger gas mass fractions in high-z galaxies naturally lead to the larger starforming complexes in comparison to the less gas-rich and less turbulent disks of low-z galaxies (e.g., Escala & Larson 2008). What drives large velocity dispersion in high-z galaxies and why such high dispersion are not observed in the local universe are both interesting questions. Galactic disks are established to be in global equilibrium, and the marginally stable state (Q ∼ 1) is suggested to be due to a self-regulating feedback loop. In the state where Q < 1, mechanisms that drive turbulent motion are required to bring the disk back to Q ∼ 1.
The picture of gas being accreted onto high-z galaxies is a possible explanation for the initial turbulent state of the disk, while processes such as clump formation/interaction and star formation feedback are suggested to maintain turbulence (e.g., Dekel et al. 2009;Ceverino et al. 2010;Krumholz & Burkhart 2016;Krumholz et al. 2018). While the relation between sSFRs and the clumpy fraction supports the picture of in-situ clump formation, we note that these results can also be explained through ex-situ mechanisms such as frequent gas-rich minor mergers. Such events can induce local starbursts and produce clumpy morphologies by increase the gas density locally, while decreasing Q to be less than one.

CONCLUSION
We use finite resolution deconvolution (FIREDEC) to resolve galaxies at the kiloparsec-scale in ground-based images of the COSMOS field. The field covers an area of ∼2 square degrees. A total of 20,185 star-forming galaxies at 0.5 < z < 2 are selected to satisfy our criteria of a color, stellar mass, and S/N cut. We deconvolve galaxies to a target angular resolution of 0.3", corresponding to a physical size of ∼2.4 kpc for the redshift range considered. Listed below are a summary of our conclusions: • We compare the multi-wavelength deconvolved images with the corresponding HST imaging (where available) in Figures 5 and 6, finding that resolved clumpy structures in deconvolved images are broadly consistent with those observed in the HST/ACS images. Some deconvolved images may have low S/N, causing deconvolved noise look similar to clumpy structures. However, noise is not correlated across the photometric filters, so the correct clump properties can still be inferred from SEDs modeling.
• Section 6.1 shows the relation between the fraction of clumpy star-forming galaxies, redshift and stellar mass. We find f Urest clump evolves with redshift, increasing from 30% to 50% between z ∼ 0.5−2 (top panel of Fig. 9). The clumpy fractions are found to be dependent on stellar mass where massive galaxies have lower clumpy fractions, regardless of the ways clumps are detected (Fig. 8). Galaxies also have higher f Urest clump compared to f Vrest clump or f mass clump , while f Vrest clump and f mass clump are similar. This is expected since the rest-frame UV radiation come from younger stellar populations, while the rest-frame optical trace older stellar populations and stellar masses within galaxies.
• In Section 6.2, we compare our clumpy fractions to the clumpy fractions derived from HST images, finding that our measurements are consistent with those from the literature (Fig. 9). Our measurements have smaller uncertainties due to a larger sample size. We also verify that there are no systematic differences between the fractions measured from the deconvolved and HST images.
• We find a larger fraction of clumpy galaxies located along the upper envelope of the star-forming main sequence as shown in Figure 10. Similarly, at fixed stellar masses, clumpy galaxies have higher SFRs compared to the average star-forming galaxies. If the sSFR is correlated with gas density, this result supports the picture of in-situ clump formation, where continuous gas accretion leads to disk instability and fragmentation.
• We trace the progenitors of two galaxy populations at z ∼ 0.7 in Section 7. Figure 13 shows the clumpy fractions decreasing as galaxies grow in mass. This is consistent with a down-sizing picture, where clumpy morphologies are found mostly in low-mass systems in the local universe.
• Following the progenitors of a galaxy population with log(M * /M ) ∼ 10.9 at z ∼ 0.7, we find that the mass residing in clumps account for ∼5% of the total mass budget of galaxies at z ∼ 1.7. The fractional mass contribution from clumps decreases to ∼1% at z ∼ 0.7. However, these galaxies are growing in mass, so the total mass within their clumps may actually remain constant through redshift.
• When comparing two different galaxy populations of similar mass, log(M * /M ) ∼ 10.2, but at different redshift (z ∼ 0.7 and z ∼ 1.7), we find that clumps in the high-z population have a higher fractional mass contribution population.
Considering the limitations of image deconvolution, one of the notable results is that our measurements are largely consistent with those from previous HST studies. In addition, we overcome the issues that previous studies (and future HST/JWST studies) have in term of limited filter coverage and sample size. The main advantage of image deconvolution is the access to a rich sample of galaxies from wide-field and multi-wavelength surveys such as COSMOS. Indeed, this study already presents that largest sample of distant galaxies with resolved imaging. It could also be very beneficial to use image deconvolution in synergy with future ground and space-based surveys (e.g., Large Synoptic Survey Telescope and Euclid Space Telescope). Considering the needs for resolved observations and given what deconvolution can provide, it is inviting to use deconvolution to address the topics of galaxy formation and evolution.

ACKNOWLEDGMENTS
We would like to thank the anonymous referee for the constructive feedback. V. Sok would also like to thank Frédéric Courbin for the interesting discussions on image deconvolution. The computations were performed on the Lesta computer cluster at the University of Geneva.
Software: EAZY (Brammer et al. 2008), FAST/IDL (Kriek et al. 2009 The first term is a chi-squared term that relates the solution (B) to the observed image (D) given a PSF (P ). The second term is a regularization term weighted by a Lagrange parameter, λ. Since a similar cost function is used to numerically model the PSF of a star (see Section 3.1), there are two different and important Lagrange parameters for deconvolution. For here on, we denote the Lagrange parameter as α when referring to the deconvolution of D to obtain B, and β when referring to the numerical fit of the star. The purpose of the second term in Equation A1 is to mitigate noise enhancement during deconvolution. Since the regularization term contains information of high frequency components for B, low Lagrange parameters mean that there is minimal penalization on noise that arise during deconvolution, which can result in "noisy" deconvolved images. On the other hand, a high value of α and β will result in smoother images as noise (and subsequently structures) are harshly penalized. In the following sections, we discuss how the best value for α and β are determined.
A.1. Numerical Fit Figure 14 shows an example of how different values of β affect the numerical fit of a star observed in the J filter. Row (a) shows the residual of the analytic fit. There is a halo-like structure that is not properly fitted by the analytic model. The goal of the numerical fit is capture this structure. Row (b) shows the residual of the star after subtracting both the analytic and numerical fits. We find that the halo-like structure is still present in these residual maps for higher values of β (β > 100), but disappear at lower β values. While the residuals of row (b) imply that a low value of β should always be used, the numerical fit itself (i.e., row c) is becoming noisier toward low values of β. This is shown in row (d) where we take the difference between the numerical fits. It is apparent that the pixel-to-pixel variation increases toward lower β. Therefore, only a range of β allows the numerical fit to maximally capture the halo-like structure in the analytic residual, while adding minimal noise to the numerical fit. The best value of β for the example in Figure 14 is β = 1. This is quantifiable by comparing the root-mean-squared errors (RMS) within row (d) of Figure 14 as a function of β. We first calculate the pixel RMS to quantify pixel-to-pixel variations in the residual maps of row (d). The pixel RMS is taken to be the standard deviation of the histogram of pixel intensities. To quantify the large-scale structures within the residual maps, we calculate the RMS of flux within apertures, where the aperture diameter is taken to be the target resolution at 0.3". These apertures are placed randomly within the residual maps of row (d). A good β value should saddle the two RMS. Figure 15 shows both RMS for different filters. We mark the best β for each filter as the vertical dotted line.

A.2. Deconvolution of Galaxies
The regularization term in Equation A1 is weighted by a different Lagrange parameter (α) when deconvolving images of galaxies. Figure 16 shows an example of how different α values affect deconvolution. Row (a) shows different deconvolution of the same galaxy, but with α increasing by a factor of ten from left to right. Row (b) shows the differences between each deconvolved image. Since the deconvolved images are obtained using the same parameters (with the exception of α), the residual maps illustrate how α affects deconvolution. In particular, the deconvolved images mainly consist of noise at lower α (α 0.01), whereas smoother images are obtained with higher values. The fact that the residual maps in row (b) have a uniform noise distribution demonstrates that a solution can converge by lowering α. However, lowering the value of α further introduces more noise in the deconvolved image, and the solution becomes degenerate.
We quantify the convergence by comparing the aperture noise within the foreground and background of the residual maps. The dividing line between the foreground and background of the galaxy is defined by the segmentation map of the galaxy. As the ratio between the two RMS converges to unity, it indicates that no new information are being fitted by deconvolution. On the other hand, the pixel RMS within the background of the deconvolved image is increasing as α decreases. This is summarized in Figure 17 for a few galaxies in the z + -band. A good α should saddle the two curves. As shown in the example of Figure 17, α = 1 provides a good deconvolution model.

B. DEPENDENCE OF THE FRACTIONAL LIGHT CONTRIBUTION OF CLUMPS ON THE SSFR
In Section 6.3, a slight dependency was observed between the fraction of clumpy galaxies and ∆ log(sSFR). However, classifying galaxies based on a luminosity cut-off (e.g., where clumpy galaxies are defined to have ≥ 8% of their luminosity in clumps) can blur out the relation between the abundance of clumps and ∆ log(sSFR). An alternative way to test this relation is to use the fractional UV contribution of clumps. Figure 18 shows two plots similar to Figures 10 and 12, but the median C UV is used instead of the clumpy fraction. C UV represents the fraction of UV light coming from clumps to the total UV light of the host galaxy. Comparing the left panel of Figure 18 to Figure 10 shows that the clumpy fraction is correlated to the median C UV . This is not surprising as clumpy galaxies are defined to have C UV ≥ 0.08. Galaxies with higher C UV are therefore generally found above the main-sequence relation.
The right panel of Figure 18 shows the median C UV as a function of ∆ log(sSFR). The darker shaded regions show the standard error of the median, which measures the discrepancy between our sample's median and the true population's median. This is calculated as σ/ √ n, where σ is the dispersion of C UV measured relative to the median value and n is the number of galaxies in each bin. In general, we find that the median C UV increases with ∆ log(sSFR). However, the median C UV decreases at ∆ log(sSFR) > 0.5, indicating that clumps within galaxies contribute less to the total UV luminosity. A slope of 0.04 ± 0.01 and 0.02 ± 0.01 is respectively obtained for the high-z and low-z bin by fitting a linear line. We test the significance of the slopes using a t-test statistic, where a t-score of 4 and 2 is obtained for the high-z and low-z bin. These t-scores correspond to p-values that are <0.05. While this result formally suggests that a correlation exists between C UV and ∆ log(sSFR), we note that the slopes are still shallow and that C UV is highly scattered. The 1-σ dispersion of C UV is illustrated as the lighter shaded region. The numerical fit as a function of β. Row (a) shows the residual of the star after subtracting the analytic fit. The PSF wings, seen as the halo-like structure, cannot be modelled by the analytic fit alone and require a numerical fit. Row (b) shows the residual of the star after subtracting both the analytic and numerical fits. At higher β, the numerical fit fails to fully capture the PSF wing, whereas lower β values allow for a more accurate description. Row (c) and (d) show the numerical fit and the difference between each numerical fit, respectively. While lower β values still allows for the PSF wings to be modelled, the numerical fit is generally noisier as a result of no noise suppression in the cost function. Figure 15. The RMS as a function for β for different filter. We plot the mean pixel and aperture RMS as the blue and red points. The error is taken as the standard deviation for all the stars in COSMOS. A good β should saddle the two RMS. The vertical dotted line indicates the value we used for each filter. Row (a) shows the deconvolved z + -band image of a galaxy, with α increasing by a factor of ten from left to right. In general, lower α result in noisier images. Row (b) shows the differences between each deconvolved image. The residuals within these maps become more uniform at lower values of α, indicating that no new information is being fitted by deconvolution. It should also be noted that the noise in the deconvolved image are enhanced as α is decreased. A good α should minimize the noise in the deconvolved image (row a), but also allow structures to be maximally fitted. The latter case is illustrated in row (b) where the residuals map become uniform between α = 0.1 and α = 1. Figure 17. The left panel shows the noise for a few deconvolved z + -band images as a function of α, while the right panel shows the ratio of the aperture noise within the foreground and background of the residual maps in Figure 16. The foreground and background of the galaxy is defined by the edge of the segmentation map for the galaxy. The grey dotted lines are values for different galaxies. The means are also plotted, where the shaded region represents the standard deviation. In general, lowering α allows finer structures to be fitted, but can cause noise to increase in the background. A good α should saddle the two curves. For the z + filter, the black dashed line at α = 1 represents a good value.  Figures 10 and 12, but we use the median CUV as opposed to the clumpy fraction. We generally find similar trends based on the clumpy fraction and CUV. In the left panel, the median CUV is typically higher for galaxies above the main sequence. In the right panel, the median CUV is observed to increase with ∆ log(sSFR) for both redshift bins. However, the line of best fit indicates that the slope of the relation is relatively shallow. The lighter shaded region shows the 1-σ dispersion of CUV, while the darker shaded region shows the standard error of the median (σ/ √ n).