© The Authors 2023

1. Introduction

Galaxy morphology has been traditionally described in a qualitative way using the Hubble sequence (Hubble 1926), which is widely recognized as a fundamental aspect in the study of galaxy formation and evolution. With its unprecedented sensitivity and resolution in the infrared, the James Webb Space Telescope (JWST) is making significant advances in our understanding of the origin of the Hubble sequence. Previous studies using the Hubble Space Telescope (HST) suggested that the majority of galaxies at z > 2 are peculiar (e.g., Conselice et al. 2008; Mortlock et al. 2013). However, early JWST studies reveal a large fraction of regular disk galaxies at high redshift (Ferreira et al. 2022a,b; Kartaltepe et al. 2023; Nelson et al. 2023; Robertson et al. 2023; Jacobs et al. 2023; Cheng et al. 2023, 2022). These high-redshift disk galaxies can have spiral arms and bars similar to those in Local Universe (Wu et al. 2023; Chen et al. 2022; Fudamoto et al. 2022; Guo et al. 2023). Galaxies with established disk and spheroidal morphologies span the full redshift range (Kartaltepe et al. 2023) and the Hubble Sequence was already in place as early as z ≈ 6 (Ferreira et al. 2022b).

Quantifying morphology is crucial for exploring galaxy evolution and it can be achieved using both non-parametric and parametric methods (e.g., Elmegreen & Elmegreen 1985; Conselice 2003; Lotz et al. 2004; Elmegreen et al. 2007; Martínez-García et al. 2014; Yu & Ho 2018, 2019, 2020; Yu et al. 2021, 2022a,b; Smith et al. 2022). The Cosmic Evolution Early Release Science (CEERS) survey (PI: Finkelstein, ID = 1345, Finkelstein et al. 2022; Bagley et al. 2023) is an early release science program in Cycle 1 that observes the Extended Growth Strip field (EGS) of the Cosmic Assembly Near-IR Deep Extragalactic Legacy Survey (CANDELS; Grogin et al. 2011; Koekemoer et al. 2011). The CEERS will observe ten pointings with the Near-Infrared Camera (NIRCam; Rieke et al. 2023), covering a total of 100 sq. arcmin. By quantifying morphology of galaxies from the first four pointings, Ferreira et al. (2022b) and Kartaltepe et al. (2023) show that spheroids exhibit a higher average Sérsic index, smaller size, and rounder shape compared to disks and peculiars, although selection effects may exist. Additionally, the average Sérsic index decreases with increasing redshift (Ferreira et al. 2022b). Despite the slightly higher average concentration in spheroids and slightly higher average asymmetry in peculiars, the concentration-asymmetry diagram does not provide a clear separation of galaxies according to their morphological types (Ferreira et al. 2022b; Kartaltepe et al. 2023).

However, these early results may not accurately reflect the intrinsic galaxy morphology as the physical resolution and signal-to-noise ratio (S/N) of the JWST images of high-redshift galaxies are limited. Furthermore, comparing galaxy morphologies at different redshifts and/or observed by different instruments is challenging as changes in resolution, noise level, and rest-frame wavelength can alter the morphology and introduce biases and uncertainties in the quantification.

A commonly used strategy for understanding measurement biases and uncertainties caused by image degradation is to use high-quality images of low-redshift galaxies to generate simulated images of high-redshift galaxies and, subsequently, to compare the measurements before and after the image simulation. The pioneering work in this area was done by Giavalisco et al. (1996). This strategy has been used to study spiral structure (e.g., Block et al. 2001; van den Bergh et al. 2002; Yu et al. 2018), bar structure (e.g., van den Bergh et al. 2002; Sheth et al. 2008), concentration-asymmetry-smoothness statistic (e.g., Abraham et al. 1996; Conselice 2003; Lee et al. 2013; Yeom et al. 2017; Whitney et al. 2021), Gini-M₂₀ statistic (e.g., Lee et al. 2013; Petty et al. 2014), and Sérsic index and size (e.g., Barden et al. 2008; Paulino-Afonso et al. 2017) of high-redshift galaxies observed by HST. Despite the widespread use of this technique, their image simulation procedures did not take into account the intrinsic galaxy size evolution and cannot be directly applied to galaxies with high redshifts. It has been found that there is a strong redshift evolution in galaxy size (e.g., Bouwens et al. 2004; Daddi et al. 2005; Trujillo et al. 2007; Buitrago et al. 2008; Oesch et al. 2010; Mosleh et al. 2012; van der Wel et al. 2014; Whitney et al. 2019). By using data from CANDELS, van der Wel et al. (2014) show that galaxies of a given stellar mass are on average smaller at higher redshifts, with fast evolution for early-type galaxies and moderate evolution for late-type galaxies.

Building on prior studies and by taking into account all relevant factors of image simulation, our goal is to use high-resolution and high-S/N nearby galaxy images to generate artificial images of galaxies located at redshift 0.75 ≤ z ≤ 3 and observed at optical rest-frame wavelength in JWST CEERS. We then go on to investigate biases and uncertainties in the quantification of galaxy morphology, focusing on six key morphological quantities: Petrosian radius, half-light radius, asymmetry, concentration, Sérsic index, and axis ratio, which are commonly used to describe the typical morphology of a galaxy. We aim to derive corrections to these biases and uncertainties to improve the accuracy and robustness of future galaxy morphology studies.

This paper is organized as follows. Section 2 outlines our sample selection of nearby galaxies and the data reduction process. Section 3 provides a detailed description of the image simulation methodology. Section 4 discusses the biases and uncertainties, and derive correction functions. Section 5 validates the effectiveness of the correction functions. Finally, a summary of the main findings is presented in Sect. 6. Throughout this work, we use AB magnitudes and assume the following cosmological parameters: (Ω_M, Ω_Λ, h) = (0.27, 0.73, 0.70).

2. Observational material

2.1. Sample of nearby galaxies

We restricted the redshift range of our image simulations to z = 0.75, 1.0, 1.25, …, and 3.0. We refrain from simulating images of galaxies at z > 3, as the of evolution galaxy optical size at these redshifts is not yet well constrained. For each target redshift, we chose the filter that observes the rest-frame optical wavelength (λ = 5000 − 7000 Å). Specifically, we used the F115W filter for z = 0.75 and 1, the F150W filter for z = 1.25, 1.5, and 1.75, and the F200W filter for z = 2.0 to 3.0. To generate artificially redshifted galaxy images, nearby (z ≈ 0) galaxy images observed in a similar rest-frame wavelength range are required. This range is covered by the g, r, and z bands (with effective wavelengths: 4796 Å, 6382 Å, and 9108 Å) provided by the Dark Energy Spectroscopic Instrument (DESI) Legacy Imaging Surveys¹ (Dey et al. 2019, and references therein).

The DESI Legacy Imaging Surveys are comprised of three public projects: the Beijing-Arizona Sky Survey (BASS; Zou et al. 2017), the Mayall z-band Legacy Survey (MzLS; Blum et al. 2016), and the Dark Energy Camera Legacy Survey (DECaLS; Blum et al. 2016). We focus on the DECaLS g-, r, and z-band images, while the MzLS z-band nearby galaxy images are found to be significantly affected by pattern noise and are therefore excluded from our analysis along with their corresponding BASS g- and r-band images. We defined our sample of nearby galaxies using the Siena Galaxy Atlas (SGA)², which is constructed based on the DESI and includes 383 620 galaxies. The SGA does not include small objects with D₂₅ < 20 arcsec, where D₂₅ is the B-band isophotal diameter at μ_B = 25 mag arcsec⁻², as a high fraction of them are spurious sources.

We selected galaxies with available Hubble type classification and best-estimated luminosity distance (D_L) from HyperLeda extragalactic database³ (Makarov et al. 2014, and references therein). The best-estimated distance is weighted average between distance derived from spectroscopic redshift and published redshift-independent distances, and it provide an homogeneous distance estimate over the whole redshift range (Makarov et al. 2014). We select galaxies with 12.88 ≤ D_L ≤ 65.01 Mpc, corresponding to cosmological z = 0.003 − 0.015. We exclude galaxies in the Galactic plane (−20° ≤b ≤ 20°) to avoid any possible severe photometric problems caused by crowd foreground stars. We excluded galaxies that are covered by masks of nearby bright sources provided by the SGA, as the emission from these galaxies is severely contaminated by the bright sources. This process will remove merging systems where the centers of two galaxies are close but have not yet merged into one. The SGA catalog also includes mid-infrared photometry from the Wide-field Infrared Survey Explorer (WISE; Wright et al. 2010). We use the DECaLS g-, r-, and z-band flux and WISE W1 flux to estimate stellar mass (M_⋆) through SED fitting using CIGALE⁴ (Boquien et al. 2019), assuming the Chabrier stellar initial mass function (Chabrier 2003), double exponential star formation history, simple stellar population of Bruzual & Charlot (2003), and attenuation law of Calzetti et al. (2000). We selected galaxies with M_⋆ of 10^{9.75 − 11.25} M_⊙. The mass cut was applied because we use the galaxy size evolution derived by van der Wel et al. (2014), which is only available in this mass range for both early- and late-type galaxies. Our final sample of nearby galaxies consists of 1816 galaxies.

Figure 1 summarizes some of the basic parameters of the sample. Most of the galaxies are nearby (median D_L = 52.2 Mpc; Fig. 1a), luminous (median M_r = −20.4 mag, corrected for Galactic extinction; Fig. 1b), massive (median M_⋆ = 10^10.2 M_⊙; Fig. 1c), and angularly large (median D₂₅ = 1.4 arcmin; Fig. 1d). The sample spans the full range of Hubble types in the nearby Universe (Fig. 1e), comprising 600 (33%) early-type galaxies (ellipticals, S0, or S0/a), and 1216 (67%) late-type galaxies (Sa–Irr). We acquired g-, r-, and z-band mosaic images and point spread functions (PSFs) from SGA. The Galactic extinction is corrected using the map of dust reddening (Schlegel et al. 1998; Schlafly & Finkbeiner 2011). We use the galaxy center, ellipticity, and position angle provided by the SGA catalog when removing foreground stars (as described in Sect. 2.2). We set R₂₅ = D₂₅/2. We measured the sky background and noise using AutoProf⁵ (Stone et al. 2021). The sky background is then subtracted from the image.

Fig. 1. Basic properties of our sample. Distribution of (a) luminosity distances, (b) absolute r-band magnitude, corrected for Galactic extinction, (c) stellar mass, (d) B-band isophotal diameter at μB = 25 mag arcsec−2, and (e) Hubble types.

2.2. Removal of foreground stars

Contamination from sources other than the target galaxy, such as foreground stars and projected close galaxies, should be removed or minimized prior to image simulation. This process was not done in some previous studies on simulating high-redshift galaxy images observed by HST (e.g., Barden et al. 2008; Paulino-Afonso et al. 2017, but also see Yu et al. 2018), rendering their results less robust. We first mask out the contamination. For each galaxy image, the SGA provides a catalog of sources that are identified using The Tractor⁶, a forward-modeling approach to performing source extraction and model fitting to the sources. We masked out each identified star, using a mask size determined by the star-centered radius at which the r-band flux start to fluctuate due to galaxy structure or background noise. We then masked out each identified projected nearby galaxy using an ellipse with semi-major axis of 1.5 R₂₅. Finally, we performed a visual inspection and manually masked out any residual stars or small background galaxies that were missed in the above process. We adopted the same mask for g-, r-, and z-band images.

