© The Authors 2025

1 Introduction

In the standard Lambda cold dark matter (ΛCDM) paradigm, galaxies form as gas accretes and cools within dark matter halos, which in turn shapes their sizes and morphologies (e.g., Fall & Efstathiou 1980; Mo et al. 1998; Somerville et al. 2018). Over their lifetimes, galaxies undergo various physical processes that regulate or disrupt their growth, such as star formation, feedback from stars and active galactic nuclei (AGNs), mergers, and gas accretion, all of which leave imprints on their structural properties. Consequently, studying galaxy morphology across cosmic time is key to understanding the physical transformations that accompany galaxy evolution, as well as the interplay between galaxies and their host dark matter halos (e.g., Conselice 2014, for a review).

Historically, galaxy morphology has been studied through a range of techniques that have evolved alongside observational capabilities. Early classifications followed the Hubble-de Vaucouleurs visual classification scheme (Hubble 1936; de Vaucouleurs 1959). These visual schemes relied on manual inspection of photographic plates to identify structural features such as spiral arms, bars, and bulges, forming the basis for the first physical interpretations of galaxy evolution. With the advent of digital imaging, quantitative approaches became more common. Parametric methods, such as fitting analytic functions to galaxy light profiles, became standard, beginning with the de Vaucouleurs r^1/4 law for elliptical galaxies and generalized by the Sérsic profile (Sérsic 1963), which remains widely used to characterize galaxy structure through parameters such as concentration, effective radius, and ellipticity. In parallel, nonparametric methods were developed to capture structural diversity without assuming specific functional forms, using parameters such as concentration, asymmetry, clumpiness (CAS), Gini coefficient, and M₂₀ to quantify galaxy morphology in a model-independent way (Conselice 2003; Lotz et al. 2004).

A series of large surveys from both ground- and space-based observatories has systematically advanced our understanding of galaxy morphology across cosmic time. In the local Universe, the Sloan Digital Sky Survey (SDSS) enabled statistical studies of galaxy structure, revealing, for instance, that the size distribution at fixed luminosity follows a log-normal form at z ≲ 0.3 (Shen et al. 2003). The Galaxy Zoo project leveraged SDSS imaging to produce visual classifications for hundreds of thousands of galaxies (Lintott et al. 2008), enabling studies of morphological diversity and its connection to environment and star formation (e.g., Schawinski et al. 2014). At intermediate redshifts, the Hubble Space Telescope (HST) played a transformative role through legacy programs such as GEMS (Rix et al. 2004), COSMOS (Koekemoer et al. 2007), and CANDELS (Grogin et al. 2011; Koekemoer et al. 2011), which enabled high-resolution imaging and morphological measurements out to z ~ 3. These datasets allowed for extensive morphological analyses, including evidence that quiescent galaxies are more compact, spheroidal, and concentrated than star-forming galaxies at fixed mass (e.g., Daddi et al. 2005; Trujillo et al. 2007; Barro et al. 2013, 2017; Schreiber et al. 2016), and that they undergo size growth primarily via dry mergers (e.g., van der Wel et al. 2014; van Dokkum et al. 2015; Whitaker et al. 2017). However, the mechanisms that drive galaxy quenching remain debated, for example, whether morphological transformation and the cessation of star formation are causally linked, or whether they occur independently as consequences of other processes (e.g., Tacchella et al. 2015). Disentangling this connection requires large, statistically robust samples of galaxies with reliable morphological measurements spanning a wide redshift range, in order to trace structural evolution across cosmic time. Furthermore, at higher redshifts (z ≳ 3), HST imaging becomes limited to the rest-frame ultraviolet, which predominantly traces young stellar populations and can bias morphological interpretations (Shibuya et al. 2015). Nevertheless, these surveys laid the groundwork for the construction of large morphological catalogs (e.g., York et al. 2000; van der Wel et al. 2012; Kawinwanichakij et al. 2021), providing critical insights into the growth and structural transformation of galaxies over more than 10 billion years of cosmic history.

The James Webb Space Telescope (JWST) has opened a new era in the study of galaxy morphology, enabling rest-frame optical measurements out to z ≳ 7 with unprecedented depth and resolution. Studies based on JWST observations have revealed the evolution of galaxy morphology and size for both star-forming and quiescent galaxies out to unprecedented redshifts (e.g., Varadaraj et al. 2024; Ormerod et al. 2024; Allen et al. 2025; Matharu et al. 2024; Wright et al. 2024; Ito et al. 2024; Yang et al. 2025). These studies show that quiescent galaxies remain structurally distinct from star-forming ones, with smaller effective radii and higher Sérsic indices, even at early times. Nevertheless, most JWST morphological analyses to date have focused on specific galaxy selections, such as mass- or color-selected subsamples, or have been limited to individual fields. As a result, there is no consistent, flux-limited morphological catalog across different major JWST fields and redshifts. This limits our ability to conduct comparative and statistical studies of galaxy structure, particularly for rarer populations such as quiescent galaxies at the full depth and breadth enabled by JWST.

In this work, we build one of the largest and most uniform morphological catalogs to date using JWST imaging, enabling consistent structural measurements across multiple fields and redshifts. We utilize publicly available mosaics from the DAWN JWST Archive (DJA, Valentino et al. 2023), which compiles deep NIRCam imaging from major extragalactic surveys, including CEERS, PRIMER (UDS and COSMOS), and GOODS. These are accompanied by matched photometric red-shifts and physical parameters. To measure galaxy morphology, we carry out two-dimensional surface brightness profile fitting using SOURCEXTRACTOR++ (Bertin et al. 2020; Kümmel et al. 2020), a modern, scalable tool optimized for catalog-level model fitting in large multiband datasets. We fit each galaxy with both a single-component Sérsic profile (Sérsic 1963) and a two-component Bulge+Disk (B+D) model consisting of an exponential disk (n_S = 1) and a de Vaucouleurs bulge (n_S = 4). This dual-model approach allows us to probe structural diversity more flexibly and to derive key parameters, such as the Sérsic index (n_S ), the effective radius, the axis ratio, and the bulge-to-total ratio (B/T). We apply this modeling to all sources above a flux and signal-to-noise threshold, yielding a homogeneous morphological catalog spanning approximately 340 000 sources across 0 < z < 12. The aim of this catalog is to enable the community to carry out detailed statistical studies of galaxy morphology evolution throughout cosmic time.