We removed contamination through the following steps. The masked regions outside the galaxy (R > 1.5 R₂₅) are set to zero. For small masked regions (area < 50 arcsec²) inside the galaxy (R ≤ 1.5 R₂₅), we estimated the intrinsic galaxy light affected by foreground stars and/or projected close galaxies using interpolation. For each masked region, we cut out a square region containing twice as many unmasked pixels as masked pixels. The interpolation was then performed by approximating the values of the unmasked pixels by a polynomial function. We used the interpolated values to fill in the small masked regions.

For large masked regions (area ≥50 arcsec²) inside the galaxy, interpolation may fail to reproduce the intrinsic galaxy light, as the complex galaxy structures may cause overfitting and lead to catastrophic results. Instead, we replaced the masked region with values from their 180° rotational symmetric pixels if the symmetric portion did not have a large masked region. The rotational symmetric images were originally used to highlight spiral structure (Elmegreen et al. 1992, 2011). We caution that this approach may not restore the flux in prominent three-armed structures, which are 120° rotational symmetric, but this effect is small since the fraction of three-armed structures is small (Elmegreen & Elmegreen 1987). In cases where the 180°-rotational symmetric portion also has a large masked region or the galaxy is highly inclined (i > 70°), we filled in the masked region with values from their mirror-symmetric pixels reflected over the galaxy major axis. In a few instances where neither of the above criteria are satisfied, we reverted to the interpolation method and used a low-order polynomial function to perform fitting and interpolation. Finally, we added Poisson noise to the cleaned regions to simulate real observations.

The above process makes use of galaxy symmetry, but it does not significantly affect the calculation of galaxy asymmetry described in Sect. 4, as the masked regions are small relative to the galaxy size. We have to sacrificed a small degree of precision in computing the galaxy asymmetry to generate cleaned images for image simulation.

The removal of foreground stars is done for each galaxy to generate their star-cleaned g-, r-, and z-band images. The effectiveness of the cleaning is illustrated in Fig. 2, showing the r-band images of two galaxies, ESO 121−026 (a Sbc galaxy) and ESO 251−004 (an elliptical galaxy), before and after undergoing the removal process. Our process successfully eliminates the vast majority of contamination surrounding the galaxies, revealing clearer and more accurate images of the galaxy structures. The median PSF FWHM of DECaLS images is ∼1.3, 1.2, and 1.1 arcsec in the g, r, and z bands, respectively (Dey et al. 2019). To facilitate the pixel-by-pixel K correction described in Sect. 3.4, we matched the g-, r-, and z-band images to a common PSF for each galaxy. We did not use Fourier transformation to find a convolution kernel, because this would introduce significant high frequency noise, as the PSF FWHMs at different bands are very close. Instead, the matching was done by searching for a kernel of Moffat function, which is convolved with the PSF of smaller FWHM to get a broadened PSF that has almost the same best-fitted Moffat function with the PSF of larger FWHM. The star-cleaned g-, r-, and z-band images with a common PSF were used to compute artificial high-redshift galaxy images observed in JWST CEERS.

Fig. 2. Demonstration of our cleaning process. Left and right columns showcase the r-band observation of two galaxies, ESO 121−026 (a Sbc galaxy) and ESO 251−004 (an elliptical galaxy). Top and bottom rows display the image before and after undergoing star removal, respectively. The vast majority of emission from foreground stars or background galaxies in the vicinity of the galaxies was successfully removed.

3. Artificially redshifting galaxies

The nearby DESI galaxy images are of fairly high quality, which allows for accurate determinations of morphological measurements. Compared to DESI images, JWST CEERS images of high-redshift galaxies have lower physical resolution and lower S/N, as the galaxies become less resolved and fainter at higher redshifts. The limited data quality may bias measurements and make them more uncertain. To understand the how well we can quantify galaxy morphology using the JWST NIRCam images, we simulated nearby DESI galaxies with respect to how they would appear at various high redshifts, as observed in JWST CEERS. Section 3.1 presents the derivation of the formulae used in the redshifting procedure. Subsequently, Sects. 3.2 and 3.3 respectively discuss the evolution of galaxy size and luminosity. The redshifting procedure is outlined in Sect. 3.4.

3.1. Formulae

We begin by assuming the existence of an extended source located at redshift, z_local. Its luminosity distance and flux density we observed are denoted as D_L(z_local) and f_local, respectively. We adopted erg s⁻¹ cm⁻² Hz⁻¹ as the unit of f because we use AB magnitude. The energy emitted by the source in 1 s at a narrow range of wavelength from λ − Δλ/2 to λ + Δλ/2 is given by:

$$\begin{array}{c}\hfill E_{\mathrm{local}}=4\pi ·\frac{c}{{\lambda }^2}·f_{\mathrm{local}}·D_{\mathrm{L}}^2(z_{\mathrm{local}})·\mathrm{\Delta }\lambda ,\end{array}$$(1)

where c is the light speed. The source is manually moved to a higher redshift of z_high with the luminosity distance of D_L(z_high). The observed flux density is denoted as f_high. The observed wavelength becomes λ′=λ ⋅ (1 + z_high)/(1 + z_local) and the wavelength width becomes Δλ′=Δλ ⋅ (1 + z_high)/(1 + z_local). The energy emitted by the source in 1 second at a narrow range of wavelength from λ′−Δλ′/2 to λ′+Δλ′/2 is:

$$\begin{array}{c}\hfill E_{\mathrm{high}}=4\pi ·\frac{c}{{\lambda }^{\prime 2}}·f_{\mathrm{high}}·D_{\mathrm{L}}^2(z_{\mathrm{high}})·\mathrm{\Delta }{\lambda }^{\prime }.\end{array}$$(2)

We assume that the source becomes intrinsically brighter at z_high, with the energy emitted given by:

$$\begin{array}{c}\hfill E_{\mathrm{high}}={\left(\frac{1+z_{\mathrm{high}}}{1+z_{\mathrm{local}}}\right)}^{\alpha }·E_{\mathrm{local}}.\end{array}$$(3)

Combining Eqs. (1)–(3), we obtain:

$$\begin{array}{c}\hfill f_{\mathrm{high}}=\frac{D_L^2(z_{\mathrm{local}})}{D_L^2(z_{\mathrm{high}})}·\frac{1+z_{\mathrm{high}}}{1+z_{\mathrm{local}}}·{\left(\frac{1+z_{\mathrm{high}}}{1+z_{\mathrm{local}}}\right)}^{\alpha }·f_{\mathrm{local}}.\end{array}$$(4)

The first term on the right-hand side of the equation is caused by the distance effect and the second term is caused by the cosmological compression of the frequency (or, equivalently, the cosmological dilation of the wavelength). The second term occurs as we consider monochromatic luminosity, while it would disappear if we considered bolometric luminosity. The third term is caused by the assumed luminosity evolution. Equation (4) is used to rescale the flux in the redshifting procedure (Sect. 3.4). Hence, we can relate the absolute magnitude of the source at local and high redshift using:

$$\begin{array}{c}\hfill M_{\mathrm{high}}=-2.5log{\left(\frac{1+z_{\mathrm{high}}}{1+z_{\mathrm{local}}}\right)}^{\alpha }+M_{\mathrm{local}}.\end{array}$$(5)

We assume the source has a physical radius of R_local at z_local and R_high at z_high. The solid angle it spans at z_local and z_high is:

$$\begin{array}{c}\hfill {\mathrm{\Omega }}_{\mathrm{local}}=\frac{\pi R_{\mathrm{local}}^2}{D_{\mathrm{ang}}^2(z_{\mathrm{local}})},\end{array}$$(6)

and

$$\begin{array}{c}\hfill {\mathrm{\Omega }}_{\mathrm{high}}=\frac{\pi R_{\mathrm{high}}^2}{D_{\mathrm{ang}}^2(z_{\mathrm{high}})},\end{array}$$(7)

respectively. D_ang(z) is the angular-diameter distance at redshift z. We assume the physical size becomes smaller at higher redshift:

$$\begin{array}{c}\hfill R_{\mathrm{high}}={\left(\frac{1+z_{\mathrm{high}}}{1+z_{\mathrm{local}}}\right)}^{\beta }·R_{\mathrm{local}}.\end{array}$$(8)

Combining Eqs. (6)–(8), we obtain:

$$\begin{array}{c}\hfill {\mathrm{\Omega }}_{\mathrm{high}}=\frac{D_{\mathrm{ang}}^2(z_{\mathrm{local}})}{D_{\mathrm{ang}}^2(z_{\mathrm{high}})}·{\left(\frac{1+z_{\mathrm{high}}}{1+z_{\mathrm{local}}}\right)}^{2\beta }·{\mathrm{\Omega }}_{\mathrm{local}}.\end{array}$$(9)

We denote the pixel scale as p and number of pixels occupied by the source as N. We therefore have:

$$\begin{array}{c}\hfill \frac{N_{\mathrm{high}}}{N_{\mathrm{local}}}=\frac{D_{\mathrm{ang}}^2(z_{\mathrm{local}})}{D_{\mathrm{ang}}^2(z_{\mathrm{high}})}·{\left(\frac{1+z_{\mathrm{high}}}{1+z_{\mathrm{local}}}\right)}^{2\beta }·\frac{p_{\mathrm{local}}^2}{p_{\mathrm{high}}^2},\end{array}$$(10)

which is the binning factor we use to downscale image size in the redshifting procedure (Sect. 3.4). The flux density per unit solid angle can be written as:

$$\begin{array}{c}\hfill \frac{f_{\mathrm{high}}}{{\mathrm{\Omega }}_{\mathrm{high}}}={\left(\frac{1+z_{\mathrm{high}}}{1+z_{\mathrm{local}}}\right)}^{-3}·{\left(\frac{1+z_{\mathrm{high}}}{1+z_{\mathrm{local}}}\right)}^{\alpha }·{\left(\frac{1+z_{\mathrm{high}}}{1+z_{\mathrm{local}}}\right)}^{-2\beta }·\frac{f_{\mathrm{local}}}{{\mathrm{\Omega }}_{\mathrm{local}}}·\end{array}$$(11)

Denoting the surface brightness as μ, we have:

$$\begin{array}{cc}\hfill {\mu }_{\mathrm{high}}=& {\mu }_{\mathrm{local}}+2.5log{\left(\frac{1+z_{\mathrm{high}}}{1+z_{\mathrm{local}}}\right)}^3-2.5log{\left(\frac{1+z_{\mathrm{high}}}{1+z_{\mathrm{local}}}\right)}^{\alpha }\hfill \\ \hfill & -2.5log{\left(\frac{1+z_{\mathrm{high}}}{1+z_{\mathrm{local}}}\right)}^{-2\beta }.\hfill \end{array}$$(12)

If we adopt z_local = 0 and z_high = z, we have:

$$\begin{array}{c}\hfill {\mu }_z={\mu }_0+2.5log{(1+z)}^3-2.5log{(1+z)}^{\alpha -2\beta }.\end{array}$$(13)

The term of 2.5log(1 + z)³ is the well-known cosmological dimming, while the term of 2.5log(1 + z)^{α − 2β} describes the evolution of surface brightness. In the redshifting procedure, we assumed that the average flux density observed in a specific filter is consistent with the flux density at the effective wavelength of the filter, so that we can ignore the difference between bandpass width of the DESI filter and JWST filter.

3.2. Size evolution

Galaxies at a fixed mass are physically smaller at higher redshift (e.g., Bouwens et al. 2004; Daddi et al. 2005; Trujillo et al. 2007; Buitrago et al. 2008; Oesch et al. 2010; Mosleh et al. 2012; van der Wel et al. 2014; Whitney et al. 2019). Using structural parameters derived from CANDELS imaging (Grogin et al. 2011; Koekemoer et al. 2011), van der Wel et al. (2014) found a significantly different rate (β) of average size evolution for early-type and late-type galaxies. The average effective radius of early-type galaxies, calculated over a range of stellar masses, evolves rapidly, following R_eff ∝ (1 + z)^{β = −1.48}, while that of late-type galaxies evolves moderately, following R_eff ∝ (1 + z)^{β = −0.75}. Therefore, when simulating high-redshift galaxies, it is important to properly account for intrinsic size evolution; otherwise, the simulated galaxies would be larger than the true galaxies of the same mass.

The evolution rate is dependent on stellar mass, with more massive galaxies showing a higher rate (more negative β), as reported in van der Wel et al. (2014). The dependence of β on mass is weak for late-type galaxies, but strong for early-type galaxies, and less massive early-type galaxies evolve at a similar rate to late-type galaxies of the same mass. We adopted the β measured by van der Wel et al. (2014, as given in their Table 2) We focused on galaxies with a stellar mass of 10^{9.75 − 11.25} M_⊙, so that β is available for both late-type and early-type galaxies. To estimate β for our nearby DESI galaxies, we fit two polynomials to the β as a function of stellar mass for late-type and early-type galaxies, respectively. These best-fit functions were used to compute β using stellar mass and galaxy type as input. The estimated β was used for generating artificially redshifted galaxy images (see Sect. 3.4). Additionally, galaxies tend to be larger at bluer wavelength than at redder wavelength (e.g., Kelvin et al. 2012). We performed a pixel-by-pixel K correction in Sect. 3.4 to correct for this effect.

The rate of size evolution may differ when using different definitions of galaxy size. As shown by Whitney et al. (2019), the outer radius of galaxies evolves at a faster rate than the inner radius, suggesting inside-out growth. While we do not consider the effect of inside-out growth in this study, we aim to compare our simulated images with real observations to investigate this effect in the future.

3.3. Monochromatic luminosity evolution

The intrinsic galaxy surface brightness has been found to brighten on average with increasing redshift (e.g., Schade et al. 1995, 1996; Lilly et al. 1998; Roche et al. 1998; Labbé et al. 2003; Barden et al. 2005; Sobral et al. 2013; Whitney et al. 2020). By definition, the evolution of surface brightness is partly attributed to intrinsic size evolution and partly attributed to monochromatic luminosity evolution. While the intrinsic size evolution is discussed in Sect. 3.2, in this section we study the monochromatic luminosity evolution by assuming a function form of L_λ ∝ (1 + z)^α, where L_λ is the monochromatic luminosity at rest-frame λ filter with the effect of cosmological dimming corrected.

To stay consistent with van der Wel et al. (2014), we used the rest-frame flux, magnitude, color index, and redshift from the 3D-HST catalog (Brammer et al. 2012; Skelton et al. 2014). Following the strategy in van der Wel et al. (2014), we selected CANDELS galaxies with F160W apparent magnitude brighter than 25.5, a flag for good model fittings, and stellar mass above the mass of completeness limit at each redshift range. We classified the selected CANDELS galaxies into early-type and late-type galaxies using the demarcation lines proposed by Williams et al. (2009) in the diagram of U − V versus V − J color index. We separated the early-type and late-type galaxies to various bins of redshift (Δz = 0.5) and stellar mass (Δlog M_* = 0.5 dex), as used in van der Wel et al. (2014). Then we calculated the median rest-frame absolute magnitude M_λ, where λ denotes U, V, B, R, I, or J filter. We fit the following function:

$$\begin{array}{c}\hfill M_{\lambda }=M_{\lambda ,0}-2.5log{(1+z)}^{\alpha },\end{array}$$(14)

to the M_λ as a function of z for each bin of stellar mass viewed at each filter to determine α. The results are shown in Figs. 3 and 4. We then fit two 2D polynomials to the measured α as a function of filter effective wavelength and stellar mass for early-type and late-type galaxies, respectively. Figures 5 and 6 show the best-fit 2D functions. The color bar encodes best-fit α for a given wavelength and mass. Both early-type and late-type galaxies get brighter on average at higher redshifts across all wavebands. The rate of evolution is more rapid in bluer wavebands than in redder ones, owing to more intense star formation in the past (Speagle et al. 2014; Scoville et al. 2023). Interestingly, the low-mass late-type galaxies have particularly high α at blue wavelength. Our results are consistent with the evolution of the characteristic magnitude of luminosity function, where the characteristic magnitude becomes brighter at higher redshifts and evolves faster in the UV than in the V band (e.g., Arnouts et al. 2005; Marchesini et al. 2012). We used the rest-frame wavelength, galaxy type, and stellar mass as inputs for estimate α for our nearby DESI galaxies by using the best-fit polynomials.

Fig. 3. Redshift evolution of intrinsic absolute magnitude at rest-frame waveband of U, B, V, R, I, and J for early-type galaxies. The measured parameter α (Eq. (14)) describing the evolution rate for a given mass range and a given rest-frame waveband is shown in the bottom-right corner.

Fig. 4. Same details as in Fig. 3, but for late-type galaxies.

Fig. 5. Best-fit α as a function of filter effective wavelength (λeff) and stellar mass (M*) for early-type galaxies. The color bar encodes α.

Fig. 6. Same details as in Fig. 5, but for late-type galaxies.

3.4. Redshifting procedure

Our redshifting procedure is based on the method of Giavalisco et al. (1996), which we updated by incorporating the K correction, galaxy size evolution, and luminosity evolution. This approach ensures that the artificially redshifted galaxy images closely match the size and brightness of true galaxies viewed at high redshifts. We compute artificial images of galaxies at z = 0.75, 1.0, 1.25, 1.5, 1.75, 2.0, 2.25, 2.5, 2.75, and 3.0, using the F115W filter for z = 0.75 and 1, the F150W filter for z = 1.25, 1.5, and 1.75, as well as the F200W filter for z = 2.0 to 3.0. To determine the typical noise level for the simulated images, we use the science data, error map, and source mask from the CEERS Data Release Version 0.5 provided by Bagley et al. (2023). In order to simulate background noise at each filter, we cut out a patch with a few sources from the CEERS science data, and then replaced these sources with background regions next to the sources to get a clean fragment of a real background image. We selected 80 galaxies, consisting of the top 20 brightest galaxies from each of the four pointings, to calculate a median ratio of galaxy flux variance to galaxy flux, used to estimate the noise from galaxy flux. This quantity depends on filter, exposure time, sky brightness, and system throughput, and varies only slightly from location to location in CEERS. We adopted a pixel scale of 0.03 arcsec/pixel, same as in Bagley et al. (2023). We generated a two-time oversampling PSFs at the F115W, F150W, and F200W filters using WebbPSF (Perrin et al. 2014). The two-time oversampling PSFs were used because the F115W and F150W PSFs on the pixel scale of 0.03 arcsec/pixel were undersampled.

Our Python algorithm for generating artificially redshifted galaxy images is summarized in the following steps: (1) We downscaled the size of g-, r-, and z-band DESI images so that the PSF occupy two pixels while preserving their total flux. This is meant to reduce the computation time of step 2, while retaining Nyquist sampling; (2) For each pixel with S/N ≥ 3, we performed a K correction using the python code developed by Blanton & Roweis (2007)⁷ to calculate the expected flux at the rest-frame wavelength through interpolation. Later, a median ratio of the rest-frame flux derived above to the DESI image flux was calculated and multiplied with the value of pixels with S/N < 3 to obtain a K-corrected DESI image; (3) We rescaled the flux using three factors (Eq. (4)). The first one is a dimming factor, caused by the longer luminosity distance at higher redshift. The second one is a brightening factor, caused by the cosmological dilation of wavelength. The first two factors lead to the well-known cosmological dimming. The third factor comes from the monochromatic luminosity evolution, such that the monochromatic luminosity is scaled with (1 + z)^α; (4) We downscaled the image size to match the half pixel size and the galaxy angular size as if it appears at high redshifts, while preserving the total flux. In addition to the change of angular size due to longer angular-diameter distances, we took into account the intrinsic galaxy size evolution that is scaled with the physical size (1 + z)^β; (5) We used photutils (Bradley et al. 2022) to calculate, by Fourier transformation, a kernel that transforms input PSF to the two-time oversampling PSF. Since the input PSF is much smaller, no high frequency noise occurs. The image is convolved with the kernel and was then downscaled in size by 50% to match the target pixel scale of 0.03 arcsec/pixel to obtain a resolution-matched image; (6) We calculated the variance map by multiplying the galaxy flux with the median ratio of variance to flux, generated noise using the map, and added the resulting flux noise to the image. Next, we overlaid the resolution-matched image on top of the clean real background image to produce the final simulated CEERS image.

We performed the six steps above for each galaxy to generate its simulated CEERS images at high redshifts. As outlined in step 2, we obtained rest-frame images by interpolating multi-band data at the pixel level, avoiding extrapolation to minimize the risk of introducing significant errors in output flux. To evaluate the uncertainty of this process, we performed a test by using 100 time bootstrap resamplings for one galaxy, which resulted in 100 more simulated images. We found that the resulting uncertainty in the output flux is small, considerably smaller than the typical background noise and galaxy flux noise. Therefore, this source of uncertainty has negligible impact on our results (see also Barden et al. 2008). Figure 7 illustrates artificially redshifted galaxy images of three example galaxies. The first column shows the DESI r-band image, while the second to forth columns illustrate the artificial JWST CEERS images at z = 1, z = 2, and z = 3. The final column plots real images of three CEERS galaxies, each with stellar masses comparable to the galaxies in the same row. Although some small-scale structures are suppressed by noise or PSF smoothing effects, the galactic-scale structures such as strong bars and prominent spirals are still present. Overall, the galaxies remain visible across all redshifts, although they become noisy and blurry. Compared to the simulated CEERS images, real CEERS images of disk looks clumpier. A detailed and rigorous comparison between them is crucial to study the evolution of structure in the future.