This paper is organized as follows. The data and point spread function (PSF) reconstruction are presented in Section 2. In Section 3, we describe our method for measuring the morphology and sizes of galaxies across such a large field. In Section 4, we present our results on the correlation between morphology and the UVJ diagram, as well as the size evolution in the Sérsic and B+D models. We summarize our results in Section 5. We adopt a Planck Collaboration VI (2020) ΛCDM cosmology with H₀ = 67.4 km s⁻¹ Mpc⁻¹, Ω_m,0 = 0.315. All magnitudes are expressed in the AB system (Oke 1974).

2 Data

This work uses NIRCam images (Rieke et al. 2005) from public surveys of the JWST, processed as part of the DAWN JWST Archive (DJA), mosaic release v7. The DJA is an online repository containing reduced images, photometric catalogs, and spectroscopic data from public JWST data, and is described in more detail in Valentino et al. (2023) and Heintz et al. (2025). In this work, we carry out morphological measurements in several major extragalactic surveys by JWST, covering a total area of ~500 arcmin². We focus on the following fields:

1. EGS from CEERS (DD-ERS 1345 Finkelstein et al. 2023);

2. GOODS from JADES (GTO 1180, 1181, 1210, 1287 Eisenstein et al. 2023);

3. UDS from PRIMER UDS (GO 1837 Donnan et al. 2024);

4. COSMOS from PRIMER COSMOS (GO 1837 Donnan et al. 2024).

Table 1 presents the area and photometric bands used in these fields and their 5σ depths (computed from empty apertures by Weibel et al. 2024).

2.1 Images

The images used in this work are drawn from the DJA mosaic release v7. These are mosaics of the whole or split fields: ceers-full-grizli-v7.2, gds-grizli-v7.2, gdn-grizli-v7.3, primer-uds-north-grizli-v7.2, primer-uds-south-grizli-v7.2, primer-cosmos-west-grizli-v7.0, and primer-cosmos- east-grizli-v7.0.

The DJA provides reduced and calibrated images, in photo-metric units (10 nJy), which we used to run SOURCEXTRACTOR++, alongside the weight maps that provide the inverse variance distribution. For source detection, we used the inverse-variance weighted stack of the long-wavelength filters (specified with ir in place of the filter in the image name). We did not apply any PSF convolution, since SOURCEXTRACTOR++ convolves the source models with the corresponding filter PSF (see Sect. 2.3).

2.2 Catalogs

The DJA also provides photometric and photo-z catalogs produced with SEXTRACTOR (Bertin & Arnouts 1996) and EAZY (Brammer et al. 2008). We used the DJA SEXTRACTOR catalogs as our principal list of sources for each field. Although we performed a separate detection with SOURCEXTRACTOR++ (Sect. 3.1), we cross-matched with the DJA sources to ensure model-fitting-based measurements for the same sources. This approach also aligns with the DJA’s provision of photo-z and physical parameters from spectral energy distribution (SED) fittings with EAZY on the SEXTRACTOR catalogs. In this work, we used the EAZY output to investigate morphology evolution as a function of redshift and galaxy type (Sect. 4).

2.3 Point spread function reconstruction

Accurate model fitting requires precise characterization of the instrument’s PSF. The PSF results from light diffraction at the aperture of the telescope and defines the resolution limit of JWST. We empirically modeled the PSF from the final mosaics using PSFEX (Bertin 2011), which builds a model for the PSF by fitting a set of basis functions to point sources provided in an input catalog. Recent work by Berman et al. (2024) shows that PSFEX accurately models the NIRCam PSF on mosaics.

To select point sources, we developed a method inspired by Leauthaud et al. (2007). This method uses the MU_MAX and MAG_AUTO measurements from SEXTRACTOR. The parameter MU_MAX represents the surface brightness of the brightest pixel of a source, while MAG_AUTO denotes the Kron-like elliptical aperture magnitude. In the MU_MAX/MAG_AUTO plane (shown in Fig. 1), the point-like sources follow a line of slope equal to one (referred to in this paper as the “starline”), and extended sources form a distinctive cloud above it. We added thresholds for MAG_AUTO: a minimal value to avoid selecting saturated sources that would show a truncated PSF, and a maximal value to avoid including sources in the extended sources cloud.

Detection of the starline was performed using linear regression. However, to avoid bias from the extended sources cloud, we first applied a clustering algorithm to detect this point cloud and remove it. We used DBSCAN (Ester et al. 1996), implemented in the SCIKIT-LEARN Python package. This algorithm identifies cores of high density (the extended sources cloud) while excluding points considered as noise (the starline). We then applied a threshold on the MU_MAX – MAG_AUTO value to further remove outliers from the extended sources cloud. Finally, we detected the starline using the RANSAC linear regression algorithm (Fischler & Bolles 1981), implemented in SCIKIT-LEARN. The threshold and width around the starline used to select point-like sources were chosen empirically on the basis of additional tests. Fig. 1 shows the selection of point-like sources in the MU_MAX/MAG_AUTO plane.

The point-like source selection was expected to physically yield the same sample, independent of the band used. However, slight differences arose due to noise or differences in the threshold used for the MU_MAX – MAG_AUTO value. We chose the F200W band to select point-like sources, as it provided the most consistent starline detection across fields. Finally, we visually inspected a randomly selected sample (N ~ 100) to verify that the selected sources were point-like. Figure A.1 shows examples of the PSFs in the GOODS-S field.

Fig. 1 Point-like sources selection for PSF estimation. The plot shows the distribution (log scale) of sources from the GOODS-S field (F200W) in the MU_MAX/MAG_AUTO plane. The red box shows the starline selection of point-like sources, with the starline indicated by red dots.

3 Morphological measurements

3.1 SOURCEXTRACTOR + +

In this work, we present one of the largest catalogs of galaxy morphologies measured from JWST images. To accomplish this, we used SOURCEXTRACTOR++ (hereafter abbreviated as SE++; Bertin et al. 2020; Kümmel et al. 2020), the successor of the classic SEXTRACTOR Bertin & Arnouts (1996), which performs brightness profile fitting to multiple sources in the multiband images. The principal advantage of SE++ is that it is optimized exactly for the purpose of large catalog creation – it can simultaneously fit all sources in arbitrarily large mosaics and in different bands, without the need for prior sample selection, cutout creation, and masking. Additionally, its flexible model fitting allows the user to define any model using a simple Python configuration file. In the following, we describe our models and catalog-making procedure.