Fig. 7. Examples of artificially redshifted galaxy images and the comparison between simulated and real CEERS images. From left to right: first four panels in each row show the DESI r-band image and the artificial JWST CEERS images at z = 1, z = 2, and z = 3. The fifth column presents real galaxy images in JWST CEERS for comparison. Apparent scales are indicated in each panel.

4. Correction of biases and uncertainties

In this section, we quantify the measurement biases and uncertainties for six commonly used morphological quantities: Petrosian radius, half-light radius, concentration, asymmetry, Sérsic index, and axis ratio. The high-quality, K-corrected DESI images enable us to accurately determine the intrinsic morphological measurements. We used the z = 2.0 high-quality K-corrected DESI images as the training set to derive functions for correcting the biases and uncertainties arising from the resolution effects present in the measurements obtained from degraded images. Using K-corrected images at other redshifts yields nearly identical results, as the only difference lies in the slight variation in observed rest-frame wavelengths. We reduced their image resolution so that the intrinsic Petrosian radius is N times the PSF FWHM. We then performed measurements and compared the resulting values with their intrinsic counterparts. We adopted exponentially-growing values of N, which are 1.98, 3, 4.55, 6.89, 10.45, 15.83, and 24. The FWHM values for the F115W, F150W, and F200W PSFs are 0.037, 0.049, and 0.064 arcsec, respectively. These images are denoted as N-FWHM images. The typical CEERS noise predominantly impacts the computation of asymmetry. We used the z = 3 simulated CEERS images, which are the noisiest in our dataset of simulated galaxies, to improve the method for removing noise contribution from the computation of asymmetry.

4.1. Galaxy size

We estimated the flux-weighted center and apparent projection parameters for each image, and measure the Petrosian radius (R_p; Petrosian 1976), defined as the radius at which the surface brightness is 20% of the average surface brightness within R_p. We re-measured the center by minimizing galaxy asymmetry (see Sect. 4.3) and used it to re-measure R_p as our final measurement; R_p encompasses at least 99% of the light within a given galaxy (Bershady et al. 2000).

The half-light radius is another indicator of galaxy size, used to study galaxy evolution (e.g., Bouwens et al. 2004; Buitrago et al. 2008). By measuring the total flux within an elliptical aperture with a radius of 1.5 R_p, we derived the fraction of light that the radius enclose as a function of radius, known as the curve of growth (cog). We then determine R₂₀, $R_{50}^{\mathrm{cog}}$, and R₈₀, the radius containing 20%, 50%, and 80% of the total galaxy light, respectively. In additional to the non-parametric approach, we fit a Sérsic model to the galaxy to obtain the half-light radius, denoted as $R_{50}^{\mathrm{fit}}$ (see Sect. 4.4).

We defined the level of image resolution as R_p, True/FWHM. We measured R_p, $R_{50}^{\mathrm{cog}}$, and $R_{50}^{\mathrm{fit}}$ on the F200WN-FWHM images and plot the differences between them and their intrinsic values (R_p, True, $R_{50,\mathrm{True}}^{\mathrm{cog}}$, and $R_{50,\mathrm{True}}^{\mathrm{fit}}$) as a function of R_p, True/FWHM in the first row of Fig. 8. The mean value of (R_p − R_p, True)/FWHM, denoted by Δ_p, is 0.770, and that of $(R_{50}^{\mathrm{cog}}-R_{50,\mathrm{True}}^{\mathrm{cog}})/FWHM$, denoted by ${\mathrm{\Delta }}_{50}^{\mathrm{cog}}$, is 0.352. Thus, R_p and $R_{50}^{\mathrm{cog}}$ are systematically slightly overestimated due to the lower resolution, and the biases should be corrected. Interestingly, the mean difference does not significantly depends on R_p, True/FWHM. The detected bias in measuring R_p is consistent with Whitney et al. (2019), who show that the measured R_p before and after image blurring correlate well with slop of ∼1 and with a systematic offset. In contrast, the mean value of $(R_{50}^{\mathrm{fit}}-R_{50,\mathrm{True}}^{\mathrm{fit}})/FWHM$, denoted by ${\mathrm{\Delta }}_{50}^{\mathrm{fit}}$, is so small: −0.084, indicating that the Sérsic fitting can extract the half-light radius without statistical bias, and no correction is needed.

Fig. 8. Evaluation of biases and uncertainties in measuring Petrosian radius (Rp; left), half-light radius using a non-parametric method ($R_{50}^{\mathrm{cog}}$; middle), and half-light radius through Sérsic model fitting ($R_{50}^{\mathrm{fit}}$; right). The top row illustrates the difference between the measured and intrinsic values as a function of Rp, True/FWHM, with the solid horizon line marking the mean difference. The middle row displays the correlation between the corrected and intrinsic values, with the one-to-one relation indicated by a dashed line. The bottom row presents the fractional uncertainties, together with their best-fit functions plotted as solid lines.

The measurements of R_p and $R_{50}^{\mathrm{cog}}$ are corrected for the biases using the functions:

$$\begin{array}{c}\hfill R_{\mathrm{p},\mathrm{cor}}/\mathit{FWHM}=R_{\mathrm{p}}/\mathit{FWHM}-{\mathrm{\Delta }}_{\mathrm{p}},\end{array}$$(15)

and

$$\begin{array}{c}\hfill R_{50,\mathrm{cor}}^{\mathrm{cog}}/\mathit{FWHM}=R_{50}^{\mathrm{cog}}/\mathit{FWHM}-{\mathrm{\Delta }}_{50},\end{array}$$(16)

respectively. In Table 1, we list Δ_p and Δ₅₀ for all three filters. The second row of Fig. 8 displays the correlation between the bias-corrected sizes and the intrinsic sizes. No correction was done for $R_{50}^{\mathrm{fit}}$. The data points lie closely around the one-to-one relation, marked by a dashed line, indicating that no obvious residual biases exist.

Fractional uncertainty of the size measurement would be more meaningful than the absolute uncertainty, as the measured size is larger in larger galaxies. We calculated the fraction uncertainty as the standard deviation (σ) of difference between measured and intrinsic values divided by the intrinsic values. We plot the fractional uncertainty as a function of R_p, True/FWHM in the third row of Fig. 8. The size measurement becomes more uncertain at lower resolutions. The following exponential functions, respectively, were fitted to the data:

$$\begin{array}{cc}& {\delta }_{R_{\mathrm{p}}}/R_{\mathrm{p}}={ϝ}_1^p·exp({ϝ}_2^p·x),\mathrm{where}x=R_{\mathrm{p},\mathrm{True}}/\mathit{FWHM},\hfill \end{array}$$(17)

$$\begin{array}{cc}& {\delta }_{R_{50}^{\mathrm{cog}}}/R_{50}^{\mathrm{cog}}={ϝ}_1^{50,\mathrm{cog}}·exp({ϝ}_2^{50,\mathrm{cog}}·x),\mathrm{where}x=R_{\mathrm{p},\mathrm{True}}/\mathit{FWHM},\hfill \end{array}$$(18)

and

$$\begin{array}{c}\hfill {\delta }_{R_{50}^{\mathrm{fit}}}/R_{50}^{\mathrm{fit}}={ϝ}_1^{50,\mathrm{fit}}·exp({ϝ}_2^{50,\mathrm{fit}}·x),\mathrm{where}x=R_{\mathrm{p},\mathrm{True}}/\mathit{FWHM}.\end{array}$$(19)

The solid curve in each bottom panel marks the best-fit function. The best-fit parameters for the results based on F115W, F150W, and F200W PSF are listed in Table 1.

The functions used to fit the data, as the above functions and those in the rest of this paper, are chosen so that the function can fit the data without obvious residuals and simultaneously meet our expectations regarding the y-axis values when R_p is sufficiently small or large. These functions are empirical. Our primary objective is not to develop a deep understanding of the underlying physics, but rather to describe and characterize the data itself. Empirical functions are used to create a smooth curve that highlights trends or patterns in the data and then used to correct biases and uncertainties, even if the curve itself does not directly represent any physical reasons.

Although our primary focus is on studying the optical morphology of simulated galaxies at z ≤ 3 observed with filters in the 1.15−2.0 μm range, our data set can also shed light on the optical morphology of galaxies at higher redshifts when observed with redder filters, such as F277W, F356W, and F444W. We used a pixel scale of 0.06 arcsec, created PSFs using WebbPSF, and generated N-FWHM images to carry out the same analysis to comprehend biases and uncertainties involved in quantifying morphology observed through these redder filters. The FWHM values for the F277W, F356W, and F444W PSFs are 0.088, 0.114, and 0.140 arcsec, respectively. However, we did not perform a Sérsic fitting, as the parameters derived from the fitting exhibited no biases. The parameters for deriving the correction functions are shown in Table 2.

4.2. Concentration

The concentration (C) measures the degree to which a galaxy’s light distribution is centrally concentrated. Following the definition in Conselice (2003), we compute C as

$$\begin{array}{c}\hfill C=5·{log}_{10}\left(\frac{R_{80}}{R_{20}}\right)·\end{array}$$(20)

Higher values indicate more centralized light distributions. We denote intrinsic concentration measured in the K-corrected DESI image as C_True. In the first row of Fig. 9, we plot C measured from the N-FWHM images against C_True. The first, second, third, and final columns display the results of image resolution levels R_p, True/FWHM = 1.98, 4.55, 10.45, and 24, respectively. The results show that C is systematically underestimated, with a more significant underestimation for galaxies with higher C_True. This effect is more pronounced in angularly smaller galaxies with lower resolution (lower R_p, True/FWHM). The correlations between C and C_True are nearly linear and, hence, we fit the data with a linear function:

$$\begin{array}{c}\hfill C=k·C_{\mathrm{True}}+b.\end{array}$$(21)

Fig. 9. Evaluation of bias and uncertainty in measuring concentration (C). Top row shows the correlation between C values obtained from downsized images (N-FWHM images) and their intrinsic values, with the solid line representing the best-fit straight line. Middle row displays the correlation between bias-corrected concentration (Ccor) and their intrinsic values. Bottom row presents the slop (k) and intercept (b), and the statistical uncertainty (σC), as a function of Rp, True/FWHM, with the curve marking their best-fit functions.

The best-fit parameters are listed at the top of each panel in the first row and are plotted as a function of R_p, True/FWHM in the bottom row. Toward lower resolutions, the best-fit slop k decreases, while the best-fit intercept b increases. The underestimation of measured C is attributed to the PSF smoothing effect, which overestimates R₂₀ more than R₈₀. Our findings are consistent with Whitney et al. (2021) and Yeom et al. (2017), who showed that C values are more underestimated at higher redshift, where galaxies are smaller and resolutions are lower. Nevertheless, the literature still lacks an accurate correction function to address the issue, which will be further explored in the following discussion.

We use the best-fit straight lines to correct for the bias and uncertainty in measuring C. The correction function is given by:

$$\begin{array}{c}\hfill C_{\mathrm{cor}}=(C-b)/k.\end{array}$$(22)

The correlations between C_cor and C_True are shown in the middle row of Fig. 9. The mean value (Δ_C) and standard deviation (σ_C) of the difference between C_cor and C_True are presented at the top of each panel in the middle row. Also, Δ_C are zero at all examined image resolution levels. σ_C serves as the statistical measurement uncertainty, which increases with lower resolution. The large uncertainty at low resolution stems from the flatten relationship between C and C_Truen or, in other words, the low k value. σ_C as a function of R_p, True/FWHM is plotted in the bottom row.

In practical applications, we observe galaxies of various sizes, which necessitates the development of functions to derive the k and b for a given resolution level. To achieve this goal, we fit k versus R_p, True/FWHM using the following function:

$$\begin{array}{c}\hfill k=\frac{-2·(x+1)}{1+{(x+1)}^g}+1,\mathrm{where}x=R_{\mathrm{p},\mathrm{True}}/\mathit{FWHM}.\end{array}$$(23)

This function converges to 1 when R_p, True is sufficiently large. We fit b versus R_p, True/FWHM using the following function:

$$\begin{array}{c}\hfill b={\xi }_1^b·{(x+{\xi }_2^b)}^{{\xi }_3^b},\mathrm{where}x=R_{\mathrm{p},\mathrm{True}}/\mathit{FWHM}.\end{array}$$(24)

This function converges to 0 when R_p, True is sufficiently large. Therefore, for a given R_p, True, the correction function can be derived using Eqs. (22)–(24). To estimate the statistical uncertainty for a given resolution level, we fit σ_C versus R_p, True/FWHM using the following function:

$$\begin{array}{c}\hfill {\sigma }_C={\eta }_1^C·x^{{\eta }_2^C},\mathrm{where}x=R_{\mathrm{p},\mathrm{True}}/\mathit{FWHM}.\end{array}$$(25)

The best-fit functions are plotted as solid curves. The best-fit parameters g, ${\xi }_1^{\hspace{0.17em}b}$, ${\xi }_2^{\hspace{0.17em}b}$, ${\xi }_3^{\hspace{0.17em}b}$, ${\eta }_1^C$, and ${\eta }_2^C$ for the results based on F115W, F150W, and F200W PSF are listed in Table 1. Those based on F277W, F356W, and F444W are provided in Table 2.

4.3. Asymmetry

4.3.1. Definition

The asymmetry (A) measures the degree to which a galaxy’s light distribution is 180° rotationally symmetric. Originally, A is determined by rotating the galaxy image by 180° about the galaxy center, set as the pixel location of the maximum value, and subtracting the rotated image from the galaxy image to obtain a difference image. Then, A is calculated as 0.5 times the ratio of the summation of the absolute pixel values in the difference image to the summation of the absolute pixel values in the original image (Abraham et al. 1996). The background asymmetry is calculated in the same vein for a portion of sky and subtracted. To solve the problem of substantial variation in A due to uncertainty in the center determination, the algorithm is improved by including a process for searching for a new galaxy center to minimize the asymmetry; background asymmetry is also minimized in the same vein and then subtracted (Conselice et al. 2000; Lotz et al. 2004). In addition, the factor of 0.5 in the original definition is dropped. This is the most conventional method to calculate galaxy asymmetry. The calculation is given by:

$$\begin{array}{c}\hfill \begin{array}{cc}\hfill A_{\mathrm{C}00}& =\frac{min(\sum |I_0-I_{180}|)}{\sum |I_0|}-\frac{min(\sum |B_0-B_{180}|)}{\sum |I_0|}\hfill \\ \hfill & =A_{\mathrm{noisy}}-A_{\mathrm{bkg},\min },\hfill \end{array}\end{array}$$(26)

where I represents the galaxy image, and B represents the sky background. The subscript of A_C00 marks Conselice et al. (2000). The summation is done over all pixels within a 1.5 R_p elliptical aperture, which has the apparent projection parameters determined in Sect. 4.1, centered on the galaxy.

Although the minimization of background asymmetry has been proposed for two decades, there are some studies stick to the background asymmetry without minimization, calculated as follows:

$$\begin{array}{c}\hfill A_{\mathrm{bkg},\mathrm{no}\min }=\frac{\sum |B_0-B_{180}|}{\sum |I_0|}·\end{array}$$(27)

For example, Rodriguez-Gomez et al. (2019) developed the PYTHON package statmorph to quantify galaxy morphology and calculate background asymmetry without minimization within a sky box located beside the galaxy segmentation. Similarly, Tohill et al. (2021) constructed a 10-pixel by 10-pixel grid over the image area outsides the galaxy segmentation to perform the calculation.

However, it has been already shown that A_{bkg,  min} overestimates the contribution of noise to the calculation of galaxy asymmetry (Shi et al. 2009; Wen & Zheng 2016), and so the A_{bkg,  no min} even more severely overestimates the noise contribution. Shi et al. (2009) showed that asymmetry values of galaxies in the field of shallower observations are systematically lower than those of the same galaxies in the field of deeper observations. These authors proposed measuring the distribution of noise asymmetries in randomly selected regions surrounding the target galaxy and calculate the 15% probability low-end tail as the final background asymmetry measurement. Although this method could remove the bias on average, but significant scatter still persists. Wen & Zheng (2016) studied asymmetries of galaxies before and after making them noisy and proposed a new noise correction, which is given by:

$$\begin{array}{c}\hfill A=\frac{min(\sum |I_0-I_{180}|)-F_2·min(\sum |B_0-B_{180}|)}{\sum |I_0|-F_1·\sum |B_0|},\end{array}$$(28)

where

$$\begin{array}{c}\hfill F_1=\frac{N_{I_0<f_1·{\sigma }_{\mathrm{bkg}}}}{N_{\mathrm{all}}},\end{array}$$(29)

and

$$\begin{array}{c}\hfill F_2=\frac{N_{|I_0-I_{180}|<f_2·{\sigma }_{\mathrm{bkg}}}}{N_{\mathrm{all}}},\end{array}$$(30)

where N_all represents the number of pixels encloses by the 1.5 R_p elliptical aperture. N_{I0 < f1 ⋅ σbkg} is the number of pixels dominated by noise in the galaxy image, selected as those with values less than f₁ ⋅ σ_bkg; N_{|I0 − I180|< f2 ⋅ σbkg} is the number of pixels dominated by noise in the difference image, selected as those with values less than f₂ ⋅ σ_bkg. This correction is based on the fact that only noisy pixels are affected. Wen & Zheng (2016) suggested using f₁ = 1 and $f_2=\sqrt{2}$, which, however, are not the optimal choices, as discussed below.

4.3.2. Improved noise correction

To understand the noise contribution, we measured the asymmetry for the resolution-match images (simulated noise not added yet), which are obtained from step 5 of the redshifting procedure and have such high S/N that the noise contribution is almost negligible. We denote the result as A_noise-free. We note that A_noise-free is still biased due to resolution effects and will be addressed later. We then measured the asymmetry for the simulated CEERS images using statmorph, in which no minimization of the background asymmetry was done, and we denote the result as A_noisy − A_{bkg,  no min}. We performed the measurement with minimization using Eq. (26) and denote the result as A_noisy − A_{bkg,  min}. We performed the measurement using Eq. (28), adopting f₁ = 1 and $f_2=\sqrt{2}$, and denote the result as A_WZ16. We plot these as a function of A_noise-free in Fig. 10.

Fig. 10. Comparison of different noise correction methods applied to asymmetry. The four panels present the results for (a) calculating background asymmetry without minimization (statmorph is used), (b) calculating background asymmetry with minimization, (c) Eq. (28) with f1 = 1 and $f_2=\sqrt{2}$, proposed by Wen & Zheng (2016), and (d) Eq. (28) with f1 = 2.25 and f2 = 2.1, our proposed values. Our improved noise correction exhibits enhanced performance in correcting for noise effects and reproducing the noise-free asymmetry.

In a perfect noise correction, A_noise-free should be accurately reproduced. We first confirm in the panel a that A_{bkg,  no min} severely overestimates noise contribution, leading to an underestimation of asymmetry and even physically meaningless negative values. As shown in panel b, A_{bkg,  min} works better than A_{bkg,  no min}, but the overestimation of noise contribution is still significant, especially at higher A_noise-free values. The mean difference between A_noisy − A_{bkg,  min} and A_noise-free gives −0.032. This correlation is sublinear. Panel c shows that A_WZ16 underestimate noise contribution by 0.037, while the main improvement is that the correlation between A_WZ16 and A_noise-free is brought close to linear.

The problem related to the underestimation of noise contribution in using Eq. (28) may be solved if a larger fraction of pixels are defined as noisy pixels. We searched for an optimal solution by testing different value of f₁ and f₂. We calculated the asymmetry using Eq. (28) by adopting f₁ = 0.9, 0.95, 0.1, …, 3.0 and f₂ = 0.9, 0.95, 0.1, …, 3.0. For each pair of f₁ and f₂, we calculated the Pearson correlation coefficient, slop, and scatter between the resulting asymmetry and A_noise-free. The coefficient, |slop − 1|, and scatter as a function of f₁ and f₂ are plotted in Fig. 11.

Fig. 11. Pearson correlation coefficient (ρ), absolute difference between the slop and unity (|slop − 1|), and scatter (σ) between asymmetry calculated using Eq. (28) and Anoise-free as a function of f1 and f2. The values are color-encoded using histogram equalization scale. The black circle marks the optimal f1 = 2.25 and f2 = 2.1 we determine empirically, while the red triangle marks the f1 = 1 and $f_2=\sqrt{2}$, as suggested in Wen & Zheng (2016).

The results show that the coefficient and scatter reach their maximum and minimum values at f₂ ≈ 2.3, while the slop reaches 1 at f₂ ≈ 1.7. The values suggested by Wen & Zheng (2016) are marked by a red triangle, which would get a nearly linear relation, but the relation will be dispersive. Constructing the tightest relation will make the relation sublinear. It is thus not feasible to use Eq. (28) to obtain corrected asymmetry that has a relationship with A_noise-free that is very tight and linear simultaneously. As a compromise, we derive the optimal f₁ and f₂ by requiring that the slop is greater than 0.9 and less than 1.1, and the scatter achieves the minimum value, resulting in f₁ = 2.25 and f₂ = 2.1. Minor changes in the criterion do not significantly impact our results. The optimal choice is marked by black circle in Fig. 11. These values are used to derive noise-removed asymmetry, denoted as A_Improve, which is plotted against A_noise-free in the panel d of Fig. 10. It is shown that the correlation between A_Improve and A_noise-free is tight (ρ = 0.93), almost linear (slop = 0.9), and has a small residual bias (Δ = 0.011) and small scatter (σ = 0.025). The scatter is also smaller than the result obtained by the noise correction proposed by Shi et al. (2009), which yield 0.044 (see their Fig. 15). Our improved noise correction exhibits enhanced performance in correcting for noise effects and reproducing the noise-free asymmetry. Nevertheless, we would like to point out that there is still a minor underestimation of noise contribution present, that is, an overestimation of asymmetry, when A_noise-free ≲ 0.1.

4.3.3. Correction of resolution effects