3.2 Brightness profile models

To measure galaxy morphology, we modeled the brightness profile of each source. The SOURCEXTRACTOR++ tool provides full control over the model. To maximize scientific applications, we used the two most common models:

1. The Sérsic model (Sérsic 1963), parameterized by a single shape parameter, the Sérsic index n_s ; the effective radius R_e; the axis ratio (a/b); the position angle θ; and the total flux f_tot. In this work, we fit for the two components of the ellipticity, e₁ and e₂, which are directly linked to (a/b) and θ. The priors used for these parameters are presented in Fig. B.1.

2. The Bulge + Disk (B+D) model, a composite of an exponential disk (n_s = 1) and a de Vaucouleurs bulge (n_s = 4) Sérsic profile. It is parametrized by the effective radii of the disk, R_disk, and bulge, R_buige; the axis ratios (a/b)_D and (a/b)_B; a common angle θ_BD; the total flux of both components, f_BD,tot; and the bulge-to-total ratio B/T = f_B,tot/f_BD,tot. The B/T prior is a bell curve ranging from 5 × 10⁻⁵ to 1, with a mean and spread that increase as a function of wavelength.

The B+D model generally provides a better fit than the Sérsic model, as it more accurately describes the central region of galaxies, particularly the presence of a bulge or a dimmer center compared to a Sérsic profile. However, because it has more parameters, it is more computationally intensive and can, in some cases, lead to degeneracy in the model parameters, especially for fainter sources where the two components may not be distinguishable.

The SE++ tool performs model fitting simultaneously on multiple images, typically across different bands. In this work, we constrain the morphological parameters (R_eff, n_s, (a/b), and θ) to be identical across all bands for a given source. This approach yields final values that represent a weighted average over the entire wavelength range used (0.8–5 µm). The SE++ configuration file defining the models is available on GitHub¹.

3.3 Tiling

It is theoretically possible to run SOURCEXTRACTOR++ directly on the full mosaics from the DJA. However, because of speed and computing power limitations, we chose to tile the full images. This approach makes sources modeling faster, thanks to parallelization. We divided the mosaics into 2′ × 2′ tiles, with 0.5′ overlap, to ensure that sources fragmented by one tile were fully present in the neighboring tile.

3.3.1 Catalog merging

To merge the sub-catalogs, we iterated over each tile to append them and create a catalog covering the whole field. At each iteration step, we cross-matched with the previously appended sub-catalogs and discarded the matched duplicates from the overlap regions. For sources that are matched (in the overlap region), we chose the one with the smallest uncertainty in the F200W magnitude, as measured by SOURCEXTRACTOR++. As the model fitting was performed independently for the Sérsic and B+D models, this step was performed separately to produce one catalog per model. The F200W band was chosen because it has better resolution compared to the LW channel, and galaxies are generally brighter in it than in other bands of the SW channel.

Finally, we cross-matched the two model-fitting catalogs with the DJA catalogs. This step allowed us to remove false detections (which were frequent near the edges of images), retain the same list of sources, and provide additional morphological measurements to the DJA catalogs. We chose a threshold of 0.3″ as the cutting distance to validate a match. This threshold was set manually using the histogram of angular distances produced by the cross-match. This value corresponded to 5 pixels in the LW channel and 10 pixels in the SW channel of NIR-Cam, which was acceptable and could be the result of differences between the source centroids estimated in SEXTRACTOR and SOURCEXTRACTOR++.

3.3.2 Merging of the model and residual images

The SE++ tool produces model and residual images. We merged only the model sub-images and generated a mosaic residual image afterward. To merge the sub-images, we used the reproject_and_coadd function from the reproject Python package (Robitaille 2018), which reprojected and co-added the images on a frame specified by the World Coordinate System (WCS) of the native DJA images. This ensured that the merged full model images had the same pixel scale, center, and orientation as the DJA images. We used this to generate mosaic residual images by subtracting the mosaic model images from the source DJA images.

3.4 Flagging and completeness

As mentioned previously, our aim was to add morphological measurements to the DJA catalogs. Therefore, by cross-matching, we retained the same number of rows (sources) in our catalogs as in the DJA catalogs. Furthermore, we added a flag keyword with four possible values as follows:

• 0 : the source was not fitted (no morphological data).

• 1 : a potential artifact occurred during model fitting.

• 2 : fitting was performed successfully. This value gives science-ready data.

• 3 : the source has S/N < 3 or a magnitude fainter than the 5σ depth of its corresponding survey. The S/N was computed by SOURCEXTRACTOR++ on the detection image.

The distribution of the flag values is shown in Fig. 2. Seventy to eighty percent of the sources were successfully fitted (flag=2). This indicates that the SOURCEXTRACTOR++ minimization algorithm converged and the parameter values were not equal to the initial values nor at the minimum or maximum of the allowed range. Sources without morphological data (flag=0) were not detected in the SOURCEXTRACTOR++ detection step and are likely artifacts or faint sources near the image depth, which would be difficult to fit nonetheless. Potential artifacts (flag=1) are sources that have been fitted, but with some parameters remaining at their initial values or that strayed to the extreme ends of the allowed parameter range. Typically, this indicates a poor fit; therefore, we flagged such sources. For the Sérsic model, this includes sources with (a/b) > 0.99, |n_S − 0.36| < 10⁻⁴, n_S > 8.35, or n_S < 0.301. For the B+D model, this includes (a/b)_D > 0.9999, (a/b)_D < 0.10001, or |(a/b)_D − 0.5| < 10⁻⁵. These values were chosen manually by identifying artifacts in the parameter distribution and by visual inspection. These sources should be handled with care by the user.

Since the list of detected and measured sources by SOURCEXTRACTOR++ is slightly different from the initial DJA list, we quantified the completeness of our morphological catalogs. Table 2 presents the number of galaxies in each field used in the DJA, and the number of successfully fitted galaxies (flag=2) in each model (Sérsic and B+D) from our work.

The completeness was measured relative to the number of galaxies in the DJA catalogs, using the same quality cut on S/N and magnitude as S/N >3 and mag <5σ depth (flag=3) on both. In total, our catalog consists of 342 892 sources with reliable model fitting and science-ready, making it one of the largest morphological catalogs based on JWST observations.