In addition to the noise contribution, the resolution effect also plays a significant role in affecting asymmetry measurements. This effect becomes particularly crucial when studying high-redshift galaxies, which are less spatially resolved. Focusing on changes in the apparent galaxy size relative to a fixed pixel size, Conselice et al. (2000) demonstrated that measured galaxy asymmetry is increasingly reduced, as the apparent size of a galaxy decreases. We measured galaxy asymmetries using re-binned images without PSF convolution, and found that the binning effect is negligible compared to the PSF effect. It has also been found that galaxy asymmetries are underestimated at higher redshifts (Conselice 2003; Yeom et al. 2017; Whitney et al. 2021); however, this should be partly attributed to the overestimation of noise contribution using the conventional method, as discussed in Sect. 4.3.2, and partly to resolution effects. An accurate correction method for biases caused by resolution effects has not been developed yet.

To account for resolution effects, we first measured intrinsic asymmetries from the high-quality, K-corrected DESI images and denote them as A_True. We then measured A from the N-FWHM images and present A as a function of A_True in the top row of Fig. 12. The first, second, third, and final columns display results obtained at the image resolution levels R_p, True/FWHM = 1.98, 4.55, 10.45, and 24, respectively. The results show that when galaxies are the most intrinsically symmetric (A_True ≲ 0.05) and are at low resolution (R_p, True/FWHM ≲ 4.55), the measured A slightly overestimates A_True. This overestimation is attributed to the asymmetry of the JWST PSF. Figure 13 displays the two-time oversampling F200W PSF on the left and the difference between the PSF and its 180°-rotational image on the right. The PSF asymmetry causes the convolution to make intrinsically symmetric galaxies appear asymmetric. In contrast, when galaxies are intrinsically asymmetric (A_True ≳ 0.05), the measured A values are underestimated, with a greater extent at higher A_True. This underestimation is due to PSF convolution smoothing out asymmetric structures, particularly the most asymmetric ones. The smoothing effect is more efficient than the effect caused by PSF asymmetry in asymmetric galaxies.

Fig. 12. Evaluation of bias and uncertainty due to resolution effects in measuring asymmetry (A). Top row shows the correlation between A values obtained from the N-FWHM images and their intrinsic values (ATrue), with the solid curve representing the best-fit function. Middle row displays the correlation between bias-corrected asymmetries (Acor) and ATrue. Bottom row presents the parameters (T, D, and A0) and the uncertainty in correcting A (σA), as a function of Rp, True/FWHM, with the curve marking their best-fit functions.

Fig. 13. Illustration of the asymmetry of two-time oversampling JWST F200W PSF. The PSF is plotted in logarithmic scale on the left, while the difference between the PSF and its 180°-rotational image is plotted in linear scale on the right.

The A_True − A relations are nonlinear, especially at low resolutions. We fit the data at each resolution level with the following function:

$$\begin{array}{c}\hfill A=A_{\mathrm{True}}·{(A_{\mathrm{True}}/D)}^T+A_0,\end{array}$$(31)

where D and T control the flatness of the function and A₀ is the y-intercept representing the asymmetry of a fully symmetric galaxy after convolving a asymmetric PSF. The best-fit functions are plotted as black solid curves. We reverse Eq. (31) to estimate the intrinsic value for a given measured A value. However, this operation is only available for A ≥ A₀. To estimate the intrinsic value for A < A₀, we replace A with 2A − A₀, which is the symmetric value of A with respect to y = A₀, to do the estimation. The correction function is as follows:

$$\begin{array}{c}\hfill A_{\mathrm{cor}}=D·{\left(\frac{\mathrm{\Delta }A}{D}\right)}^{\frac{1}{1+T}},\mathrm{where}\mathrm{\Delta }A=\{\begin{array}{c}A-A_0,\mathrm{if}A\ge A_0\hfill \\ 2A_0-A,\mathrm{if}A<A_0\hfill \end{array}.\end{array}$$(32)

The second row in Fig. 12 plots the correlations between A_cor and A_True. The mean difference (Δ_A) and scatter (σ_A) between them are presented at the top of each panel. At high resolution, Δ_A is small (|Δ_A|< 0.01). However, At low resolution, a non-negligible Δ_A value (∼0.04) is observed, indicating that a residual small bias still exists when the resolution is sufficiently low. Meanwhile, σ_A is significant (∼0.1) at R_p, True/FWHM < 4.55, which stems from the flatness of the A_True − A relation.

To obtain the correction function for a given resolution level, we study the parameters T, D, and A₀ as a function of R_p, True/FWHM in the bottom row of Fig. 12. We respectively fit the correlations between T, D, and A₀ with R_p, True/FWHM using the functions:

$$\begin{array}{cc}& T={\xi }_1^T·{(x+{\xi }_2^T)}^{{\xi }_3^T},\mathrm{where}x=R_{\mathrm{p},\mathrm{True}}/\mathit{FWHM},\hfill \end{array}$$(33)

$$\begin{array}{cc}& D={\tau }_1·exp({\tau }_2·x)+{\tau }_3·ln(x+{\tau }_4),\mathrm{where}x=R_{\mathrm{p},\mathrm{True}}/\mathit{FWHM},\hfill \end{array}$$(34)

and

$$\begin{array}{c}\hfill A_0={\xi }_1^{A_0}·{(x+{\xi }_2^{A_0})}^{{\xi }_3^{A_0}},\mathrm{where}x=R_{\mathrm{p},\mathrm{True}}/\mathit{FWHM}.\end{array}$$(35)

The best-fit functions are marked with black solid curves. The correction function Eq. (32) can therefore be derived for a given galaxy size through Eqs. (33)–(35). To estimate the statistical uncertainty for a given resolution level, we plot σ_A against R_p, True/FWHM in Fig. 12 and fit them using the following function:

$$\begin{array}{c}\hfill {\sigma }_A={\eta }_1^A·x^{{\eta }_2^A},\mathrm{where}x=R_{\mathrm{p},\mathrm{True}}/\mathit{FWHM}.\end{array}$$(36)

The best-fit function is plotted as a black solid curve. The best-fit parameters ${\xi }_1^T$, ${\xi }_2^T$, ${\xi }_3^T$, τ₁, τ₂, τ₃, τ₄, ${\xi }_1^{A_0}$, ${\xi }_2^{A_0}$, ${\xi }_3^{A_0}$, ${\eta }_1^A$, and ${\eta }_2^A$ for the results based on F115W, F150W, and F200W PSF are listed in Table 1. Those based on F277W, F356W, and F444W are shown in Table 2.

4.4. Sérsic index and axis ratio

We use IMFIT (Erwin 2015) to carry out a two-dimensional (2D) fitting using a single Sérsic model (Sersic 1968) on the images. This allows us to determine the half-light radius ($R_{50}^{\mathrm{fit}}$), Sérsic index (n), and axis ratio (q). We use a two-time oversampling PSFs in our analysis. IMFIT finds the optimal model by adjusting the 2D function parameters through nonlinear minimization of total χ². The Levenberg-Marquardt algorithm is used for the χ² minimization. We set the lower and upper bounds of n for the fitting to 0.5 and 6, respectively. The bias and uncertainty in measuring $R_{50}^{\mathrm{fit}}$ have been discussed in Sect. 4.1.

We denote intrinsic n and q values measured from the high-quality, K-corrected DESI images as n_True and q_True, respectively. We plot the difference between n measured from the N-FWHM image and n_True, that is n − n_True, as a function of R_p, True/FWHM in the first row of Fig. 14. The mean difference, Δ_n = −0.11, indicates that the impact of resolution degradation is small, properly because IMFIT already accounts for PSF effects. As we go on to show in Sect. 5, there is no obvious difference on average between intrinsic values and those measured on the simulated CEERS images. Our findings are in line with previous studies based on model galaxies (Barden et al. 2008; Davari et al. 2016; Euclid Collaboration 2023; but see Paulino-Afonso et al. 2017). We thus refrain from developing a correction function for n. The second row plots the profile of measurement uncertainty σ_n, which we fit with the function:

$$\begin{array}{c}\hfill {\sigma }_n={\eta }_1^n·x^{{\eta }_2^n},\mathrm{where}x=R_{\mathrm{p},\mathrm{True}}/\mathit{FWHM}.\end{array}$$(37)

Fig. 14. Evaluation of biases and uncertainties due to resolution effects in measuring Sésic index (n) and axis ratio (q) using IMFIT.

The black curve represents the best-fit function. The measurement of n becomes more uncertain with lower resolutions.

Similarly, we calculate q − q_True and plot it as a function of R_p, True/FWHM in the third row of Fig. 14. The mean difference is very small (Δ_q = −0.005), indicating little or no measurement bias. The bottom panel presents the profile of statistical uncertainty σ_q, which we fit with the function:

$$\begin{array}{c}\hfill {\sigma }_q={\eta }_1^q·x^{{\eta }_2^q},\mathrm{where}x=R_{\mathrm{p},\mathrm{True}}/\mathit{FWHM}.\end{array}$$(38)

The black curve marks the best-fit function. The measurement of q becomes more uncertain at lower resolutions, yet these uncertainties remain negligible when compared to the broad dynamical range of q. The best-fit parameters ${\eta }_1^n$, ${\eta }_2^n$, ${\eta }_1^q$, and ${\eta }_2^q$ for the results based on F115W, F150W, and F200W PSFs are provided in Table 1.

5. Application to simulated CEERS images

In this section, we describe how we apply the correction functions to the morphological quantities measured from the simulated CEERS images to understand the effectiveness of these functions. We start by measuring R_p, $R_{50}^{\mathrm{cog}}$, $R_{50}^{\mathrm{fit}}$, C, A_C00, n, and q on the simulated CEERS images. The corrections are summarized as follows. We correct size measurements using Eqs. (15) and (16) to obtain R_p, cor and $R_{50,\mathrm{cor}}^{\mathrm{cog}}$, respectively. The resolution level is calculated as R_p, cor/FWHM. We obtain the parameters k and b through Eqs. (23) and (24), respectively, and then calculate the bias-corrected concentration (C_cor) using the correction function Eq. (22). We calculate the noise-removed asymmetry using Eq. (28), with f₁ = 2.25 and f₂ = 2.1, obtain the parameters T, D, and A₀ through Eqs. (33)–(35), respectively, and, finally, we calculate the bias-corrected asymmetry (A_cor) using the correction function Eq. (32). There is no correction for the parameters ($R_{50}^{\mathrm{fit}}$, n, and q) measured through Sérsic fitting using IMFIT. The statistical uncertainties are computed using the uncertainty function for each morphological measurement, but the derived uncertainties are not plotted in the figures of the following discussion for the sake of clarity.

In Fig. 15, we plot R_p and R_p, cor as a function of R_p, True in the top two rows, $R_{50}^{\mathrm{cog}}$ and $R_{50,\mathrm{cor}}^{\mathrm{cog}}$ as a function of $R_{50,\mathrm{True}}^{\mathrm{cog}}$ in the middle two rows, and $R_{50}^{\mathrm{fit}}$ as a function of $R_{50,\mathrm{True}}^{\mathrm{fit}}$ in the bottom row. Results for z = 0.75, 1.5, 2.25, and 3.0 are presented in columns one through four, respectively, and the difference and scatter between y-axis and x-axis values are indicated in the top of each panel (Figs. 16 and 17 follow the same strategy). Values of R_p and $R_{50}^{\mathrm{cog}}$ slightly overestimate their intrinsic values by 0.06 arcsec and 0.03 arcsec, respectively. After bias correction, these overestimations are reduced to 0.01 arcsec, which is negligibly small. The correlations between $R_{50}^{\mathrm{fit}}$ and $R_{50,\mathrm{True}}^{\mathrm{fit}}$ exhibit a very small average offset, but they display a slightly larger scatter compared with the non-parametric measures of galaxy sizes.

Fig. 15. Comparison of galaxy sizes measured from simulated CEERS images, bias-corrected sizes, and intrinsic sizes. The top two rows present results for Petrosian radius (Rp). The middle two rows present results for half-light radius obtained through non-parametric method ($R_{50}^{\mathrm{cog}}$). The bottom row presents results for half-light radius obtained through Sérsic fitting ($R_{50}^{\mathrm{fit}}$), with no correction applied. The mean difference (Δ) and scatter (σ) between the x-axis and y-axis values are indicated at the top of each panel. Columns one through four show results for z = 0.75, 1.5, 2.25, and 3.0, respectively.

The top two rows of Fig. 16 presents the plots of C and C_cor as a function of C_True. CEERS galaxies with higher C_True tend to have measured C to be underestimated, and since the galaxies are angularly smaller at higher redshifts, the underestimation becomes more severe. Scatters are especially large at high redshifts. As early-type galaxies are more centrally concentrated than late-type galaxies, the C values of early-type galaxies will be more underestimated without correction. After bias correction, C_cor correlates well with C_True, exhibiting a small average offset of approximately −0.02 and with an average scatter of around 0.22.

Fig. 16. Evaluation of the effectiveness of the correction functions in reproducing the intrinsic concentration and asymmetry. The top two rows present comparisons of the concentration measured from simulated CEERS images (C), bias-corrected values (Ccor), and intrinsic values (CTrue). The bottom two rows illustrate comparisons of asymmetry measured from simulated CEERS images using the conventional method (AC00), bias-corrected values (Acor), and intrinsic values (ATrue). Columns one through four show results for z = 0.75, 1.5, 2.25, and 3.0, respectively. The mean difference (Δ) and scatter (σ) between the y-axis and x-axis values are indicated at the top of each panel. In the ATrue − Acor relations, data are divided into two groups: angularly small galaxies with Rp, True/FWHM < 5 (red circles) and angularly large galaxies with Rp, cor/FWHM ≥ 5 (blue crosses).

The bottom two rows of Fig. 16 presents the plots of A_C00 and A_cor as a function of A_True. Our results show that A_C00 underestimate A_True on average, especially at high A_True values or at high redshifts, due to the overestimation of noise contribution and the effects of PSF smoothing. Since late-type galaxies have higher intrinsic asymmetry than early-type galaxies (Abraham et al. 1996), the measured asymmetry for late-type galaxies (if not corrected) will be more severely underestimated compared to that of early-type galaxies. Our findings suggest that employing the conventional algorithm for calculating galaxy asymmetry, as proposed in Conselice et al. (2000), is unsuitable and inadequate to measure the intrinsic asymmetry of high-redshift galaxies observed in CEERS. In a recent study, Kartaltepe et al. (2023) visually classified the morphological type of 850 galaxies at z > 3 observed in JWST CEERS and used Statmorph to perform morphological measurements. They used images from the filter corresponding to the rest-frame optical emission at the redshift of the galaxy. Their results reveal that the Concentration–Asymmetry diagram does not clearly differentiate between morphological types, although peculiar galaxies exhibit slightly higher asymmetry, on average. Our findings provide a possible explanation to the results in Kartaltepe et al. (2023): the CEERS PSF and noise introduce significant bias and uncertainty to the measurements of concentration and asymmetry, causing the resulting concentration–asymmetry diagram to be degenerate across galaxy types.

The A_True − A_cor relations demonstrate that our improved method can accurately reproduce the intrinsic galaxy asymmetry when z ≲ 1.5. However, when z ≳ 1.5, the relations become increasingly dispersive with higher redshifts. We divide the simulated CEERS galaxies into two groups: angularly small galaxies with R_p, True/FWHM < 5, marked with red circles, and angularly large galaxies with R_p, cor/FWHM ≥ 5, marked with blue crosses. The results show that, for angularly large galaxies, A_cor reproduces A_True well with a mean difference of 0.01 across all redshifts, and a scatter increasing from 0.03 at z = 0.75 to 0.06 at z = 3.0. In contrast, for angular small galaxies, the A_cor still overestimates A_True with a mean difference of approximately 0.1 across all redshifts and with scatter increasing from 0.06 at z = 0.75 to 0.09 at z = 3.0. The latter is caused by the incomplete removal of noise correction for the most symmetric galaxies, even when using our improved noise correction (Fig. 10), and by the flatness of the correction function (Fig. 12). Therefore, our asymmetry correction function is only efficient for angularly large galaxies (R_p, cor/FWHM ≥ 5) observed in JWST CEERS.

In Fig. 17, we plot n against n_True in the top row and q against q_True in the bottom row. The average difference between n and n_True is small, approximately −0.03, with an average scatter of around 0.5. This suggest that the Sérsic fitting can extract the index without significant bias for galaxies observed in JWST CEERS, but the statistical uncertainty should be appropriately considered. The q − q_True relations are quite tight with negligible offset and scatter, indicating that the Sérsic fitting can robustly extract the axis ratio.

Fig. 17. Evaluation of the effectiveness of the correction functions in reproducing the intrinsic Sérsic index and axis ratio. Top: comparison between Sérsic index measured from the simulated CEERS images (n) and the intrinsic values (nTrue). Bottom: comparison between axis ratio measured from the simulated CEERS images (q) and the intrinsic values (qTrue). The mean difference (Δ) and scatter (σ) between the y-axis and x-axis values are indicated at the top of each panel.

6. Summary and conclusions

Early JWST studies have shown that galaxies with established disk and spheroidal morphologies span the full redshift range, revealing that the Hubble sequence was already in place at the early universe (e.g., Ferreira et al. 2022a,b; Kartaltepe et al. 2023; Nelson et al. 2023; Robertson et al. 2023). However, potential biases and uncertainties in characterizing the morphologies of these high-redshift galaxies may significantly impact the results. To address this issue, we defined a sample of nearby galaxies based on the DESI survey, conststing of 1816 galaxies with stellar masses ranging from 10^9.75 M_⊙ to 10^11.25 M_⊙. High-quality DESI images of these nearby galaxies allow for an accurate determination of the galaxy morphology. We removed the contamination from foreground stars or background galaxies in the DESI images and used the resulting cleaned images of the g, r, and z bands as input to compute artificial images of galaxies located at 0.75 ≤ z ≤ 3 and observed at rest-frame optical wavelengths in CEERS. We used the F115W filter for z = 0.75 and 1, the F150W filter for z = 1.25, 1.5, and 1.75, and the F200W filter for z = 2.0 to 3.0. For the artificially redshifting process, we take into account angular size changes due to distance and intrinsic evolution, flux changes due to cosmological dimming and intrinsic evolution, spectral change, and changes in resolution and noise level. The rest-frame flux is obtained by pixel-by-pixel K correction. The simulated images are binned onto the pixel scale of 0.03 arcsec/pixel and corrected to the JWST PSF. A patch of real CEERS background and poisson noise from galaxy light are then added.

We focus on quantifying the biases and uncertainties for six widely used morphological quantities: Petrosian radius (R_p), half-light radius (R₅₀), asymmetry (A), concentration (C), axis ratio (q), and Sérsic index (n). The resolution of CEERS images emerges as the primary factor influencing these measurements, while the typical CEERS noise predominantly impacts the computation of A. We improve the method for removing the contribution of noise from the computation of A. To correct for the resolution effects, we reduce the DESI image resolution, conduct each measurement, and compare these re-measured values with their intrinsic values to understand resolution effects and to derive formulae as a function of resolution level, defined as R_p/FWHM, for correcting the biases and uncertainties. Finally, we apply the correction functions to our artificially redshifted CEERS images to validate the methods. Our main results are as follows.

1. R_p and R₅₀, measured using non-parametric approaches, are slightly overestimated due to PSF smoothing, and this overestimation does not significantly depends on the resolution level. We derived Eqs. (15) and (16) to correct these biases. In comparison, R₅₀ measured using 2D image fitting proves to be unbiased. The functions of statistical uncertainties for these three parameters are provided.

2. C is underestimated due to PSF smoothing, with the effect being more pronounced for higher C values or lower resolutions. The bias can be corrected using correction function of Eq. (22) and the statistical uncertainty is given by Eq. (25).

3. By incorporating a more accurate noise effect removal procedure, we improved the computation of A over existing methods, which can often overestimate, underestimate, or lead to significant scatter of noise contributions. We show that A of the most intrinsically symmetric galaxies are overestimated due to the PSF asymmetry. In contrast, A of intrinsically asymmetric galaxies are underestimated owing to smoothing, particularly for large A values and at lower resolutions. For angularly large CEERS galaxies where R_p/FWHM ≥ 5, the biases can be robustly corrected using Eq. (32). However, for smaller galaxies, these biases cannot be completely removed. When studying asymmetry, the statistical uncertainty given by Eq. (36) should be taken into account.

4. The measurements of n and q through 2D image fitting have negligible biases. The statistical uncertainty for axis ratio is negligible, whereas the uncertainty for the Sérsic index is more significant and should be properly considered.

5. Although our primary focus is on studying the optical morphology of simulated galaxies at z ≤ 3 observed with F115W, F150W, and F200W filters, we also provide the correction for F277W, F356W, and F444W filters, which may be useful for studying optical morphology of galaxies at higher redshifts. The parameters of the correction functions vary across different filters and are provided in Tables 1 and 2.

These tests establish a solid foundation for future quantitative statistical studies aimed at comprehending the cosmological evolution of galaxy morphology. The dataset of artificially redshifted images also holds significant value for additional studies, such as those focused on spiral arms and bars.

Acknowledgments

We are grateful to the anonymous referee for their invaluable feedback and insightful comments that greatly improved the quality of this paper. S.Y.Y. acknowledges the support by the Alexander von Humboldt Foundation. S.Y.Y. thanks Luis C. Ho for inspiring him to pursue research on high-redshift galaxies. We thank the fruitful discussion with John Moustakas. We thank Song Huang for creating the Slack workspace that enabled SYY to find collaboration with CC and FS. This research made use of Photutils, an Astropy package for detection and photometry of astronomical sources (Bradley et al. 2022). We acknowledge the usage of the HyperLeda database (http://leda.univ-lyon1.fr). This work is based on observations taken by the 3D-HST Treasury Program (GO 12177 and 12328) with the NASA/ESA HST, which is operated by the Association of Universities for Research in Astronomy, Inc., under NASA contract NAS5-26555. The Legacy Surveys consist of three individual and complementary projects: the Dark Energy Camera Legacy Survey (DECaLS; Proposal ID #2014B-0404; PIs: David Schlegel and Arjun Dey), the Beijing-Arizona Sky Survey (BASS; NOAO Prop. ID #2015A-0801; PIs: Zhou Xu and Xiaohui Fan), and the Mayall z-band Legacy Survey (MzLS; Prop. ID #2016A-0453; PI: Arjun Dey). DECaLS, BASS and MzLS together include data obtained, respectively, at the Blanco telescope, Cerro Tololo Inter-American Observatory, NSF’s NOIRLab; the Bok telescope, Steward Observatory, University of Arizona; and the Mayall telescope, Kitt Peak National Observatory, NOIRLab. Pipeline processing and analyses of the data were supported by NOIRLab and the Lawrence Berkeley National Laboratory (LBNL). The Legacy Surveys project is honored to be permitted to conduct astronomical research on Iolkam Du’ag (Kitt Peak), a mountain with particular significance to the Tohono O’odham Nation. NOIRLab is operated by the Association of Universities for Research in Astronomy (AURA) under a cooperative agreement with the National Science Foundation. LBNL is managed by the Regents of the University of California under contract to the US Department of Energy. This project used data obtained with the Dark Energy Camera (DECam), which was constructed by the Dark Energy Survey (DES) Collaboration. Funding for the DES Projects has been provided by the US Department of Energy, the US National Science Foundation, the Ministry of Science and Education of Spain, the Science and Technology Facilities Council of the United Kingdom, the Higher Education Funding Council for England, the National Center for Supercomputing Applications at the University of Illinois at Urbana-Champaign, the Kavli Institute of Cosmological Physics at the University of Chicago, Center for Cosmology and Astro-Particle Physics at the Ohio State University, the Mitchell Institute for Fundamental Physics and Astronomy at Texas A&M University, Financiadora de Estudos e Projetos, Fundacao Carlos Chagas Filho de Amparo, Financiadora de Estudos e Projetos, Fundacao Carlos Chagas Filho de Amparo a Pesquisa do Estado do Rio de Janeiro, Conselho Nacional de Desenvolvimento Cientifico e Tecnologico and the Ministerio da Ciencia, Tecnologia e Inovacao, the Deutsche Forschungsgemeinschaft and the Collaborating Institutions in the Dark Energy Survey. The Collaborating Institutions are Argonne National Laboratory, the University of California at Santa Cruz, the University of Cambridge, Centro de Investigaciones Energeticas, Medioambientales y Tecnologicas-Madrid, the University of Chicago, University College London, the DES-Brazil Consortium, the University of Edinburgh, the Eidgenossische Technische Hochschule (ETH) Zurich, Fermi National Accelerator Laboratory, the University of Illinois at Urbana-Champaign, the Institut de Ciencies de l’Espai (IEEC/CSIC), the Institut de Fisica d’Altes Energies, Lawrence Berkeley National Laboratory, the Ludwig Maximilians Universitat Munchen and the associated Excellence Cluster Universe, the University of Michigan, NSF’s NOIRLab, the University of Nottingham, the Ohio State University, the University of Pennsylvania, the University of Portsmouth, SLAC National Accelerator Laboratory, Stanford University, the University of Sussex, and Texas A&M University. BASS is a key project of the Telescope Access Program (TAP), which has been funded by the National Astronomical Observatories of China, the Chinese Academy of Sciences (the Strategic Priority Research Program “The Emergence of Cosmological Structures” Grant # XDB09000000), and the Special Fund for Astronomy from the Ministry of Finance. The BASS is also supported by the External Cooperation Program of Chinese Academy of Sciences (Grant # 114A11KYSB20160057), and Chinese National Natural Science Foundation (Grant # 12120101003, # 11433005). The Legacy Survey team makes use of data products from the Near-Earth Object Wide-field Infrared Survey Explorer (NEOWISE), which is a project of the Jet Propulsion Laboratory/California Institute of Technology. NEOWISE is funded by the National Aeronautics and Space Administration. The Legacy Surveys imaging of the DESI footprint is supported by the Director, Office of Science, Office of High Energy Physics of the US Department of Energy under Contract No. DE-AC02-05CH1123, by the National Energy Research Scientific Computing Center, a DOE Office of Science User Facility under the same contract; and by the US National Science Foundation, Division of Astronomical Sciences under Contract No. AST-0950945 to NOAO. The Siena Galaxy Atlas was made possible by funding support from the US Department of Energy, Office of Science, Office of High Energy Physics under Award Number DE-SC0020086 and from the National Science Foundation under grant AST-1616414.

References

All Tables

All Figures

	Fig. 1. Basic properties of our sample. Distribution of (a) luminosity distances, (b) absolute r-band magnitude, corrected for Galactic extinction, (c) stellar mass, (d) B-band isophotal diameter at μB = 25 mag arcsec−2, and (e) Hubble types.
In the text

Fig. 2. Demonstration of our cleaning process. Left and right columns showcase the r-band observation of two galaxies, ESO 121−026 (a Sbc galaxy) and ESO 251−004 (an elliptical galaxy). Top and bottom rows display the image before and after undergoing star removal, respectively. The vast majority of emission from foreground stars or background galaxies in the vicinity of the galaxies was successfully removed.

In the text

	Fig. 3. Redshift evolution of intrinsic absolute magnitude at rest-frame waveband of U, B, V, R, I, and J for early-type galaxies. The measured parameter α (Eq. (14)) describing the evolution rate for a given mass range and a given rest-frame waveband is shown in the bottom-right corner.
In the text

	Fig. 4. Same details as in Fig. 3, but for late-type galaxies.
In the text

	Fig. 5. Best-fit α as a function of filter effective wavelength (λeff) and stellar mass (M*) for early-type galaxies. The color bar encodes α.
In the text

	Fig. 6. Same details as in Fig. 5, but for late-type galaxies.
In the text

	Fig. 7. Examples of artificially redshifted galaxy images and the comparison between simulated and real CEERS images. From left to right: first four panels in each row show the DESI r-band image and the artificial JWST CEERS images at z = 1, z = 2, and z = 3. The fifth column presents real galaxy images in JWST CEERS for comparison. Apparent scales are indicated in each panel.
In the text

Fig. 8. Evaluation of biases and uncertainties in measuring Petrosian radius (Rp; left), half-light radius using a non-parametric method ($R_{50}^{\mathrm{cog}}$; middle), and half-light radius through Sérsic model fitting ($R_{50}^{\mathrm{fit}}$; right). The top row illustrates the difference between the measured and intrinsic values as a function of Rp, True/FWHM, with the solid horizon line marking the mean difference. The middle row displays the correlation between the corrected and intrinsic values, with the one-to-one relation indicated by a dashed line. The bottom row presents the fractional uncertainties, together with their best-fit functions plotted as solid lines.

In the text

Fig. 9. Evaluation of bias and uncertainty in measuring concentration (C). Top row shows the correlation between C values obtained from downsized images (N-FWHM images) and their intrinsic values, with the solid line representing the best-fit straight line. Middle row displays the correlation between bias-corrected concentration (Ccor) and their intrinsic values. Bottom row presents the slop (k) and intercept (b), and the statistical uncertainty (σC), as a function of Rp, True/FWHM, with the curve marking their best-fit functions.

In the text

Fig. 10. Comparison of different noise correction methods applied to asymmetry. The four panels present the results for (a) calculating background asymmetry without minimization (statmorph is used), (b) calculating background asymmetry with minimization, (c) Eq. (28) with f1 = 1 and $f_2=\sqrt{2}$, proposed by Wen & Zheng (2016), and (d) Eq. (28) with f1 = 2.25 and f2 = 2.1, our proposed values. Our improved noise correction exhibits enhanced performance in correcting for noise effects and reproducing the noise-free asymmetry.

In the text

Fig. 11. Pearson correlation coefficient (ρ), absolute difference between the slop and unity (|slop − 1|), and scatter (σ) between asymmetry calculated using Eq. (28) and Anoise-free as a function of f1 and f2. The values are color-encoded using histogram equalization scale. The black circle marks the optimal f1 = 2.25 and f2 = 2.1 we determine empirically, while the red triangle marks the f1 = 1 and $f_2=\sqrt{2}$, as suggested in Wen & Zheng (2016).

In the text

Fig. 12. Evaluation of bias and uncertainty due to resolution effects in measuring asymmetry (A). Top row shows the correlation between A values obtained from the N-FWHM images and their intrinsic values (ATrue), with the solid curve representing the best-fit function. Middle row displays the correlation between bias-corrected asymmetries (Acor) and ATrue. Bottom row presents the parameters (T, D, and A0) and the uncertainty in correcting A (σA), as a function of Rp, True/FWHM, with the curve marking their best-fit functions.

In the text

	Fig. 13. Illustration of the asymmetry of two-time oversampling JWST F200W PSF. The PSF is plotted in logarithmic scale on the left, while the difference between the PSF and its 180°-rotational image is plotted in linear scale on the right.
In the text

	Fig. 14. Evaluation of biases and uncertainties due to resolution effects in measuring Sésic index (n) and axis ratio (q) using IMFIT.
In the text

Fig. 15. Comparison of galaxy sizes measured from simulated CEERS images, bias-corrected sizes, and intrinsic sizes. The top two rows present results for Petrosian radius (Rp). The middle two rows present results for half-light radius obtained through non-parametric method ($R_{50}^{\mathrm{cog}}$). The bottom row presents results for half-light radius obtained through Sérsic fitting ($R_{50}^{\mathrm{fit}}$), with no correction applied. The mean difference (Δ) and scatter (σ) between the x-axis and y-axis values are indicated at the top of each panel. Columns one through four show results for z = 0.75, 1.5, 2.25, and 3.0, respectively.

In the text

Fig. 16. Evaluation of the effectiveness of the correction functions in reproducing the intrinsic concentration and asymmetry. The top two rows present comparisons of the concentration measured from simulated CEERS images (C), bias-corrected values (Ccor), and intrinsic values (CTrue). The bottom two rows illustrate comparisons of asymmetry measured from simulated CEERS images using the conventional method (AC00), bias-corrected values (Acor), and intrinsic values (ATrue). Columns one through four show results for z = 0.75, 1.5, 2.25, and 3.0, respectively. The mean difference (Δ) and scatter (σ) between the y-axis and x-axis values are indicated at the top of each panel. In the ATrue − Acor relations, data are divided into two groups: angularly small galaxies with Rp, True/FWHM < 5 (red circles) and angularly large galaxies with Rp, cor/FWHM ≥ 5 (blue crosses).

In the text

Fig. 17. Evaluation of the effectiveness of the correction functions in reproducing the intrinsic Sérsic index and axis ratio. Top: comparison between Sérsic index measured from the simulated CEERS images (n) and the intrinsic values (nTrue). Bottom: comparison between axis ratio measured from the simulated CEERS images (q) and the intrinsic values (qTrue). The mean difference (Δ) and scatter (σ) between the y-axis and x-axis values are indicated at the top of each panel.

In the text