A total of 80.4% of the DJA S/N >3 and mag <5σ depth sources were successfully fitted with both a Sérsic and a B+D model. The completeness is higher for the B+D model than for the Sérsic model. This is likely because the B+D model is a better description for some sources whose Sérsic model parameters tend toward the minimum or maximum allowed values (therefore classified as flag=1). Finally, we analyzed the distribution of certain physical parameters (F277W magnitude, z_phot, mass, and Kron radius) for the non-detected sources (flag=0) and potential artifacts (flag=1). We did not find any significant correlation between these parameters, indicating that the incompleteness due to flag=1 is not biased with magnitude, redshift or stellar mass.

Therefore, we consider our morphological catalogs to have a relatively high completeness compared to the DJA, and, as such, they provide highly valuable information for studying the morphology of galaxies and its evolution through cosmic time, as demonstrated in Sect. 4.

Fig. 2 Histogram of the flag values for the full catalogs. The percentages are relative to the total number of sources in the DJA catalogs. Two distinctive histograms for the Sérsic and Bulge+Disk models are shown because their fitting success rates differ. The figure shows a vast majority (70–80%) of sources are correctly fitted and provides science-ready morphological data.

3.5 Comparisons with previous work

To validate our measurements, we compared them with previous morphological catalogs in the same fields, particularly van der Wel et al. (2012). These catalogs contain Sérsic modeling performed using Hubble Space Telescope (HST) observations in the same fields used in this work: GOODS, EGS, UDS, COS-MOS. By cross-matching them with our catalog and selecting only early-type galaxies (see Eq. (2)), we find 3263 matches with S/N >10.

The size and morphology of galaxies are known to differ in different wavelength ranges. To ensure that we compared morphologies measured at similar observer-frame wavelengths, we used the measurements in F160W from van der Wel et al. (2012). However, since our measurements correspond to the averaged morphology over NIRCam’s wavelength range (~1–5 µm), we scaled the effective radii R_eff in van der Wel et al. (2012) to 2.5 µm using the following scaling relation: $$\frac{\Delta \log R_{\text{eff}}}{\Delta \log \lambda }=-0.25,$$(1)

which was found by van der Wel et al. (2012) to describe the wavelength dependence of R_eff for early-type galaxies. We note that a more complex version of this relation has been proposed by the same authors, which integrates the effect of redshift and stellar mass. However, we adopted this simpler and more general relation for comparison purposes.

Figure 3 shows the one-to-one comparisons between the two works for R_eff (top middle panel) and the Sérsic index n_S (top right panel), as well as the histograms of R_eff, n_S, and the axis ratio (a/b) (bottom panels). In general, there is good agreement in the effective radii, with a slight trend toward larger radii from our measurements. One reason for the slight offset could be the scaling using Eq. (1) that we applied. Without this scaling, the trend inverses and our effective radii are smaller than those of van der Wel et al. (2012).

Fig. 3 Comparison between Sérsic model fitting in our work and van der Wel et al. (2012) for UVJ-selected quiescent galaxies. The figure compares effective radii Reff and Sérsic index n. We corrected the Reff values from van der Wel et al. (2012), originally measured at 1.6 µm, to 2.5 µm using Eq. (1), in order to match the average wavelength of our JWST images. The color of the hexbin cells indicates the counts in each cell in log scale. Gray dotted lines show identical values between the two works. This figure includes 3263 matched galaxies.

4 Morphological evolution with redshift and galaxy type

To demonstrate the scientific application of our morphological catalogs, we investigated statistical distributions of several morphology indicators as a function of redshift and galaxy type. Our measurements cover a wide range of redshifts, enabling studies of morphological evolution across cosmic time. For this analysis, we used the photo-z, rest-frame colors, and other physical parameters from the DJA EAZY catalogs. In particular, we investigated morphological evolution as a function of galaxy type (quiescent vs. star-forming).

4.1 UVJ diagram as a function of morphology and redshift

To investigate the relationship between galaxy morphology, type, and redshift, we used the UVJ color diagram to select star-forming and quiescent galaxies (Williams et al. 2009). In principle, the UVJ selection criteria for quiescent galaxies change with redshift (Williams et al. 2009). However, for our qualitative analysis, we adopted a redshift-independent selection corresponding to the 0.7 < z < 1.3 range used in Schreiber et al. (2016), who also investigated the morphology dependence of the UVJ diagram from HST observations: $$UV{J}_{\text{quiescent}}\{\begin{array}{c}U-V>1.3,\text{\hspace{0.17em}and}\\ V-J<1.6,\text{and}\\ U-V>0.88(V-J)+0.49.\end{array}$$(2)

In addition, we focused our analysis on log M_★/M_⊙ > 10, flag=2, and S/N >10, resulting in 13 685 galaxies. This ensured that our sample had sufficiently high S/N and robust morphological estimates from the model fitting.

First, we investigate the location of galaxies in the UVJ diagram as a function of the Sérsic index, n_S, and redshift. Fig. 4 shows the UVJ diagram in six redshift bins at 0 < z < 6 (in different panels), color coded by the Sérsic index. The dotted lines mark the regions separating star-forming and quiescent galaxies (Eq. (2)). Fig. 4 shows a clear correlation between the Sérsic index and UVJ colors. There is a gradient of the Sérsic index, n_S , in a direction roughly orthogonal to the boundary between quiescent and star-forming galaxies, such that galaxies with a higher Sérsic index preferentially populate the redder and quiescent UVJ region. This trend holds for the different redshift bins. For z > 4, the samples are too small to draw statistically robust conclusions, especially for quiescent galaxies.

Second, because we also used B+D models to fit the galaxies in our catalogs, we could analyze how their B/T values populate the UVJ diagram. For this analysis, we used the B/T measured in the F200W band because it offers better resolution, being in the SW channel. Figure 5 shows UVJ diagrams for three broad redshift ranges and three B/T ranges. The color indicates the B/T value, with the color scale indicated at the top. We also show the contours estimated by kernel density. Fig. 5 shows that bulge-dominated galaxies (B/T > 0.6) preferentially occupy the quiescent region, whereas disk-dominated galaxies (B/T < 0.2) occupy the star-forming region. Intermediate galaxies that show both bulges and disks (0. 2 < B/T < 0. 6) form a diverse population that can be classified as both star-forming and quiescent based on their UVJ colors. However, there is a trend showing that galaxies with a higher B/T preferentially occupy the quiescent UVJ region. At z > 3 there are more galaxies with higher B/T that are star-forming than quiescent. This results from the relative rarity of quiescent galaxies at these epochs, as well as from the presence of a population of compact star-forming galaxies (e.g., blue nuggets, Barro et al. 2013; Dekel & Burkert 2014, which we discuss further in the next part of this section). These qualitative trends from the independently fitted B+D model are in good agreement with those from the Sérsic model. This consistency attests to the robustness of the morphological measurements.

To further investigate the correlation between bulge- and disk-dominated morphology and star formation activity, we defined two classes based on the two independent morphology estimates, n_S and B/T, as follows:

1. Bulge-dominated galaxies: n_S > 1 and B/ T > 0. 5.

2. Disk-dominated galaxies: n_S < 4 and B/ T < 0. 5.

With this, we attempted to provide a robust binary classification of whether a galaxy is disk- or bulge-dominated and explore how they populate the UVJ color space.

In Fig. 6, we show the UVJ diagrams in the same red-shift ranges as in Fig. 4, with color indicating the two different classes: red for bulge-dominated galaxies and blue for disk-dominated galaxies. We plotted kernel density contour lines for each class. This classification is consistent with the broad conclusion that quiescent galaxies are bulge-dominated and star-forming galaxies are disk-dominated.

The disk-dominated sample still exhibits a degree of bimodality, with a small population of quiescent galaxies. Indeed, while star formation is typically associated with disk galaxies, quiescent disks do exist. At lower masses, quenching is often environmentally driven, particularly via strangulation (Larson et al. 1980; Moran et al. 2007), ram pressure stripping (Gunn & Gott 1972), and galaxy harassment (Moore et al. 1996, 1998) in dense environments (Peng et al. 2010; Cortese et al. 2021). At higher masses, some fast-rotating quiescent disks can form through a combination of mild dissipative contraction and secular evolution (Toft et al. 2017; D’Eugenio et al. 2024). These quiescent disks may also result from gas exhaustion following a compaction event or mergers that preserve disk kinematics but suppress star formation (Omand et al. 2014; Zolotov et al. 2015).

Bulge-dominated galaxies preferentially occupy the quiescent UVJ region out to z ~ 3, with a small bimodality appearing in the star-forming region at higher redshifts. This indicates that bulges can still be actively forming stars, especially at earlier cosmic times, but by z < 1, most have migrated toward quiescence. This trend may indicate a quenching pathway, in which bulge growth coincides with or precedes quenching, consistent with theoretical and observational evidence for compaction-driven quenching (e.g., Barro et al. 2013, 2014; Tacchella et al. 2015, 2016). The location of bulge-dominated galaxies in the UVJ diagram – toward bluer colors compared to the main disk-dominated population – supports a fast quenching scenario, in which lower-mass, young, star-forming galaxies undergo rapid quenching (t_Q < 2 Gyr e.g., Moutard et al. 2018; Belli et al. 2019). This is consistent with the blue nugget phase, where galaxies experience a burst of compact star formation before quenching into red nuggets (e.g., Dekel & Burkert 2014; Zolotov et al. 2015; Barro et al. 2017). Finally, the lower number of bulge-dominated galaxies in the higher redshift ranges is consistent with the inverted Hubble classification evolution: galaxies form as disks and evolve to bulges (e.g., Quilley & de Lapparent 2022).

This qualitative analysis of the correlation between quenching and morphology is consistent with the current picture from both theory and observations. Importantly, by providing morphological measurements from JWST for such a large sample, our work paves the way for more in-depth and quantitative population studies that can unveil the details of galaxy quenching and the accompanying morphological transformations.

Fig. 4 Distribution of Sérsic indices, nS, in the UVJ color space for different redshift ranges, for galaxies with log M★/M⊙ > 10. The color of each hexagonal bin represents the median value of nS, as indicated by the colorbar. The dotted line shows the quiescent vs. star-forming separation using Eq. (2). The quiescent region is predominantly populated by galaxies with high Sérsic indices (nS > 3).

4.2 Size evolution

In this section, we show the evolution of the size with redshift for galaxies classified by their morphology type (namely, bulge-and disk-dominated). We analyzed the radii from the two independent models: the effective radius r_e from the Sérsic model, and the disk effective radius r_disk and bulge effective radius r_bulge from the B+D model. We computed size in kpc by converting the measured angular radius, using the z_phot given in the DJA EAZY catalogs.

In Fig. 7, we show the size evolution with z_phot for the Sérsic effective, disk, and bulge size measures for galaxies with log M_★/M_⊙ > 10, which corresponds to a mean mass of ≈4 × 10¹⁰ M_⊙. We show the trend for all galaxies, disk- and bulge-dominated as defined in Sect. 4.1. We note that our sizes were computed as the average size over the ~1–5 µm range, whereas size-redshift evolution studies typically use the size measured in a given band that probes the same rest-frame optical range over the studied redshift range. This means that our measurements probe the rest-frame optical at z > 0.5; however, this approach may introduce differences when comparing our results with studies that base their size measurements on a single band for a defined redshift range.

Our measurements show that the overall size (r_e, from the Sérsic model) of all galaxies with log M_★/M_⊙ > 10 increases with redshift from about 1 kpc at z ~ 5 to ~2.5 kpc at z ~ 0.5. Disk-dominated galaxies show larger sizes by about 0.1–0.2 kpc compared to the whole sample, and increase with time, while bulge-dominated galaxies are significantly smaller (~1 kpc) with little to no evolution out to z ~ 1.7, but a relatively steep increase in size out to ~2 kpc at z ~ 0.6.

To quantify the size-redshift evolution, we fitted a r_e = R(1 + z)^a model, commonly used in the literature, to our measurements of the Sérsic effective radius. The fit was performed using the least squares method through the lmfit Python package, which returns the uncertainty of each parameter. Fitting the model to all galaxies yields r_e = (3.66 ± 0.04)(1 + z)^{−0.59±0.01} kpc. By restricting the fit to the disk-dominated galaxies population, we obtain r_e = (4.02 ± 0.04)(1 + z)^{−0.61±0.01} kpc. For the bulge-dominated population, the fit is r_e = (3.34 ± 0.10)(1 + z)^{−0.87±0.04} kpc.

We compared our results with the literature based on JWST measurements, namely, Varadaraj et al. (2024), who measured size evolution for Lyman break galaxies (LBGs) selected at log M_★/M_⊙ > 9, and Ormerod et al. (2024), who measured the size evolution for log M_★/M_⊙ > 9.5 galaxies. We also compared with the results from Allen et al. (2025) who measured the size evolution of = 5 × 10¹⁰ M_⊙ galaxies, which was obtained by scaling all galaxies to the same mass using the mass-size relationship. This comparison shows a good consistency with the literature. However, we note that, because of the strong dependence of galaxy size on stellar mass, the different mass limits and sample selection used in the literature complicate the comparison and are likely the reason for the differences.

The middle and right panels of Fig. 7 show the evolution of the disk and bulge sizes of galaxies from the two-component (B+D) model fit. This shows that the disk size increases from ~1.2 kpc to ~3.4 kpc from z ~ 5 to z ~ 0. 5, with the disk-dominated population exhibiting larger disks by approximately 0.2 kpc, as expected. The bulge-dominated population may still host some disk components, which remain small (~1 kpc) out to z ~ 2.5 and increase to ~3.0 kpc by z ~ 0.6. This suggests that the disk components of bulge-dominated galaxies begin to grow their r_disk at later times. Bulges, on the other hand, show little-to-no growth with redshift with the total population exhibiting bulge sizes of about 1 kpc. However, the bulge component of the bulge-dominated population displays a mild increase in r_bulge of ~ 0.8 kpc from z ~ 5 to z ~ 0.6, reaching r_bulge ~ 1.6 kpc. The associated 1 σ dispersion is of comparable magnitude, making it difficult to draw robust conclusions.

We also investigated the size evolution of quiescent and star-forming galaxies selected using their UVJ colors in Fig. 8. Overall, as expected, quiescent galaxies have smaller sizes than star-forming galaxies. The sizes of quiescent galaxies increase with redshift, from ~0.2 kpc for the earliest quiescent galaxies at z ~ 5 to ~2 kpc by z ~ 0.6. Star-forming galaxies are larger and increase in size from ~1.2 kpc to ~2.6 kpc by z ~ 0.5. This shows that the size growth of quiescent galaxies is steeper than that of star-forming galaxies. This size evolution is consistent with previous studies showing that both populations experience growth over time, driven by minor mergers and continuous star formation in star-forming galaxies, whereas quiescent galaxies primarily grow via dry mergers (e.g., van der Wel et al. 2014; Whitaker et al. 2017). While high-redshift quiescent galaxies are compact at formation, often referred to as “red nuggets” (e.g., Damjanov et al. 2009), Fig. 8 indicates that they continue to grow in size down to z ~ 0.5, a trend that can be driven by minor mergers and other accretion processes (e.g., Newman et al. 2012; Belli et al. 2015). However, part of the observed size growth of quiescent galaxies may be attributed to progenitor bias – as new, larger galaxies continue to quench and enter the quiescent population at later times, they increase the average size of the population, making the size evolution appear more dramatic than if individual galaxies were tracked over time (e.g., Barro et al. 2013; Carollo et al. 2013).

The disk sizes of star-forming galaxies increase from ~1.4 kpc at z ~ 5 to 3.4 kpc by z ~ 0.5, consistent with expectations from inside-out growth models, where gas accretion and star formation preferentially occur in the outskirts (e.g., Patel et al. 2012; Morishita et al. 2014; Matharu et al. 2024). The disk sizes of quiescent galaxies remain small and are similar to, or smaller than, their bulges, at least out to z ~ 2, where the disk size increases to about 2.6 kpc. However, the scatter is significantly higher, making it difficult to draw meaningful conclusions about the disk sizes of quiescent galaxies. The bulge sizes for both star-forming and quiescent galaxies exhibit a mild evolution with redshift. Star-forming galaxies grow by about 0.6 kpc, reaching ~1.2 kpc by z ~ 0.5. Quiescent galaxies are slightly more compact and, at high redshift, have bulges of about 0.4 kpc, increasing to about 1.2 kpc by z ~ 0.6. The persistence of small bulge sizes in both populations suggests that bulges reach their final configuration early, while the surrounding disks continue to evolve, particularly in star-forming galaxies. This is consistent with models in which early compaction events, such as mergers or disk instabilities, form a central bulge, after which the fate of the galaxy depends on the availability of fresh gas for continued star formation or quenching mechanisms (e.g., Barro et al. 2017; Tacchella et al. 2018).

Fig. 5 Distribution of B/T in the UVJ color space for different redshift and B/T ranges, for galaxies with log M★/M⊙ > 10. Each point is colored by its B/T value in the F200W band. The dotted line shows the quiescent versus star-forming separation using Eq. (2). The contours indicate 25, 50, and 75% of the density. Bulge-dominated galaxies (high B/T) predominantly occupy the quiescent region, whereas disk-dominated galaxies (low B/T) are generally in the star-forming region of the UVJ space.

Fig. 6 UVJ color diagram for bulge- and disk-dominated log M★/M⊙ > 10 galaxies. The classification combines criteria based on both nS and B/T, estimated independently. Contour lines show kernel density estimates corresponding to 25, 50, and 75% of the density.

Fig. 7 Evolution of galaxy sizes as a function of redshift. The plots show the Sérsic effective radius re (left), the disk effective radius rdisk (center), and the bulge effective radius rbulge (right) as a function of zphot for log M★/M⊙ > 10 galaxies. This corresponds to a mean mass of ≈ 4 × 1010 M⊙. The points with error bars show the mean and 1 σ at a given redshift for all galaxies (black), disk-dominated galaxies (blue), and bulge-dominated galaxies (red). The re evolution is modeled in the left panel for all galaxies, as well as for the bulge- and disk-dominated populations, following re = R(1 + z)a (see Sect. 4.2). We compare the size evolution for all galaxies with previous JWST-based work from Varadaraj et al. (2024), Ormerod et al. (2024), Allen et al. (2025) and Yang et al. (2025).

Fig. 8 Evolution of the sizes of quiescent and star forming galaxies with log M★/M⊙ > 10. The plots show the Sérsic effective radius re, the disk effective radius rdisk, and the bulge effective radius rbulge as a function of zphot. Galaxies are classified as quiescent or star-forming based on their UVJ colors. The points and error bars indicate the mean and 1 σ dispersion in the corresponding zphot bin.

Fig. 9 Evolution of the stellar mass surface density, Σe. Shown are masses and redshifts from the DJA, calculated using EAZY, and effective radii, re, measured by SOURCEXTRACTOR++ in this work. Galaxies are classified as quiescent or star-forming according to their position in the UVJ diagram. The contour lines result from kernel density estimation using seaborn.kdeplot. The data for this figure includes 45 201 galaxies such that log M★/M⊙ > 9.

4.3 Evolution of the stellar mass surface density for quiescent and star-forming galaxies

To investigate the relationship between the compactness of a galaxy and its star formation activity, we calculated the stellar mass surface density, Σ_e, which is the mass contained within the Sérsic effective radius, given by ${\Sigma }_e=M_{★}/2\pi r_e^2$. Fig. 9 shows Σ_e as a function of stellar mass for quiescent and star-forming galaxies of log M_★/M_⊙ > 9, to facilitate comparison with the existing literature, in six redshift bins at 0.5 < z < 6. This figure demonstrates that quiescent galaxies are more compact, with higher surface mass densities than those of star-forming galaxies, and that quiescent galaxies become increasingly dense at earlier times. This finding is in good agreement with previous work, e.g., Barro et al. (2017), which shows this relation out to z ~ 3. Our results indicate that this relation extends out to z ~ 5 and that quiescent galaxies increase in compactness the earlier they form. This is consistent with the observed compactness of some of the earliest quiescent galaxies found by JWST (Carnall et al. 2023; de Graaff et al. 2025; Weibel et al. 2024; Ito et al. 2024; Wright et al. 2024).

5 Conclusions

This work presents a catalog of galaxy morphologies measured from JWST imaging of the major extragalactic surveys CEERS, GOODS, and PRIMER, processed through the DAWN JWST Archive. The catalog contains morphology information for more than 340 000 sources, detected on NIRCam LW stacks, with photometric redshifts extending to z > 10. We fitted each source with two independent models of brightness profiles: a single Sérsic model and a Bulge+Disk model, using SOURCEXTRACTOR++.

To validate our measurements, we compared our results with those from the literature obtained using independent methods and software, and we find good consistency. To demonstrate the scientific application, we used our morphological measurements in combination with the photo-z and physical parameters from the DJA EAZY catalogs. We investigated the relation between galaxy morphology, UVJ colors as an indicator of star formation activity, and redshift. We summarize our results as follows:

• Using UVJ color selection, we find that early-type, bulge-dominated galaxies predominantly reside in the quiescent region, whereas late-type, disk-dominated galaxies exhibit blue, star-forming colors. We also observe a bimodality in these distributions, with some populations of quiescent galaxies that are disk-dominated, as well as star-forming galaxies that are bulge-dominated;

• The two-component fits reveal that low B/T galaxies preferentially occupy the star-forming UVJ region, while high B/T galaxies populate the quiescent region. At z > 3, however, we observe a population of high B/T and bulge-dominated galaxies, consistent with a blue nugget phase;

• The Sérsic effective radius (r_e), disk effective radius (r_disk), and bulge effective radius (r_bulge) all show a decreasing trend with increasing redshift. Star-forming galaxies exhibit systematically larger sizes compared to quiescent galaxies at all redshifts, consistent with prior studies. Quiescent galaxies, while smaller than star-forming ones, show a steeper increase in their effective radius with time;

• Quiescent galaxies are significantly more compact than their star-forming counterparts, leading to high stellar mass surface densities (Σ_e). We find that Σ_e for quiescent galaxies is nearly an order of magnitude higher at z ~ 4 compared to z ~ 1, consistent with the observed compactness of some of the earliest quiescent galaxies observed by JWST.

This morphological catalog is a valuable addition to the DJA, enabling a range of in-depth studies of the morphological transformations associated with galaxy evolution.

Data availability

Our catalog is available at the CDS via anonymous ftp to cdsarc.cds.unistra.fr (130.79.128.5) or via https://cdsarc.cds.unistra.fr/viz-bin/cat/J/A+A/699/A343. It is also available from the DAWN JWST Archive. The code used to run the model fitting and the notebooks used to analyze the results are available on GitHub under the GNU GPL v3.S license.

Acknowledgements

The data products presented herein were retrieved from the DAWN JWST Archive (DJA). DJA is an initiative of the Cosmic Dawn Center (DAWN), which is funded by the Danish National Research Foundation under grant DNRF140. This work has received funding from the Swiss State Secretariat for Education, Research and Innovation (SERI) under contract number MB22.00072 This work has been conducted during the research internship of Aurélien Genin at the Cosmic DAWN Center, under the supervision of Marko Shuntov. This internship has also been made possible thanks to the financial support of the “Space: Science and Challenges of Space” chair at Ecole polytechnique, financed by ArianeGroup and Thales Alenia Space, as well as an Erasmus+ internship grant.

Appendix A PSF results

In Fig. A.1, we show the point spread functions (PSF) calculated over the GOODS-S field in different bands. As presented in Sect. 2.3, these empirical PSFs have been generated using PSFEX and a point-like sources selection based on MU_MAX and MAG_AUTO measurements.

Fig. A.1 PSFs calculated with PSFEx for different bands of JWST on the GOODS-S field. The point-like sources used are selected using the starline method. The name of the band (resp. channel) is written in the top right (resp. left) of each image. The images are displayed in logarithmic scale and color-inverted.

Appendix B SOURCEXTRACTOR++ model priors

In Fig. B.1, we show the priors used in our Sérsic and Bulge+Disk models (Sect. 3.2) for SOURCEXTRACTOR++. The first three panels correspond to the priors on the effective radius (for both models), Sérsic index and ellipticities for the Sérsic model, while the fourth panel shows the priors on the axis ratio that we adopt for the B+D fits.

Fig. B.1 Priors that we adopt for the source model parameters.

Appendix C Usage of Amazon Web Services

This work made use of Amazon Web Services (AWS), in particular the EC2² and S3³ services. AWS EC2 offers on-demand cloud computing, with the possibility to choose very different machines (refers as instances) with various resources in terms of computing power or memory. AWS S3 is a cloud storage service, already used to store the DJA. Here, we provide more details on the computing and use of AWS.

C.1 AWS EC2, VS Code, Jupyter interfacing

For simplicity of use, we developed some basic bash scripts to allow for easy interfacing between AWS EC2, VS code and Jupyter. These scripts enable to start an EC2 instance, launch a Jupyter server in it, and connect it to VS Code. The user can then use VS Code seamlessly to code, while running Jupyter notebooks directly on the EC2 instance. This has the added benefit of keeping the notebooks running, even if VS Code and the user computers are shutdown, which proves very useful for SOURCEXTRACTOR++ as it takes many hours to run. This code has been made available under the GNU GPL v3.0 license on GitHub⁴.

C.2 Automation

This section details the different steps used in this work to generate a morphological measurements catalog on one field. Step 1 is the pre-processing with PSF estimation and tiling. Step 2 is the result of the tiling process to run SOURCEXTRACTOR++ in parallel. Steps 3 and 4 are the merging of catalogs between all the tiles, the Sérsic and B+D models and the DJA. These steps are detailed in the following process:

1. Main instance (m5d.4xlarge for memory)

   • (a)

     Download all the frames of the field from the DJA;

   • (b)

     Pre-process them (decompress and save to a S3 bucket);

   • (c)

     Find point-like sources, estimate the PSFs in all the bands and save them to a S3 bucket;

   • (d)

     Tile the images and save them to a S3 bucket.

2. For every tile

   • (a)

     Start an instance (by default, c6a.4xlarge);

   • (b)

     Download the necessary files (frames, PSFs, configuration files, code);

   • (c)

     Run SOURCEXTRACTOR++ with the selected model (Sérsic or Bulge+Disk);

   • (d)

     Save the resulting catalog and images to a S3 bucket.

3. Main instance (m5d.4xlarge for memory)

   • (a)

     Merge tile catalogs and images;

   • (b)

     Save the full catalog and images to a S3 bucket.

4. Main instance (m5d.4xlarge for memory)

   • (a)

     Merge the full SOURCEXTRACTOR++ catalogs with the two models (Sérsic and Bulge+Disk) to the DJA catalogs;

   • (b)

     Save the complete catalog to the DJA.

Each of the sub-steps in this process are coded as bash and/or Python scripts. The dja_sepp Python package has also been developed to help with many of the different operations involved. It is available through pip and on GitHub⁵. This repository also holds many Jupyter notebooks used for the pre-and post-processing steps, as well as for the analysis of results (Sect. 4).

References

All Tables

All Figures

	Fig. 1 Point-like sources selection for PSF estimation. The plot shows the distribution (log scale) of sources from the GOODS-S field (F200W) in the MU_MAX/MAG_AUTO plane. The red box shows the starline selection of point-like sources, with the starline indicated by red dots.
In the text

	Fig. 2 Histogram of the flag values for the full catalogs. The percentages are relative to the total number of sources in the DJA catalogs. Two distinctive histograms for the Sérsic and Bulge+Disk models are shown because their fitting success rates differ. The figure shows a vast majority (70–80%) of sources are correctly fitted and provides science-ready morphological data.
In the text

Fig. 3 Comparison between Sérsic model fitting in our work and van der Wel et al. (2012) for UVJ-selected quiescent galaxies. The figure compares effective radii Reff and Sérsic index n. We corrected the Reff values from van der Wel et al. (2012), originally measured at 1.6 µm, to 2.5 µm using Eq. (1), in order to match the average wavelength of our JWST images. The color of the hexbin cells indicates the counts in each cell in log scale. Gray dotted lines show identical values between the two works. This figure includes 3263 matched galaxies.

In the text

Fig. 4 Distribution of Sérsic indices, nS, in the UVJ color space for different redshift ranges, for galaxies with log M★/M⊙ > 10. The color of each hexagonal bin represents the median value of nS, as indicated by the colorbar. The dotted line shows the quiescent vs. star-forming separation using Eq. (2). The quiescent region is predominantly populated by galaxies with high Sérsic indices (nS > 3).

In the text

Fig. 5 Distribution of B/T in the UVJ color space for different redshift and B/T ranges, for galaxies with log M★/M⊙ > 10. Each point is colored by its B/T value in the F200W band. The dotted line shows the quiescent versus star-forming separation using Eq. (2). The contours indicate 25, 50, and 75% of the density. Bulge-dominated galaxies (high B/T) predominantly occupy the quiescent region, whereas disk-dominated galaxies (low B/T) are generally in the star-forming region of the UVJ space.

In the text

	Fig. 6 UVJ color diagram for bulge- and disk-dominated log M★/M⊙ > 10 galaxies. The classification combines criteria based on both nS and B/T, estimated independently. Contour lines show kernel density estimates corresponding to 25, 50, and 75% of the density.
In the text

Fig. 7 Evolution of galaxy sizes as a function of redshift. The plots show the Sérsic effective radius re (left), the disk effective radius rdisk (center), and the bulge effective radius rbulge (right) as a function of zphot for log M★/M⊙ > 10 galaxies. This corresponds to a mean mass of ≈ 4 × 1010 M⊙. The points with error bars show the mean and 1 σ at a given redshift for all galaxies (black), disk-dominated galaxies (blue), and bulge-dominated galaxies (red). The re evolution is modeled in the left panel for all galaxies, as well as for the bulge- and disk-dominated populations, following re = R(1 + z)a (see Sect. 4.2). We compare the size evolution for all galaxies with previous JWST-based work from Varadaraj et al. (2024), Ormerod et al. (2024), Allen et al. (2025) and Yang et al. (2025).

In the text

Fig. 8 Evolution of the sizes of quiescent and star forming galaxies with log M★/M⊙ > 10. The plots show the Sérsic effective radius re, the disk effective radius rdisk, and the bulge effective radius rbulge as a function of zphot. Galaxies are classified as quiescent or star-forming based on their UVJ colors. The points and error bars indicate the mean and 1 σ dispersion in the corresponding zphot bin.

In the text

Fig. 9 Evolution of the stellar mass surface density, Σe. Shown are masses and redshifts from the DJA, calculated using EAZY, and effective radii, re, measured by SOURCEXTRACTOR++ in this work. Galaxies are classified as quiescent or star-forming according to their position in the UVJ diagram. The contour lines result from kernel density estimation using seaborn.kdeplot. The data for this figure includes 45 201 galaxies such that log M★/M⊙ > 9.

In the text

	Fig. A.1 PSFs calculated with PSFEx for different bands of JWST on the GOODS-S field. The point-like sources used are selected using the starline method. The name of the band (resp. channel) is written in the top right (resp. left) of each image. The images are displayed in logarithmic scale and color-inverted.
In the text

	Fig. B.1 Priors that we adopt for the source model parameters.
In the text