© The Authors 2026

1 Introduction

One of the most fundamental objectives of astronomy is to map and understand the contents of the Universe. Our knowledge of the cosmos has grown in lock-step with technological progress. Extragalactic surveys have closely followed this explosive growth in technological capability, from the first deep surveys with photographic plates in the 1970s to ground-based “deep-field” imaging with electronic detectors, which routinely reached B ~ 28 mag at the end of the twentieth century. Each technological advance opened a new window in wavelength, depth, or resolution. Space-based observatories, starting with the Hubble Space Telescope (HST) pioneered a series of surveys starting with the Hubble Deep Field North (HDF-N; Williams et al. 1996) and continuing with the Great Observatories Origins Deep Survey (GOODS; Giavalisco et al. 2004), the Hubble Ultra-Deep Field (HUDF; Beckwith et al. 2006; Ellis et al. 2013; Illingworth et al. 2013; Teplitz et al. 2013) and the Cosmic Assembly Near-infrared Deep Extragalactic Legacy Survey (CANDELS; Grogin et al. 2011; Koekemoer et al. 2011). These data revealed a distant Universe that was much more complex than had been seen in low-resolution ground-based observations. These surveys, often combined with ground-based spectroscopic measurements, revealed the cosmic history of star formation, reionization, and stellar mass assembly, together with the transformation of the physical and morphological properties of galaxies (e.g., Madau & Dickinson 2014; Conselice 2014; Robertson 2022, for reviews).

However, the majority of these surveys, although reaching great depths, only covered small areas, totalling only 0.2 deg² (800 arcmin²). Despite the discovery of some remarkable objects (e.g., Oesch et al. 2016) these surveys were limited in their ability to probe a large range of galaxy populations and environments. Moreover, their volume at lower redshifts is small, meaning that it is impossible to trace the evolution of a well-defined population of objects over a significant range in cosmic time.

The Cosmic Evolution Survey (COSMOS; Scoville et al. 2007) was designed to bridge the gap between shallow, wide-area ground-based surveys and deep, pencil-beam space-based surveys. Thanks to one of the largest ever allocations of HST time, a contiguous 2 deg² patch of the sky was covered with the Advanced Camera for Surveys (ACS) in the F814W band (Koekemoer et al. 2007). At the same time, a series of complementary observing programs were started to provide the necessary wavelength coverage across the electromagnetic spectrum. COSMOS demonstrated the feasibility of weak lensing measurements in tomographic redshift bins to map the distribution of dark matter over large areas (Massey et al. 2007), enabled the mapping of large-scale structures at intermediate redshifts (Scoville et al. 2013), and demonstrated how mass and environment drive the evolution of galaxies across cosmic time (Peng et al. 2010).

The James Webb Space Telescope (JWST) represents the next major technological step forward, providing orders-of-magnitude increases in sensitivity in 1–5 μm with the Near Infrared Camera (NIRCam; Rieke et al. 2023a) and 6–30 μm from the Mid Infrared Instrument (MIRI; Wright et al. 2022) delivering an unprecedented combination of sensitivity, spatial resolution, and field of view. One of the largest contiguous JWST programs, in both area and time allocated, is the COSMOS-Web survey (GO#1727, PIs: Kartaltepe & Casey, Casey et al. 2023). COSMOS-Web covers the central part of the COSMOS field, covering ~0.54 deg² (1920 arcmin²) in four NIRCam filters (F115W, F150W, F277W, F444W) and ~0.2 deg² (720 arcmin²) in one MIRI filter (F770W). This combination of depth and area, together with existing multiwavelength data, makes the COSMOS-Web a unique survey to study galaxy evolution across a large range of cosmic history, capturing some of the rarest and most extreme objects in the Universe.

To optimally exploit the rich multiwavelength data in COSMOS with JWST imaging from COSMOS-Web, traditional fixed-aperture photometric techniques are impractical. The point-spread function (PSF) in COSMOS ranges from 1″ in ground-based bands to 0.1″ in NIRCam. However, techniques where a model profile is fitted to each source after convolution with the PSF are particularly well suited to COSMOS data. These methods enable consistent measurement of total fluxes for sources at different wavelengths observed with different PSFs. In the COSMOS2020 catalog (Weaver et al. 2022), these methods for deep multiwavelength photometry were validated by comparing measurements made using the profile fitting code The Farmer (Weaver et al. 2023) with the CLASSIC aperture photometry. These showed excellent agreement with the advantage in the profile fitting method that the data do not need to be PSF-homogenized before measurement.

Today, building on the legacy of previous COSMOS catalogs, we present a new ultra-deep COSMOS catalog derived from our unique JWST NIRCam data. Sources are detected on the four NIRCam bands using the dual-extraction “hot and cold” detection method (Rix et al. 2004; Leauthaud et al. 2007) that achieves high completeness and purity. Measurements are performed using aperture photometry for the space-based HST/ACS F814W, JWST/NIRCam, and JWST/MIRI bands and a novel model-fitting approach based on SourceXtractor++ (Bertin et al. 2020; Kümmel et al. 2020; Kümmel et al. 2022) to measure total photometry in 37 bands as well as NIR morphology for over 700 000 galaxies. We also provide additional morphological measurements from bulge-disk decomposition, from the independent model-fitting code GaLight (Ding et al. 2020; Birrer et al. 2021), as well as from machine learning-based morphological classification (Huertas-Company et al. 2024). We used the rich multiband coverage at ~0.3–8 μm to derive photometric redshifts and physical parameters from SED fitting using LePHARE (Arnouts et al. 2002; Ilbert et al. 2006), including novel non-parametric star formation histories (SFHs) using CIGALE (Boquien et al. 2019).

This paper is organized as follows. In Sect. 2 we describe the data. In Sect. 3 we describe the methods used for source detection and photometric and morphological measurements. In Sects. 4 and 5 we describe, compare, and validate these measurements. In Sect. 6 we present the photometric redshifts derived from SED fitting with LePHARE. Finally, in Sect. 7 we present the physical properties measured from both LePHARE and CIGALE. Our conclusions are presented in Sect. 8.

We adopt a standard ΛCDM cosmology with H₀ = 70 km s⁻¹ Mpc⁻¹, Ω_m,0 = 0.3, and Ω_Λ,0 = 0.7. All magnitudes are expressed in the AB system (Oke 1974), for which a flux f_ν in μJy (10⁻²⁹ erg cm⁻² s⁻¹ Hz⁻¹) corresponds to AB_ν = 23.9–2.5 log₁₀(f_ν/μJy).

Fig. 1 Survey footprint in the COSMOS field. The regions show the footprints of different imaging instruments color coded accordingly. The JWST footprint from NIRCam and MIRI is shown in the colored areas. The gray regions mark the area affected by HSC bright star masks.

2 Observations and data reduction

In this section we provide a brief overview of the imaging data we used to construct the galaxy catalog. Figure 1 shows the footprint of the different surveys in the COSMOS field. The COSMOS-Web catalog is built from the JWST NIRCam imaging whose footprint is shown in red, while the MIRI imaging is marked in yellow. The footprint of legacy imaging surveys included in the catalog are shown with the regions that fully encompass COSMOS-Web. In Fig. 2 we show the depths of all bands as a function of wavelength and in Table 1 we provide information about the bands, including their 5σ empty aperture depths.

Fig. 2 Depths in all bands corresponding to 5σ. These are measured from the variance of the background measured in empty apertures averaged over the field of view. These apertures are 0.15″ in diameter in the NIRCam and HST bands, 0.5″ in MIRI, and 1″ for the ground-based data. No aperture correction is applied. The segment length corresponds to the filter width.

2.1 JWST data

COSMOS-Web is a 255 h JWST Cycle 1 Treasury Program (#1727) covering 0.54 deg² of contiguous four-filter NIRCam imaging (F115W, F150W, F277W, F444W) and 0.2 deg² of MIRI imaging (F770W) acquired in parallel (Casey et al. 2023). The mosaic is roughly square, measuring 46° × 46°, and is divided into 152 visits that use the 4TIGHT dither. The NIRCam mosaic has different number of exposures (one to four), and therefore different depths, at different positions, while for MIRI this ranges from two to eight. These are detailed in Tables 1 and 2 of Casey et al. (2023). In Table 1 of our paper we show the 5σ empty aperture depths averaged over the entire mosaic, while they vary on the mosaic based on the number of exposure, for example ranging from 27.74 to 28.34 mag AB for F 277 W, as detailed in Table 4 in the data reduction paper Franco et al. (2025, see Table 1 in Harish et al. 2025 for the MIRI data). This is also discussed in Sect. 3.4.2.

Data reduction of the NIRCam and MIRI data are described in detail in two forthcoming papers (Franco et al. 2025; Harish et al. 2025), which we summarise briefly. The COSMOS-Web NIRCam observations were processed using the JWST Calibration Pipeline (Bushouse et al. 2023), with additional optimizations for image quality and astrometric precision. Raw NIRCam exposures were retrieved from the Mikulski Archive for Space Telescopes (MAST) and processed with pipeline version 1.14.0, supplemented by custom corrections appropriate for deep JWST imaging (e.g., Bagley et al. 2024). These corrections included mitigation of 1/f noise (Schlawin et al. 2020), background subtraction, artifact removal (including wisp correction and claw removal), and identification and masking of defective pixels. Calibration was performed using the Calibration Reference Data System (CRDS) pmap-1223, which corresponds to the NIRCam instrument mapping imap-0285. The final science mosaics were produced at a pixel scale of 0.03″ pixel⁻¹, ensuring optimal spatial resolution for accurate photometric measurements. Astrometric refinement was conducted using the JWST/HST Alignment Tool (JHAT; Rest et al. 2023), aligning the NIRCam images to a reference catalog constructed from HST/ACS F814W mosaics (Koekemoer et al. 2007), with astrometry calibration tied to Gaia Early Data Release 3 (EDR3; Gaia Collaboration 2023). This procedure resulted in a median absolute positional offset of <5 mas, with a median absolute deviation (MAD) of <12 mas across all filters. Similar to the NIRCam data reduction, the MIRI observations were processed using version 1.12.5 of the JWST Calibration Pipeline (Bushouse et al. 2023), with an additional custom background subtraction step to mitigate the strong sky and thermal background present in our data (see also Yang et al. 2023; Pérez-González et al. 2024). Calibration was performed using CRDS context version pmap-1130, and the astrometric alignment was based on the reference catalog from the HST/ACS F814W mosaics (Koekemoer et al. 2007). The overall positional accuracy in RA and Dec was 0.35 and 6 mas, respectively, with a MAD of 28 mas.

2.2 Ultraviolet data

We use the same U– band imaging from the Canada-France-Hawaii Telescope’s (CFHT) MegaCam instrument as in COSMOS2020. These include data for the CFHT Large Area U– band Deep Survey (CLAUDS; Sawicki et al. 2019) and for COSMOS. Processing of these images is described in Weaver et al. (2022) and Sawicki et al. (2019).

2.3 Ground-based optical data

In common with previous COSMOS catalogs, most optical data comes from the Subaru telescope, using either the Hyper Suprime-Cam (HSC) or Suprime-Cam instruments (Miyazaki et al. 2002, 2018). We use the third public data release (PDR3) of the HSC Subaru Strategic Program (HSC-SSP) comprising the g, r, i, z, y broad, as well as NB0816, NB0921, NB1010, narrow bands (Aihara et al. 2022). Compared to the PDR2 used in COSMOS2020, PRD3 is more uniform and slightly deeper with improved sky subtraction and better photometric and astrometric calibrations (Aihara et al. 2022). We used the public PDR3 data access tools¹ to retrieve the processed images and weight maps from the “UltraDeep” layer. Since we use only the ~0.54 deg² central area, the depth is uniform over the COSMOS-Web footprint.

We also include medium and narrow bands from Suprime-Cam used in COSMOS2015 (Laigle et al. 2016) and COSMOS2020 (Taniguchi et al. 2007, 2015). These include 11 medium bands (IB427, IA484, IB505, IA527, IB574, IA624, IA679, IB709, IA738, IA767, IB827), and two narrow bands (NB711, NB816). These are the same images as in COSMOS2020, which are also PSF homogenized on a tile-level on the individual images (see Sect. 3.1.2 in Weaver et al. 2022).

2.4 HST optical data

We use the full COSMOS HST ACS F814W dataset (Koekemoer et al. 2007). Specifically we use HST ACS F814W mosaics produced from data that was recalibrated using updated reference files and calibration pipelines, including improved treatments for charge-transfer efficiency, updated flat fields, biases, and dark current reference files (appropriate for the observation dates corresponding to the data). These mosaics have also had their astrometric alignment updated, with their absolute astrometry directly aligned to the Gaia-DR3 reference frame² where these improvements follow the approaches developed and described in detail by Koekemoer et al. (2011), updated as needed to accommodate current calibrations.

2.5 Ground-based near-infrared data

UltraVISTA (McCracken et al. 2012) was a large public survey carried out on the VISTA (Sutherland et al. 2015) telescope with the VIRCAM (Dalton et al. 2006) instrument. The survey took place between 2010–2023. We use the final UltraVISTA DR6 which comprises the entire UltraVISTA public survey dataset as well as guaranteed time observations which were taken in 2010. DR6 comprises almost 100 000 images taken over 2290 hours of observing time. Full details are given in the documentation³. Compared to previous releases, depths are now completely uniform over the 1.5° × 1.2° UltraVISTA survey area.

Fig. 3 Summary of the key steps in preparing the photometric catalog. First, source lists are generated using the “hot and cold” detection technique and used to perform aperture photometry on PSF-homogenized space-based images. Next, sources are grouped and a joint two-dimensional model light profile is extracted for all sources. Using the morphological parameters extracted from these fits, forced photometry is then performed on all bands. Finally, in the post-processing step, the photometric errors for each source are adjusted based on variance measurements in empty apertures and the source lists are cleaned and flagged. The resulting data products include morphological measurements, model-based photometry in 37 bands, and simple aperture photometry in 6 space-based bands, for more than 700 000 galaxies.

2.6 Other data

We include data from the Spitzer telescope’s mid-infrared (Werner et al. 2004) IRAC camera in channels 1, 2, 3, and 4. These data were reprocessed as part of the Cosmic Dawn Survey (DAWN Euclid Collaboration: McPartland et al. 2025; Euclid Collaboration: Moneti et al. 2022). Although photometry is extracted in these bands, this information is not used for photometric redshift and physical parameter estimation from SED fitting.

3 Object detection and photometry

To construct our catalog, we develop a set of techniques that make optimal use the rich multiwavelength data in COSMOS. Two sets of photometric measurements are made, the first using apertures on PSF-homogenized JWST and HST images. The second uses a two-dimensional light profile model fitted on all available data, both ground and space. In both cases, the photometric extraction is carried out on sources detected using a “hot and cold” detection technique described in Sect. 3.4. Our approach is summarized in Fig. 3.

3.1 Point-spread function reconstruction and homogenization

Accurately measuring the PSF is crucial for flux measurements, profile fitting and morphological measurements. In the following, we describe the our methodology in PSF reconstruction and image homogenization.

3.1.1 Point-spread function reconstruction

To reconstruct the PSF, we use PSFEx⁴ (Bertin 2011). PSFEx can empirically reconstruct the PSF at different positions on the detector focal plane and at various pixel scales, determining how much sky area is encoded in each pixel. Berman et al. (2024) show that PSFEx provides the best performance in reconstructing the NIRCam PSF in combined images. In contrast to empirical PSF fitters like PSFEx, forward modeling approaches such as WebbPSF (Perrin et al. 2014) produce PSF models for each single exposure. However, applying these forward modeling approaches to mosaics is non-trivial (Bushouse et al. 2023; Harvey & Massey 2024).

We first run SExtractor on all bands to construct catalogs containing vignettes of all stars. We use the MU_MAX and MAG_AUTO diagram (Leauthaud et al. 2007) to select point sources. We account for the spatial variation of the PSF in two ways. First, we construct a separate model for each of the 20 tiles of the full mosaic (cf. Franco et al. in prep.) and all bands. Second, for each tile, we use order–1 polynomial interpolation to fit spatial variations of the PSF across the field, as implemented in PSFEx. Finally, we used a sampling step of 1/4th of the PSF FWHM for each band.

For validation and quality control, we use the mean relative error diagnostics developed in Berman et al. (2024) to quantify pixel-level mismatches between sources and PSF models. Finally, we examine residuals in the second adaptive moments (Hirata & Seljak 2003; Mandelbaum et al. 2005) of input sources and PSF models. This is done using Galsim’s HSM module (Rowe et al. 2015). Specifically, the FindAdaptiveMom’s function returns the size (σ_HSM) and shape (g₁, g₂) moments of an object. In this case, the objects are sources from a training catalog or a PSF model. While the NIRCam and MIRI PSFs are not well approximated by elliptical Gaussians, these second moment statistics are still useful for detecting significant size and shape mismatches when used alongside our other quality metrics.

Finally, we found no signs of systematic modeling errors using the mean relative error diagnostics. Examples of these figures for the NIRCam PSFs can be found in Appendix A. These figures show that our PSF models adequately capture the key features in the sources we are trying to model. This is supported by the strong agreement between the second moments in sources and PSF models. The average size error was less than 10% across the NIRCam filters, which we use to model each source (c.f. Sect. 3.5.4). The errors in shape were similarly low.

3.1.2 Point-spread function homogenization

The PSF FWHMs vary significantly between different bands. To ensure consistency in our aperture photometry measurements in different bands, we construct a set of PSF-homogenized mosaics for the JWST/NIRCam and HST data. We adopt the lowest resolution NIRCam F444W PSF as our target. We followed a procedure similar to Weaver et al. (2022), using the Python tool pypher (Boucaud et al. 2016) to generate convolution kernels with a band-dependent regularization factor. This parameter was optimized to minimize residuals between the convolved PSF model and the target while preventing harmonic artifacts from Fourier transformations. All JWST/NIRCam and HST tile mosaics were convolved with their respective kernels, ensuring a uniform PSF across all images, and accounting for spatial variations between tiles. We also generate a convolution kernel to match the NIRCam F444W image to the MIRI F770W PSF, which we use to examine the fraction of flux lost when performing aperture photometry on the MIRI data, as described in more detail in Sect. 3.4.3. While convolution can introduce correlated noise in the backgrounds, potentially affecting flux measurements, this effect is accounted for in the photometric error estimates (see Sect. 3.6).

3.2 Masking

Parts of our mosaics contain bright low-z galaxies, diffuse scattered light from bright stars, and imaging artifacts. Here, photometric measurements cannot be made accurately, and these regions must be identified and masked. For bright stars, we developed a semiautomatic procedure to identify and mask regions of the NIRCam mosaics affected by diffraction spikes. This requires a sufficiently large NIRCam PSF model. Since WebbPSF is not optimized beyond 30″ and does not account for scattered light, we use a larger PSF model provided by STScI⁵. This model is convolved with a Gaussian kernel and normalized to one. It is then converted into a series of region files, capturing the extent of the PSF model down to a range of limiting thresholds. This allows us to adopt the most appropriate mask for each star, accounting for their brightness. For each star brighter than mag_F150W ≤ 17.0, we automatically select the most appropriate region based on the segmentation map of the star, translate the region to the star position, and rotate to match the position angle of the observations. From there, we manually adjust the regions around each star, scaling or rotating if necessary to ensure that the entirety of the emission is masked. Imaging artifacts and a few nearby diffuse galaxies are manually added to the mask. Figure 4 shows a cutout of the detection image in the B10 tile, with the overlaid star mask.

Fig. 4 Cutout of the NIRCam ${\chi }_{+}^2$ detection image in tile B10, with the custom star mask overlaid. The mask covers bright stars and their large diffraction spikes. At the top right is a large “Dragon’s Breath” artifact (type II), which arises from the bright (K ~ 6.5 Vega mag) star to the north. This feature, and other artifacts, are manually masked.

3.3 Detection image

We construct a detection image as the coaddition of all four COSMOS-Web NIRCam bands following the “χ²” technique (Szalay et al. 1999; Drlica-Wagner et al. 2018). This method optimally combines multiple images to maximize the detection signal-to-noise. In the first step, we generate PSF-homogenized science mosaics by matching the PSF to that of the lowest-resolution F444W band (see Sect. 3.1.2). Then, for each band and tile, we construct noise-equalized images (NSCI) by multiplying the PSF-homogenized science mosaics (SCI) by the square root of the weight map (WHT)⁶, i.e., $$\mathrm{N}\mathrm{S}\mathrm{C}\mathrm{I}=\mathrm{S}\mathrm{C}\mathrm{I}\sqrt{\mathrm{W}\mathrm{H}\mathrm{T}}.$$(1)

Next, we measure the root-mean-square (RMS) of the noise-equalized images by iteratively fitting the negative tail of the pixel distribution (i.e., below 1σ) with a Gaussian function. By dividing the PSF-homogenized, noise-equalized images by their measured RMS values, we produce signal-to-noise ratio (S/N) maps for each filter. A nominal χ² detection image would be created by adding these maps in quadrature (e.g., Szalay et al. 1999). However, marginally negative S/N pixels across multiple individual bands may lead to false detections in a χ² image. To mitigate this, we modify slightly the construction of the χ² image by truncating the individual S/N maps (setting negative pixels to zero) before combination in the χ² image. Finally, we use the square root of the χ² image as our detection image. This ensures that our first-pass measurement of galaxy shapes are calculated on a detection image with pixel values proportional to flux density and not the square of the flux density.

Figure 5 illustrates the variations in the χ² distribution for this data set. Without truncation, the pixel values in the detection image follow a χ² distribution with four degrees of freedom for the four NIRCam bands. With truncation, the distribution instead resembles a χ² distribution with two degrees of freedom because the S/N maps halve the number of contributing bands. The effective 3σ, 4σ, and 5σ detection thresholds for the χ² distribution with N = 2 are indicated. We adopt an initial threshold of χ² = 13.3, corresponding to 3σ; the equivalent $\sqrt{{\chi }_{+}^2}$ value is 3.69. Sources whose peak $\sqrt{{\chi }_{+}^2}$ equals this threshold are excluded from the catalog unless they also satisfy the additional hot+cold aperture criteria described in Sect. 3.4. Under this criterion, the limiting S/N in F277W (the deepest NIRCam filter) is approximately 4σ, barring detection in any other NIRCam filter.

Fig. 5 Distribution of pixels in the PSF-homogenized, noise-equalized NIRCam maps (top panels) and χ2 distribution for the quadrature sum of the four NIRCam S/N maps (bottom panel). We show the distributions without truncation (black), with positive truncation ${\chi }_{+}^2$ (blue), and with negative truncation ${\chi }_{-}^2$ (red). The χ2 distribution for degrees of freedom N = 4 and N = 2 are shown in dark and light dashed gray lines. The distribution of the sources’ maximum ${\chi }_{+}^2$ values that are recovered by the hot+cold catalog are shown in the green histogram. Vertical lines mark the equivalent of 3, 4, and 5σ thresholds for χ2 distributions with N = 2.

3.4 The hot+cold aperture catalog

Source detection for the COSMOS-Web catalog is performed using a hot+cold scheme, which we outline in the following section. We additionally produce an aperture photometry catalog on PSF-homogenized space-based images. As we describe below, this detection methodology allowed us greater control and optimal performance in terms of maximizing completeness and minimizing contamination.

3.4.1 Source detection

Source detection and aperture photometry is performed using SEP (Barbary 2016), a python implementation of SExtractor (Bertin & Arnouts 1996). We employ a hot and cold detection scheme, similar to catalogs produced for the CANDELS survey (e.g., Galametz et al. 2013; Guo et al. 2013; Nayyeri et al. 2017; Stefanon et al. 2017) and originally presented by Rix et al. (2004). The premise of the hot and cold detection strategy is that the optimal detection parameters are fundamentally different for ultra-deep, pencil-beam surveys and shallower, wide-field surveys. While deep and narrow surveys can use aggressive detection parameters, pushing the catalog to the detection limit, shallow and wide surveys generally must be less aggressive to properly deblend bright sources. The hot and cold strategy involves first extracting objects with a high threshold optimized for bright, extended sources (the “cold mode”) and then running a separate extraction optimized for faint, isolated sources (the “hot mode”).

For the cold mode, we adopt a detection threshold of $\sqrt{{\chi }_{+}^2}$ = 4.66, which corresponds to 4σ (N = 2). The detection image is convolved with a top-hat filter kernel with a diameter of 9 pixels, to optimize detection for extended sources. We use minarea = 15, deblend_nthresh = 64, deblend_cont = 0.001, and clean_param = 2.0. For the hot-mode, we use a detection threshold of 3.698, which corresponds to 3σ (N = 2). We use a Gaussian filter kernel with a FWHM of 3 pixels, and adopt minarea = 8, deblend_nthresh = 32, deblend_cont = 0.01, and clean_param = 0.5.

The two catalogs are then merged as follows. First, bright stars are removed from the cold-mode catalog via the star masks described in Sect. 3.2. In particular, any source with >80% of its segmentation map pixels overlapping the star mask is removed from the catalog; this mostly applies to the stars themselves or galaxies falling entirely on the diffraction spikes. Sources with between 0 and 80% of their pixels overlapping the mask are not removed from the catalog, but are flagged with the keyword flag_star. Next, an elliptical mask is defined based on the cold-mode sources from a_image and b_image, using a scale factor of 6 and a minimum radius of ten pixels. Flagged sources are not used to construct the elliptical mask, since their shape parameters can be unreliable. Any hot mode source that overlaps the elliptical mask, or intersects with the cold-mode segmentation map, is removed; otherwise, it is added to the final catalog. The result is a single catalog in which the majority of sources are taken from the cold mode, except isolated/faint sources detected in the hot mode (Fig. 6).

The left panel of Fig. 6 provides an illustration of our source detection method, where cold mode sources are generally brighter and allowed to crowd with each other more significantly. Hot mode sources can only be detected in regions outside the cold mode elliptical mask. In general, hot mode sources are on the margins of the detection limits, fainter than 28.0 (27.6) mag in four-exposure (two-exposure) depth areas of the mosaics.

Fig. 6 Illustration of our “hot and cold” detection technique. Left: a 90″ × 90″ cutout of the COSMOS-Web ${\chi }_{+}^2$ detection image, showing Kron ellipses of all detected sources found from cold mode (blue) and hot mode (orange) SEP runs. Note the star in the lower-right corner that has been manually masked. Right: catalog completeness and contamination measured using an independent catalog from PRIMER (which reaches one magnitude deeper). Completeness is measured as a fraction of PRIMER sources recovered as a function of F277W magnitude (the deepest COSMOS-Web filter) for cold-mode (blue) and hot-mode (orange) sources. Contaminants are measured as sources found in COSMOS-Web but not detected in PRIMER imaging. We split these measurements into “2 exposure” (bottom row) and “4 exposure” (top row) areas where the imaging depth varies significantly. Completeness at the 8σ (5σ) survey depth exceeds 95% (80%) and contamination at these depths does not exceed 2% (4%).

3.4.2 Detection completeness and contamination

The right panel of Fig. 6 shows the detection completeness and contamination as a function of magnitude in the deepest F277W band, split between the hot- and cold-mode detections. This is measured by cross-matching with an independent reduction of the PRIMER-COSMOS Survey (Dunlop et al. 2021, GO#1837) within our larger COSMOS-Web mosaic. Processed using the same imaging pipeline, none of the COSMOS-Web imaging was used in the construction of the PRIMER mosaics or vice versa. PRIMER-COSMOS has eight bands of NIRCam coverage over 140 arcmin²; we generate a PRIMER-COSMOS hot and cold catalog with the same technique used for COSMOS-Web, by producing a PSF-homogenized positive-truncated $\sqrt{{\chi }_{+}^2}$ image (with 8 degrees of freedom for the eight filters). The PRIMER-COSMOS catalog is generally 1–1.2 magnitudes deeper than COSMOS-Web in matched filters, so our completeness and contamination measurements here are not particularly sensitive to variations in depth within the PRIMER-COSMOS mosaics themselves.

Completeness is measured as a fraction of PRIMER sources recovered by F277W magnitude bin in overlapping regions not masked in either PRIMER or COSMOS-Web. F277W magnitudes are measured in a 0.3″ diameter circular aperture without aperture correction. The blue and orange histograms indicate the relative fraction of those sources detected in the cold versus hot mode. Completeness is measured independently in regions of the COSMOS-Web mosaics at four-exposure depth versus two-exposure depth and are remarkably consistent with expectation as a function of source S/N. Contamination is measured as fraction of COSMOS-Web sources that are not recovered in the PRIMER-COSMOS catalog above a nominal 4σ detection threshold; such contaminants are deemed to be pure instrumental noise but are exceedingly rare above the nominal COSMOS-Web detection limits.

Of 784 016 sources in the COSMOS-Web catalog, 566 521 are recovered by cold mode detection (72.3%) and 217 495 are added by hot mode detection (27.7%). Of the 580 496 sources brighter than F277W < 28.2 (in $0_{\cdot }^{\mathrm{\prime }\mathrm{\prime }}3$ diameter circular apertures), 533 861 are cold mode (92.0%) and 46 635 are hot mode (8.7%). While hot mode sources dominate the tail of very faint sources, it is clear that a significant fraction of real sources can only be recovered via this extra detection mode, even well above the nominal 5σ depth of our imaging.

3.4.3 Aperture photometry

Aperture photometry is measured for HST/ACS F814W, JWST/NIRCam F115W, F150W, F277W, F444W, and JWST/MIRI F770W. Only these bands are processed for the initial hot+cold aperture photometry catalog due to their high spatial resolution (FWHM≤0.3″). All bands except F770W are PSF-homogenized to F444W, as described in Sect. 3.1.2. We provide photometry in circular apertures (aper) and Kron elliptical apertures (auto). The former is provided for convenience and simple measurements; in most cases, the auto photometry is the more reliable measure of the total flux, except for highly blended sources.

Aperture photometry is measured in circular apertures with diameters of 0.2″, 0.3″, 0.5″, 0.75″, and 1″. No aperture or PSF corrections are applied. The auto photometry is measured in elliptical apertures computed using a Kron (Kron 1980) factor of k = 1.6 and a minimum circular radius of 1.1 pixels (kron1). These small Kron apertures are intended to capture realistic colors while maximizing S/N. We corrected measurements in these small apertures to the total flux by additionally performing photometry on the F444W image using the default Kron factor k = 2.5 (kron2). This aperture correction is applied multiplicatively to the fluxes and uncertainties for all filters. We then applied an additional correction for the fraction of PSF flux that falls outside this larger Kron aperture. For this we computed a grid of PSF corrections for a range of aperture semimajor axes and axis ratios using the F444W PSF; these are applied to the fluxes based on the Kron aperture parameters. These residual PSF corrections typically range from 1–1.2. The MIRI/F770W fluxes are handled following Finkelstein et al. (2024): we perform an additional Kron aperture measurement on a F444W image which has been PSF-matched to F770W, and compare the result to our nominal F444W photometry. We correct our F770W measurements using this ratio, which accounts for the fraction of flux lost due to the larger PSF, assuming the F444W image traces the profile of the object well.

3.5 The SE++ catalog

Photometric measurements are made in 37 bands, summarized in Table 1 and Fig. 2, using the nominal images (non-PSF-homogenized). There are significant band-to-band PSF FWHM variations between ground- and space-based images in the COSMOS-Web dataset, ranging from ~0.04 arcsec to ~1 arcsec. For this reason, we adopt a multiband model-fitting approach using SourceXtractor++ where parametric models convolved with the corresponding-band PSF (described in Sect. 3.1.1) are fitted to all detected sources in all available bands.

3.5.1 SourceXtractor++

SourceXtractor++ (Bertin et al. 2020; Kümmel et al. 2020; Kümmel et al. 2022, hereafter SE++) is a multiband, multiobject model-fitting engine developed for the Euclid mission (Laureijs et al. 2011; Euclid Collaboration: Mellier et al. 2024) and a successor to SExtractor2.0 (Bertin & Arnouts 1996). The photometric and morphological parameter recovery was tested extensively in the Euclid Morphology Challenge where it achieved the highest scores (Euclid Collaboration: Bretonniere et al. 2023; Euclid Collaboration: Merlin et al. 2023). SE++ is an optimized photometric measurement program, including flexible profile fitting (combining Sérsic, constant, and point-source models) with simultaneous coupled-parameter multiband fitting. Additionally, it directly uses the native WCS information, avoiding the need for image resampling. This makes SE++ highly efficient for accurate multiband photometric and morphological measurements. Importantly, since SE++ fits a source model convolved with the corresponding PSF in each band, we use the native images that are not PSF homogenized in contrast with the hot+cold catalog. Our catalog preparation technique is described in the following sections.

3.5.2 Grouping and iterative fitting

In SE++ neighboring sources are grouped together and fitted simultaneously. This is essential to accurately retrieve the flux and shape of the sources without contamination from neighboring sources. Additionally, incorporating lower-resolution ground-based images – where light profiles from neighboring sources overlap significantly – requires more conservative, larger groups than would be required than if only NIRCam images were included. We build the groups in two steps. First, we draw ellipses around all sources by taking the positions, Kron radius, axis ratio and position angle from the input hot+cold catalog. We scale the Kron radius (measured on the NIRCam images) with a factor that corresponds to the ratio of the NIRCam and UltraVISTA PSF FWHM, and clip it between values of 0.7 and 3.0 arcsec. Next, we construct a 0.3σ threshold segmentation map (masking the stars) on the UltraVISTA K_S band. This is necessary to account for the overlapping wings of relatively bright nearby sources (see Fig. 7) and ensure that these are grouped and fitted together. We chose the UltraVISTA K_S band because it has one of the highest FWHM and wavelength coverage that sits between the NIRCam bands. Finally, we combine the images from the two steps to create the final partition map, where groups with connecting pixels are given a unique group_id. Sources from the input hot+cold catalog whose centroids fall in the same group are assigned the same group_id.

Figure 7 illustrates this technique, showing a relatively crowded region in NIRCam F277W, UltraVISTA K_S and the resulting partition image. The red contours encircle the sources that are grouped together, while the partition image shows the groups color-coded by their group_id.

Another feature of SE++ is the iterative fitting of grouped sources. Because these are processed simultaneously, many parameters must be fitted. Iterative fitting deals with this by fitting the brightest source⁷ first, while masking the others. It then subtracts 100%⁸ of the estimated flux of the brightest object, and repeats the procedure for the following brightest source. One such iteration over all sources in the group is called a “meta iteration”. We allow three meta iterations⁹ with a threshold of 10% relative change in the χ² to consider the meta iterations converged¹⁰. This procedure deals effectively with groups of blended sources, increasing the reliability of the fitted models and the derived photometry and morphology.

3.5.3 Two-stage association mode run

We used the “association” or “assoc” mode in SE++ v0. 22 to extract the photometric catalogs using a measurement catalog derived from our hot and cold detection scheme (Sect. 3.4.1). Specifically, we used the no-detection functionality of the assoc mode available in SE++ v≥0.21, to bypass the need to provide a detection image. The input assoc file provides initial values for the centroids and the model parameters (flux, half-light radius, and axis ratio) that are subsequently fitted for each source on the measurement images.

We adopt a two-stage assoc mode run where in the first, modeling stage, we fit the structural parameters (including the total flux) of the parametric models on the four NIRCam bands simultaneously. In the second, forced photometry stage, we only fit for the photometry in the remaining bands. In the following, we describe the two stages in detail.

Fig. 7 Cutout of a relatively crowded region in the high-resolution NIRCam (left), the low-resolution ground-based UltraVISTA KS-band image (middle), and the partition image (right). Groups are defined using the partition image such that sources falling in the blobs with connecting pixels (here shown in the same color) are assigned the same group_id.

3.5.4 Modeling run

In the modeling run, we fit user-defined models on every detected source in the input assoc catalog. For each source or a group of sources, the fit is done inside a square frame whose size is defined using the source Kron size with a ~25% margin. We checked that these model frames are large enough to account for the total light of the source by imposing a minimum frame size and appropriate scaling factors to the Kron size¹¹. Models are convolved with the PSF, rasterized following the pixel grid and WCS, and fitted to the science image of each corresponding band by using a modified least-squares function that aims at minimizing the residual; we used the WHT maps for this. The 1σ uncertainty estimates of the model parameters are obtained from the covariance matrix of the fit, which is computed by inverting the approximate Hessian matrix of the loss function at the best-fit values.

To enhance the scientific utility of the catalogs, we perform two independent runs of SE++ fitting all sources with the following models:

1. Sérsic models (Sérsic 1963), defined by the Sérsic index n_S, effective radius R_S,eff, axis ratio (a/b)_S, position angle θ_S, and total flux f_S,tot. We parametrize and fit the ellipticities e₁ and e₂, which are directly related to the axis ratio. We place priors on R_S,eff, n_S, and e₁, e₂, as shown in Fig. B.1;

2. Bulge + Disk models, built as composites of an exponential disk (n_S = 1) and a de Vaucouleurs bulge (n_S = 4) light profiles. We parametrize them by the effective radii of the bulge, R_B,eff, and disk, R_D,eff; the axis ratios, (a/b)_B and (a/b)_D; a common position angle θ_BD, shared by both components, which we assume to be coaxial; the total flux of both components, f_BD,tot; and the bulge-to-total ratio, B/T = f_B,tot/f_BD,tot. From this, we also derive the total fluxes of each component, f_B,tot and f_D,tot. We apply priors to R_B,eff, R_D,eff, (a/b)_B, and (a/b)_D as shown in Fig. B.1. The B/T prior follows a bell curve ranging from 5 × 10⁻⁵ to 1, with a mean and spread that increase with wavelength.

For each source, the same set of structural parameters (R_eff, n, a/b, θ) is fitted on the four NIRCam modeling bands (F115W, F150W, F277W, F444W), and we do not let them vary with wavelength. This makes the resulting parameters an effective average of all fitted bands, weighted by the corresponding weight map. As such, the structural parameters correspond to the averaged morphology over the 1–5 μm wavelength range.

3.5.5 Forced-photometry run

In the second, forced-photometry stage, we run on MIRI F770W, HST/ACS F814W, and the remaining ground-based and IRAC data by fixing the structural parameters from the modeling run and fitting only for the flux and a small coordinate offset. These are initialized with the best-fit parameters from the modeling run; in the case of the flux, the initial value for all bands is a mean of the flux in the four modeling bands from the first run. The total flux is fitted for all bands independently, without imposing any prior, while the coordinate offset (±0.05 arcsec) is fitted for all bands jointly. However, we note that this coordinate offset is not enough to capture the proper motion of stars due to the relatively large difference between the observation epochs of different bands; this is out of the scientific scope of this work.

Figure 8 demonstrates the source extraction with SE++ on a randomly selected 20″ × 20″ area, by showing the science, model and residual images of five wide filters from both modeling and forced photometry run. The models accurately reproduce the flux distribution in the science images, yielding clean residual images. This validates our methodology in measuring photometry for every NIR-detected source across many ground- and space-based bands. However, some residuals remain in the HSC g because low-redshift galaxy sizes typically decrease with wavelength (e.g., Vulcani et al. 2014), whereas we fix sizes to measurements at 1 μm to 5 μm.

3.6 Calibration of uncertainties

Up to this point in our analysis, the photometric errors only accounted for photon noise from detected sources and not for Poisson noise from the background. Consequently, the uncertainties in the catalog are underestimated for faint objects. We compute correction factors for both the hot+cold catalog (based on PSF-homogenized mosaics) and for SE++ catalog (based on native-resolution images). For aperture photometry in each band, we place 200 000 random circular apertures in a range of diameters from 0.03″ to 1″ on the noise-equalized version of each mosaic, avoiding regions containing detected sources (as defined by the segmentation maps). For a fixed aperture size and filter, we empirically determine the background Poisson noise by fitting the negative tail of the distribution of measured fluxes with a Gaussian function. Following Labbé et al. (2003); Gawiser et al. (2006); Whitaker et al. (2011); Skelton et al. (2014); Rieke et al. (2023b); Finkelstein et al. (2024), we model the Poisson background noise as a function of aperture area (number of pixels, N) by fitting the relation $${\sigma }_N=\alpha N^{\beta /2},$$(2)

where σ_N is the noise in an aperture containing N pixels, and α and β are fit parameters.

Calibrating model-based photometric uncertainties is more complex. SE++ determines the best-fit photometric and structural parameters of a model by minimizing a loss function, which is a χ² difference between the model and observed data. The uncertainties on these parameters are derived using the approximate Hessian matrix of the loss function at the best-fit values, which approximates the covariance matrix of the fitted parameters. Since this naturally includes correlations between parameters, the uncertainties given by SE++ are the marginalized errors from this covariance matrix and includes both the pixel noise (from the variance map) and the morphology degeneracies. However, these uncertainties are typically underestimated by about a factor of two (Euclid Collaboration: Bretonniere et al. 2023; Euclid Collaboration: Merlin et al. 2023).

Model-based photometry is not measured in a fixed area, unlike circular (or elliptical) apertures, complicating the uncertainty calibration. However, that uncertainty should scale with the effective size over which photometry is measured, and so we first must generate an estimate of the effective aperture area for each source given its model characteristics. To do this, we generate models on a three-dimensional grid of R_eff, n_S and a/b with 20 points along each dimension. For each point on the grid and for each band, we generate a Sérsic model on the image pixel grid and convolve it with the corresponding PSF.

We then generate elliptical apertures which contain 90% of the model flux and count the effective pixel area in that elliptical aperture, N_eff. The adoption of 90% contained flux was calibrated empirically against circular sources in the field whose aperture photometry is quite similar to their model-based photometry. Then for each source in our catalog, we measure N_eff by interpolating numerically from the R_eff, n_S, a/b grid using scipy’s RegularGridInterpolator. We then adopt the same power law scaling measured using circular apertures (σ_Neff$\alpha N_{\text{eff\hspace{0.17em}}}^{(\beta /2)}$) to generate an estimate of the Poisson background noise for model-based photometric measurements.

All photometric uncertainties already present in the hot and cold and SE++ catalogs are then added in quadrature with this calculation of Poisson background noise (or random aperture noise per source and per band), σ_N. We note that generally, the Poisson background noise calculation dominates the error budget of sources close to the detection threshold ≲3σ).

Fig. 8 Demonstration of the performance of our cataloging methodology with SourceXtractor++ over a relatively crowded and randomly selected 20″ × 20″ area. The top, middle, and bottom rows show the science, model, and residual images, respectively, in five wide filters: HSC g, HST/ACS F814W, UltraVISTA KS, NIRCam F277W, and MIRI F770W. The science and model images are log scaled, while the residual is linearly scaled between ±6 standard deviations of the residual cutout image. The RGB science and model images are made using F444W, F277W, and F150W.

3.7 Cleaning and flagging

The catalog may contain spurious sources or objects with incorrect photometry. In this section we describe the procedures we adopted to identify these problematic sources and identify them with the warn_flag keyword.

3.7.1 Hot pixels

Despite extensive efforts to mitigate imaging artifacts, small clusters of one to four pixels with an unusually high S/N remain in our NIRCam mosaics. Although masking these “hot pixels” on the images would be the best approach, this task is impractical for COSMOS-Web; therefore, hot pixels remain in the final mosaics. We therefore flag hot pixels at the catalog level.

PSF homogenization (for the construction of the detection image) renders hot pixels less obviously identifiable; thus, many are included as sources in the initial hot+cold catalog. Nominally, one should make a concerted effort to remove them from hot+cold before proceeding to the next stage of catalog construction. This is because, in theory, known false sources present in the catalog could impact the model photometry of neighbor sources. However, in making this catalog, we have discovered that removing them at a later stage (after running SE++) produces nearly identical results. This is because the Sérsic model fits built on NIRCam imaging for hot pixels are intrinsically small (<1 pixel). Therefore, we proceed by flagging them only after the full catalog is made.

Hot pixels are identified using a curve of growth analysis, where concentric circular apertures measure the enclosed flux on native resolution NIRCam images. Sources for which concentric apertures indicate a morphology more compact than the intrinsic PSF in any of the four NIRCam bands are flagged as hot pixels. Specifically, we calculate the ratio of flux densities, F_ν, in circular apertures that are $0_{\cdot }^{\mathrm{\prime }\mathrm{\prime }}1$ and $0_{\cdot }^{\mathrm{\prime }\mathrm{\prime }}25$ in radius, dubbed apertures 0 and 1 respectively, i.e. R_[filt] ≡ F_[filt],0/F_[filt],1. Note this is applied to non-PSF homogenized aperture photometry, which is computed as described in Sect. 3.5.4. We then derive an error on R_[filt], σ_R[filt] using error propagation (neglecting covariance terms). We then empirically calibrate the threshold ratios for each filter using their native resolution, above which we only expect hot pixels to reside. These thresholds R_thresh are 0.75, 0.70, 0.39, and 0.36 for F115W, F150W, F277W, and F444W respectively. Sources are flagged as hot pixels if they are more than 1σ above this threshold in any of the four filters, i.e. R_{source,[filt]} − σ_R[filt] ≥ R_thresh[filt]. This curve-of-growth criterion results in 38916 sources flagged as hot pixels (of which 12563 are brighter than [F277]<28.2).

After iteratively visually inspecting sources flagged as hot pixels, it was determined that a subset were real galaxies and could be recovered by requiring additional criteria for flagging. For example, if best-fit Sérsic model radii are found to be smaller than a single NIRCam pixel ($0_{\cdot }^{\mathrm{\prime }\mathrm{\prime }}03$) then the hot pixel flag is retained; if the fit is larger, it is likely a true astrophysical source whose emission may manifest like a hot pixel in one of the four NIRCam filters, but modeled across all filters (Sect. 3.5.4), is substantially larger than a pixel. Modifying the hot pixel criteria with this size criterion reduces the number of sources flagged as hot pixels to 13 241 (of which 6940 are brighter than [F277] < 28.2).

To summarize, hot pixels are identified as sources fulfilling the curve of growth criteria for an NIRCam filter and have SE+++ Sérsic model fits with radii <$0_{\cdot }^{\mathrm{\prime }\mathrm{\prime }}03$. Hot pixels are assigned warn_flag = 1.

3.7.2 Assuring consistency between space- and ground-based bands

Inconsistent photometry between space- and ground-based bands is identified in the SE++ catalog in cases where the flux ground-based bands disagree with those from NIRCam bands at the similar wavelengths. This mainly happens for relatively faint sources near the detection limits which are blended with brighter neighboring sources in the low-resolution ground-based bands. Due to association confusion, SE++ can fit a source that has a biased flux several magnitudes brighter than what is measured at similar wavelength in the high-resolution NIRCam where there is no association confusion. The improved and conservatively large grouping described in Sect. 3.5 helps alleviate many such association confusions, but a fraction remains. We flag sources with inconsistent ground- and space-based photometry by identifying those where the lower flux uncertainty in UltraVISTA Y, J, or H bands exceeds twice that of the nearest NIRCam band, as follows: $$\begin{array}{rl}& (2\times F_{\mathrm{F}115\mathrm{W}}<(F_Y-\delta F_Y))\mathrm{\&}\delta {\mathrm{m}\mathrm{a}\mathrm{g}}_Y<0.5\\ & (2\times F_{\mathrm{F}115\mathrm{W}}<(F_J-\delta F_J))\mathrm{\&}\delta {\mathrm{m}\mathrm{a}\mathrm{g}}_J<0.5\\ & (2\times F_{\mathrm{F}150\mathrm{W}}<(F_J-\delta F_J))\mathrm{\&}\delta {\mathrm{m}\mathrm{a}\mathrm{g}}_J<0.5\\ & (2\times F_{\mathrm{F}150\mathrm{W}}<(F_H-\delta F_H))\mathrm{\&}\delta {\mathrm{m}\mathrm{a}\mathrm{g}}_H<0.5.\end{array}$$(3)

If a source satisfies these four conditions, then it is assigned warn_flag = 2. Unsurprisingly, many warn_flag = 2 are nearby bright stars that in the ground-based bands contaminate the flux of many neighboring sources. Outside HSC star masks, about 1.6% of the sources have warn_flag = 2 (Table 2).

Additionally, we identify a second case where the ACS/NIRCam flux is bright enough so that the sources should be detected in the HSC and UltraVISTA bands but are not. In this case, the model magnitude is zero or negative but the diameter aperture flux in the corresponding band is positive, F_[filt](1″) > 0: $$\begin{array}{rl}& F_i\le 0\mathrm{\&}\delta F_i>0\mathrm{\&}{\mathrm{m}\mathrm{a}\mathrm{g}}_{\mathrm{F}814\mathrm{W}}<28\mathrm{\&}\delta F_{\mathrm{F}814\mathrm{W}}>0\\ & F_H\le 0\mathrm{\&}\delta F_H>0\mathrm{\&}{\mathrm{m}\mathrm{a}\mathrm{g}}_{\mathrm{F}150\mathrm{W}}<26\mathrm{\&}\delta F_{\mathrm{F}150\mathrm{W}}>0\\ & F_{K_S}\le 0\mathrm{\&}\delta F_{K_S}>0\mathrm{\&}{\mathrm{m}\mathrm{a}\mathrm{g}}_{\mathrm{F}277\mathrm{W}}<26\mathrm{\&}\delta F_{\mathrm{F}277\mathrm{W}}>0.\end{array}$$(4)

If a source satisfies these three conditions, then it is assigned warn_flag = 3. Table 2 quantifies the number of sources affected by these flags.

3.7.3 Other artifacts

We identify and flag other artifacts as follows. Sources with detection in only one NIRCam or MIRI band are predominantly “snowballs” and hot pixels and are flagged with warn_flag = 4. These have a magnitude error lower than 0.2 mag in one of the NIRCam or MIRI bands, but a magnitude error larger than 0.5 mag in all other bands (or below the limiting magnitude of the band). Sources with unrealistically small radii (R_eff < 0.00047) are typically noise detections and are flagged with warn_flag = 5. Finally, sources with an unrealistic flux ratio between $0_{\cdot }^{\mathrm{\prime }\mathrm{\prime }}25$ and $0_{\cdot }^{\mathrm{\prime }\mathrm{\prime }}1$ that is smaller than the ratio identified for stars are also identified to be snowballs and hot pixels and are flagged with warn_flag = 6. More precisely, the flux ratio condition is F_F277W($0_{\cdot }^{\mathrm{\prime }\mathrm{\prime }}25$)/F_F277W($0_{\cdot }^{\mathrm{\prime }\mathrm{\prime }}1$) < 1.6 and F_F444W($0_{\cdot }^{\mathrm{\prime }\mathrm{\prime }}25$)/F_F444W($0_{\cdot }^{\mathrm{\prime }\mathrm{\prime }}1$) < 1.9.

In summary, we provide the warn_flag to flag artifacts and sources with potentially problematic photometry and derived photo-z. For the most secure sources that do not satisfy any of the criteria described above, we assign warn_flag = 0, and we advise users to use this sample for most scientific applications. Sources with warn_flag = 2, 3 are to be handled with care; they have NIRCam and ACS photometry that is not problematic, but the issues with the ground-based photometry prevents robust SED fitting. We advise users to discard sources with all the other flags warn_flag = 1, 4, 5, 6 from scientific analysis, or to carefully inspect them.

Fig. 9 Number counts for the SE++catalog. Left: counts in NIRCam and MIRI bands. Right: counts in HST/ACS and for ground-based broadband filters. In both cases, counts are computed in 0.35 magnitude bins and normalized by the effective area after masking, namely 0.43 deg2 in all bands except for MIRI, which covers 0.15 deg2.

3.7.4 HSC star mask

In our ground-based imaging, and in particular the HSC data, the larger PSF of bright stars compared to JWST data means more sources are affected by photometric contamination. We flag these sources using the HSC star masks (Coupon et al. 2018) from the COMOS2020 catalog (Weaver et al. 2022). These are conservative and flag all sources with flux contamination from stars in all ground-based bands. Table 2 shows the number of sources remaining in the catalog after applying the HSC star mask flag. We flag sources whose ground-based photometry is affected by bright stars with flag_star_hsc = 1. Compared to the total number, 17% of the sources are affected by the HSC star mask and therefore have unreliable ground-based photometry. In terms of area, about ~20% is affected by the HSC star mask. The reason why 17% of the sources are flagged is because detections are not carried out at the cores of these masks, in the area affected by the JWST masks.

4 Photometric validation and comparisons

In this section we compare our two sets of photometric catalogs. For SE++ we provide total photometric quantities in all 37 bands, while for the hot and cold catalogs we provide PSF-homogenized aperture photometry along with aperture-to-total correction in the four NIRCam, one MIRI, and the HST/F814W band. For the SE++ catalog, we also compare the ground-based photometry to the COSMOS2020 catalog. We use the total photometry derived from the Sérsic model-fitting of SE++ as the primary photometric reference.

4.1 Magnitude number counts

Figure 9 shows the galaxy number counts computed using SE++ photometry in the JWSTand HST/ACS and the ground-based broad bands. Galaxies are selected using the star-galaxy classification described in Sect. 6.3.

Overall, the slope of the magnitude number counts at intermediate magnitudes shows the expected trend (Metcalfe et al. 2001) with respect to wavelength, with bluer bands having steeper slopes. As expected (Gardner et al. 1993), at longer wavelengths, the slope of the counts flattens at intermediate magnitudes. The amplitude of the break becomes more important at redder wavelengths. In a recent paper, Manzoni et al. (2025) provide an excellent summary of the physical origins of these changes.

We also compared with number counts from COSMOS2020. In the left panel of Fig. 9 we show the F277W number counts for sources in our catalog that match and are detected in COSMOS2020 that turn over at ~25.5 mag, comparable to the K_S depth of COSMOS2020. At m_F277W ≳ 25 the COSMOS2020-matched counts become shallower, which is likely due to the deeper and redder selection function of COSMOS-Web. In the right panel of Fig. 9 we directly compare the K_S number counts from COSMOS-Web and COSMOS2020 (taken from Weaver et al. 2022). This includes all sources detected in the respective catalogs. There is a relatively good agreement with two noticeable differences. First, the COSMOS-Web K_S counts turn over at a fainter magnitude ~26.5 compared to ~25.5 for COSMOS2020 due to the deeper UltraVISTA DR6 used in COSMOS-Web versus DR4 used in COSMOS2020. Second, the COSMOS-Web K_S counts are slightly steeper, likely due to the deeper and redder selection function.

4.2 Comparison between SE++ and hot and cold catalogs

We compare the SE++ model magnitudes and auto magnitude fluxes in the hot+cold catalog in Fig. 10. The auto photometry is corrected to total as described in Sect. 3.4.3. We also show the regions corresponding to the ±1 σ and ±3 σ photometric uncertainties, obtained by adding in quadrature the errors from both photometric data sets.

Overall, there is excellent agreement between the two sets of photometric measurements, with a running median offset (Δ = mag_SE++ − mag_hot+cold) < 0.13 mag and a median scatter of ~0.3 mag. This offset and scatter is larger for the NIRCam SW F115W and F150W bands and smallest for the NIRCam F277W and F444W, remaining within the ±1σ photometric uncertainty envelope. The offset shows that SE++ magnitudes are overall slightly fainter for faint and brighter for bright sources. These differences can come both from inadequate modeling by SE++ or imperfect PSF homogenization and uncertain aperture-to-total factors applied in hot+cold. Despite these secondary differences, the comparison shows that the two photometry sets are consistent within the uncertainties.

Fig. 10 Photometric comparison between SE++ model magnitudes and auto magnitudes from the hot+cold catalog, as a function of auto magnitude. The auto photometry is corrected to total as described in Sect. 3.4.3. The panels show the following bands: F115W (top left), F150W (top center), F277W (top right), F444W (bottom left), HST F814W (bottom center), and F770W (bottom right). In each panel, the solid teal line indicates the running median offset, while the envelope denotes the 3σ-clipped standard deviation. The dotted white curves show the ±1σ and ±3σ envelopes of the combined photometric uncertainties. The vertical dashed teal line marks the 5σ depth in each band. Summary statistics report the median offset, Δ, and standard deviation, σ, for sources brighter than this limit. Sources with S/N < 3 and warn_flag > 2 are removed from the comparison.

4.3 Comparison with COSMOS2020

To validate model-fitting ground-based photometry from SE++, we compare with the COSMOS2020 catalog (Weaver et al. 2022). COSMOS2020 is a galaxy catalog covering ~2 deg² in the COSMOS field and provides both aperture photometric measurements (CLASSIC catalog) and from model-fitting (The Farmer catalog). We compare sources matched in coordinates within $0_{\cdot }^{\mathrm{\prime }\mathrm{\prime }}6$. To ensure a comparison of consistent quantities, i.e., model-derived total fluxes, we compare with The Farmer photometry.

There are two noticeable differences in the methodologies of COSMOS2020/The Farmer and SE++ that need to be kept in consideration in the comparison. The first one is the different models that are fitted to sources. The SE++ catalog fits Sérsic models to all sources, where the structural parameters are fitted on the NIRCam bands simultaneously. The Farmer, on the other hand, uses a decision tree to decide which one out of five discrete models to fit. This is described in detail in Weaver et al. (2023), but briefly, these models are a pointsource, a circularly symmetric and exponential light profile at a fixed effective radius, an exponential profile (Sérsic model with n_S = 1), a de Vaucouleurs profile (Sérsic model with n_S = 4) and a composite of an exponential and a de Vaucouleurs profile. The structural parameters are fitted on a chi-squared combination of izYJHK_s bands and the photometry is measured on all bands individually in a forced photometry approach, similar to the SE++ catalog. The second difference is the input dataset. These are presented in detail in Sect. 2 of this paper and in Sect. 2 of Weaver et al. (2022), but briefly, the main difference is that we use HSC PDR3 data processed with the HSC pipeline, as opposed to HSC PDR2 and a custom pipeline in Weaver et al. (2022). Furthermore, for the UltraVISTA data we use DR6, as opposed to DR4 in COSMOS2020. Keeping these differences in data and methodology in mind, we compare the magnitudes and colors of both independent catalogs.

4.3.1 Photometric comparisons with COSMOS2020

In Fig. 11 we compare the total magnitudes from SE++ and The Farmer in all the broad, medium and narrow ground-based bands and one IRAC (ch1) band. The density histogram shows the magnitude difference between SE++ and The Farmer as a function of magnitude. We mark the median offset with a solid yellow line, and the 3σ-clipped standard deviation with an envelope; in the legend we show the averaged offset (Δ) and standard deviation (σ) for magnitudes brighter than the 5σ depth.

Overall, there is excellent agreement between the two sets of photometric measurements with averaged offset Δ ≲ 0.08 mag, and with σ ≲ 0.3 for all bands. Additionally, there is no significant magnitude trend of the rolling offset higher than 1 σ significance. There are, however, some noticeable second-order trends that are more prominent for the HSC bands. This refers to the kink around 25–26 mag toward brighter SE++ magnitudes (but offset less than σ). This could be due to several reasons, from the different image background modeling to the difference in PSF and source modeling in The Farmer, as mentioned above.

Finally, the offset and scatter is largest for the IRAC bands at fainter SE++ magnitudes. This is likely due to our simplistic treatment of the IRAC PSF where we model the high spatial variability with a second-order polynomial (as implemented by PSFEx). Further complication is the confusion by forcing models that are significantly below (~1.5 mag) the depth of the IRAC bands (compared to the NIRCam detection) coupled with the severe source blending. We do not investigate or correct this in further detail, since the IRAC bands are largely redundant, given the JWST coverage. We only include them in the photometry for legacy value and do not use them in the photo-z and physical parameter inference.

Fig. 11 Photometric comparison with the COSMOS2020 THE FARMER catalog. The density histogram shows the difference between SE++ and THE FARMER as a function of THE FARMER magnitude. The solid yellow line shows the median offset, while the envelope marks the 3σ-clipped standard deviation. The dashed line shows the median offset computed down to the 5σ depth magnitude, shown by the vertical, dashed teal line. The legend shows the average offset and standard deviation for magnitudes brighter than the 5σ depth.

Fig. 12 Color comparison in the ground-based bands with the COSMOS2020 Farmer catalog. Similarly to Fig. 11, the density histogram shows the color difference between SE++ and the FARMER as a function of magnitude. The solid yellow lines and envelopes show the median offset and the 3σ-clipped standard deviation. The legend shows the averaged offset (also marked in a dashed yellow line) and the standard deviation for magnitudes brighter than the 5σ depth.

4.3.2 Color comparisons with COSMOS2020

Because colors provide the most stringent constraints in SED fitting, validating them is crucial. We compare colors from SE++ with The Farmer photometry in COSMOS2020 for u − g, g − r, r − i, i − z, z − J, and J − K_S bands in Fig. 12. The density histogram shows the color difference as a function of magnitude, and the solid yellow line and envelope show the rolling median offset and 3σ-clipped standard deviation.

Overall, there is excellent agreement in the colors, with an averaged color offset less than 0.08 mag and a standard deviation of <0.4. The difference is highest for the r − i color, where SE++ shows redder colors with about 1σ significance for the faint end. In the case of z − J, there is also a 1σ significant offset toward bluer SE++ colors for bright sources. Additionally, there are small second-order curvatures for the faintest magnitudes, more notably for the HSC bands, similar to the magnitude comparison, and due to the same origins. This can be due to slight differences between the HSC DR2 versus DR3 data and photometric calibrations. Nonetheless, the overall excellent agreement in both photometry and colors means that the SE++ photometry passes this test.

5 Morphological measurements

A major advantage of model fitting is that we obtain morphological measurements for all detected sources. Furthermore, to maximize the scientific applications of the catalog, we provide independent measurements from Galight and morphological classification from machine learning. We present and compare these for sources with S/N > 5 and warn_flag = 0, 2, 3 in the following.

5.1 SE++

As described in Sect. 3.5, we measure morphologies from the two independent SE++ fits: (1) Sérsic models, resulting in the Sérsic index n_S, total effective radius R_S,eff and axis ratio (a/b)_S; and (2) Bulge + Disk models, measuring bulge and disk effective radii R_B,eff, R_D,eff, axis ratios (a/b)_B,eff, (a/b)_D,eff and bulge-to-total flux ratio B/T. These structural parameters are effectively averaged over the 1–5 μm wavelength range, with B/T being measured for all bands. Figure 13 shows the distributions of several structural parameters: R_S,eff vs. F277W magnitude in the top left, R_S,eff histogram for the Sérsic, bulge and disk components in the bottom left, n_S histogram in the bottom middle and (a/b)_S in the bottom right panel.

The R_S,eff distribution peaks around 0.1″, and shows a clearly isolated locus of unresolved, point-like sources below 0.01″. Unresolved sources are expected to have their measured size independent of magnitude; however, the slight slope shown in Fig. 13 is due to the uniform prior on R_S,eff with a lower limit of 0.01× the Kron radius from hot+cold. There is also a cloud of ~2% of sources with sizes larger than the main locus, a lot of which are also scattered toward large R_S,eff compared to the independent measurements from Galight (Sect. 5.2). These are predominantly faint (mag ≳ 27) in hot+cold in many cases fainter than the depth of one or more of the four modeling bands, and may have brighter companions nearby. This can drive a large solution for the model R_S,eff.

The Sérsic index distribution peaks at n_S ~ 1 and exponentially decreases at higher values. There is an accumulation of solutions at the minimum and maximum allowed values (0.3 and 8.5), mainly because of degeneracy of the solutions for unresolved sources. The axis ratio shows a distribution peaking at (a/b)_S ~ 0.4 and skewed toward lower values. These distributions have the expected shape, agreeing with those obtained independently with Galight and serve as a validation of the morphological measurements.

Fig. 13 Distributions and comparisons of several morphological quantities measured by SE++ and Galight for sources with with S/N > 5. Top row: density histogram of the effective radius from the Sérsic fit as a function of F277W magnitude, comparison of Reff (middle) and nS (right) between SE++ and Galight. Bottom row: histograms of the Sérsic total, bulge, and disk effective radii (left), Sérsic index (middle), and axis ratio (right) from SE++ and Galight.

5.2 Galight

An independent structural measurement using Galight, which is built on Lenstronomy (Ding et al. 2020; Birrer et al. 2021), complements SE++ measurements by providing structural parameters measured in four NIRCam bands individually, rather than averaged over a broad wavelength range as in SE++. The Sérsic parameter constraints include the half-light radius between 0.01" arcsec and the cutout frame size, the Sérsic index between 0.3 and 9, and the axis ratio between 0.1 and 1, with initial guesses derived from the SE++ catalog. Galight adopts the same source detection maps and PSFs as SE++, ensuring consistency. These modeling details are described in Yang et al. (2025). Furthermore, Galight offers two additional independent morphology measurements, one using two Sérsic models and another employing a single Sérsic with a central point source model. This catalog and its release are described in Yang et al., in prep.

5.3 Machine learning

Galaxy morphological classification can be efficiently carried out by machine learning (ML) models such as supervised convolutional neural networks (CNN) by extracting image features that are correlated with galaxy morphology (Huertas-Company & Lanusse 2023, for a review). We apply the supervised ML models from Huertas-Company et al. (2024) to classify galaxies into four broad morphological classes: spheroid, disk-dominated, bulge-dominated and irregular. These are described in detail in Huertas-Company et al. (2025). Briefly, a supervised CNN model is trained with labels from the CANDELS survey and domain-adapted to JWST using adversarial domain adaptation. We used the same architecture, input image size (32 × 32 pix) and normalization as Huertas-Company et al. (2024) and train the network using COSMOS-Web images as target domain. Uncertainties in the classification are estimated with an ensemble of ten separate trainings for each of the three F150W, F277W and F444W filters, with different initial conditions and training sets. The final probability for each morphological class in each of the three filters is computed as the average of the outputs from the ten networks. For validation, Huertas-Company et al. (2024) shows that this ML model achieves about 80–90% agreement with visual classification (tested on CEERS, but with similar performance in COSMOS-Web, Huertas-Company et al. 2025).

5.4 Comparison between the methods

Comparison of the morphological measurements between the three methods is shown in Figs. 13 and 14, showing generally excellent agreement. The effective radii show a relatively tight correlation, centered at the one-to-one ratio and with no significant offsets. The striping is a result of the parameter bounds and, for Galight the fact that the stamp size scales with the R_eff solution by SE++. The large scatter in the Sérsic index highlights the difficulty of measuring this parameter, but we find no significant offset between SE++ and Galight. The parameter distributions also agree well.

In Fig. 14 we show the histograms of the Sérsic index and B/T in F277W measured from SE++ for four morphological classes from the ML technique for spheroid, disk-dominated, irregular and bulge-dominated. Galaxies are assigned to a morphological class based on the highest probability from the four categories, measured in F277W. In general, there is good agreement with the distributions of n_S and B/T and the corresponding morphological classes. For example, spheroid and bulge-dominated systems show distributions skewed toward larger n_S and B/T with medians of n_S ~ 2.5 and B/T ~ 0.55. On the other hand, irregular, and disk-dominated systems tend to have lower values with medians of n_S ~ 1.2 and B/T ~ 0.15. Additionally, as expected, irregulars tend to be more extended than disk dominated, while spheroids less compact than bulge dominated. The dispersion is, perhaps, relatively high, which can be due to both uncertain ML classifications and SE++ measurements.

Fig. 14 Distributions of the Sérsic index (top row) and B/T in F277W (bottom row) measured from SE++ for spheroid, disk-dominated, irregular, and bulge-dominated galaxy types classified with machine learning. The dashed vertical line marks the median. These include sources with S/N > 5 classified as galaxies.

6 Photometric redshifts

6.1 Method

We use the template-fitting code LePHARE (Arnouts et al. 2002; Ilbert et al. 2006) to measure photometric redshifts. While this code was used for previous COSMOS catalogs (Laigle et al. 2016; Weaver et al. 2022), we modified significantly the configuration for the COSMOS-Web catalog. The main changes are outlined below.

The photometric redshifts of galaxies in the previous versions of the COSMOS catalog were estimated with a limited set of templates from Ilbert et al. (2013) which reduces the risk of degeneracies in the color-redshift space. Given the improvement in the multicolor coverage and sensitivity from NIRCam imaging, we decided to increase the size of the galaxy template library with a more representative library generated using the Bruzual & Charlot (2003) Stellar Synthetic Population models (hereafter BC03). We select the same set of templates as the ones used to derive physical parameters in Ilbert et al. (2015). The advantage of such an approach is that physical parameters are derived simultaneously, with a library more diverse and representative than previous COSMOS catalogs. These templates are built assuming six different Star Formation Histories (SFH), including four exponentially declining and two delayed analytic forms, as well as two metallicities (solar and half-solar). To minimse computational time, we restricted the number of ages at 43, ranging from 0.05 to 13.5 Gyr. We add emission lines to the BC03 templates by adopting the recipe in Saito et al. (2020). The emission line fluxes are allowed to vary together by a factor of two around the fiducial value during the fit.

Dust attenuation is a crucial ingredient in the modeling of galaxy SEDs. Given the variety of galaxy populations present in the COSMOS-Web catalog, we decided to explore a wide range of dust attenuation applied to the templates, with a maximum E(B − V) ≤ 1. Dusty galaxies are known to be present in the field (e.g., McKinney et al. 2025) and exploring such a range of attenuation is necessary to recover their photo-z. We included three different attenuation curves (Calzetti et al. 2000; Arnouts et al. 2013; Salim et al. 2018). The high-z analog dust attenuation curve from Salim et al. (2018) includes a bump at 2175Å. We find that this bump reduces systematic differences between modeled and observed magnitudes at fixed redshift, compared to the Prevot et al. (1984) curve previously used. The BC03 templates do not include dust emission in the modeled spectra. Since the COSMOS-Web wavelength coverage probes 7.7 μm with MIRI imaging, the Polycyclic aromatic hydrocarbon (PAH) emission could already contribute to the SED of z < 0.2 galaxies. Therefore, we modeled this dust emission. We assumed that the energy absorbed in ultraviolet-to-optical range is fully remitted in IR by the dust. We rescaled the dust emission template library from Béthermin et al. (2012) based on Magdis et al. (2012) to the expected IR luminosity using energy balance. The intergalactic medium (IGM) absorption is accounted for by using the analytical correction of Madau (1995).

Before performing a fit on the full catalog, we adapted the absolute calibration of the photometry following Ilbert et al. (2006). We restricted the fit to the galaxies having a spectroscopic redshift. The median of the difference between the predicted and the observed magnitudes is used to modify the absolute calibration in each band. Then, the fit is redone iteratively after having applied the correction. The procedure converged in two iterations. The values of the offsets are shown in Fig. 15. These offsets usually remain below 0.05 mag. We note that the low values obtained in the NIRcam bands are driven by the better sensitivity in these bands relative to ground-based data, which anchors the fit at these wavelengths.

The fitting procedure outputs several quantities. For each element i in the library, we derived a probability P_i associated with the ${\chi }_i^2$. We produce a likelihood distribution which is the sum of all individual P_i at a given redshift. This distribution is considered as our posterior redshift probability distribution function (PDF_z), assuming a flat prior for all parameters. We adopt the median of the PDF_z as a fiducial redshift point-estimate. The associated uncertainties encompass 68% of the PDF_z. We also provide the redshift corresponding to the template which produces the minimal χ² over the full library. This redshift does not necessarily correspond to the median of the PDF_z in case of a multi-peaked PDF.

When running LePHARE, we fit the observations with stellar and AGN templates in parallel to the galaxy templates. For the stellar library, we adopt the same library as Kauffmann et al. (2022) which includes brown dwarf templates. This is necessary to identify degeneracies with high redshift galaxies. For the AGN library, we adopt the templates from Salvato et al. (2011) which includes hybrid between galaxies and AGN emission. However, we emphasize that a specific work needs to be carried out to measure the photometric redshifts of sources dominated by an AGN which required additional steps not performed in this paper (e.g., correction for source variability, applying prior based on the X-ray emission).

Fig. 15 Magnitude offsets as a function of band wavelength derived as the median difference between the model and observed magnitudes with LePHARE. The shaded regions mark ±0.1, ±0.05 and ±0.01 mag. These are used as an absolute calibration of the photometry when fitting with LePHARE.

6.2 Photometric redshift performance

For the absolute photometry calibration and assessing the photo-z performance, we use a sample of about 12 000 spectroscopic redshifts with a high (>97%) confidence level out to z = 8. These are compiled from most spectroscopic programs in COSMOS (both public and private) and are presented in detail in Khostovan et al. (2025).

For the performance assessment, we use the standard metrics, such as the median absolute deviation (MAD), defined as $${\sigma }_{\mathrm{M}\mathrm{A}\mathrm{D}}=1.48\times \mathrm{median}\left(\frac{|\mathrm{\Delta }z-\mathrm{median}(\mathrm{\Delta }z)|}{(1+z_{\text{spec\hspace{0.17em}}})}\right),$$(5)

where Δz = z_phot − z_spec. Sources whose photo-z solutions deviate by |Δz| > 0.15(1 + z_spec) are classified as outliers (Hildebrandt et al. 2012). Finally, the overall redshift bias is defined as the median of Δz.

We evaluate the photo-z performance for different magnitude, color, and galaxy type selected samples in Figs. 16 and 17 and Table 3. In general, the photo-z show excellent performance, with σ_MAD = 0.012, <2% outliers and b = 0.007 for galaxies brighter than the total model magnitude m_F444W = 28, the 5σ detection limit in F444W (top left panel of Fig. 16). The performance decreases as a function of magnitude, going from σ_MAD = 0.011 and η = 1.44% for m_F444W < 23, to σ_MAD = 0.030 and η = 9.52% for the faintest sample at 26 < m_F444W < 28. For the fainter samples, there is an increasing number of sources with z_phot≲ 1.5 and z_spec≳ 2, which can be attributed to a misidentification between the Lyman and Balmer breaks for sources with lower signal-to-noise ratio. The bias does not show a particular magnitude dependence. We note that typical COSMOS-Web objects are expected to be faint, e.g., 25 ≲_F444W ≲ 28, which are also the most abundant (Fig. 9). Therefore, the performance of σ_MAD ~ 0.025 and η ~ 9% can be considered as representative of these typical, faint and numerous JWST objects.

To investigate the photo-z performance for different populations, we analyze samples selected by color and type, shown in Fig. 17. We chose the F115W-F150W and F277W-F444W colors in three selections < −0.25 (red), −0.25 < color < 0.25 (middle), and >0.25 (blue) using the SE++ model photometry. As a function of the F115W-F150W color, there is no significant change, with all three samples showing excellent performance with σ_MAD ≈ 0.011 and η < 4%. However, as a function of F277W-F444W the performance degrades going from blue to red samples with σ_MAD = 0.010 for the blue, σ_MAD = 0.015 for the middle, and σ_MAD = 0.035 for the red samples. The fraction of outliers increases from ~1.5% for the blue and middle to 7.28% for the red sample. There is also a slight bias of b = 0.014 for the m_F277W − m_F444W > 0.25 sample that does not exceed the outlier definition. This red population is dominated by highly dust-obscured galaxies that can have sub-millimeter counterparts (e.g., Barrufet et al. 2023; Gottumukkala et al. 2024), which are known to be difficult cases for SED fitting (e.g., Casey 2012; Hayward & Smith 2015).

We also assess the performance for star-forming and quiescent galaxies selected using the NUV − r − J diagram, the latter satisfying (NUV − r) > 3 × (r − J) + 1 and (NUV − r) > 3.1 (Ilbert et al. 2013). Both show excellent performance with σ_MAD = 0.008 and η = 2.26% for star-forming and σ_MAD = 0.013 and η = 1.95% for quiescent galaxies. These metrics showcase the excellent quality of our photometric redshifts with stable performance for different magnitude, color and type selected populations.

To assess the performance at the highest redshifts, we compare with the spec-z of a sample of [OIII] emitters identified in the first 10% of the data (taken in December 2024) from the COSMOS-3D NIRCam/WFSS survey (ID #5893 PI: K. Kakiichi). This includes 34 high-confidence line emitters at 6.8 ≲ z ≲ 9 (Wang et al., in prep.). Additionally, we include a compilation of spec-z available on the DAWN JWST Archive (DJA¹², Heintz et al. 2024). These include high confidence (grade=3) PRISM spectra from CAPERS (#6368, PI: Dickinson) and transient (#6585, PI: Coulter) programs. The comparison is shown in Fig. 18. There is a relatively good agreement between the photo-z and spec-z with σ_MAD = 0.038 and η = 7.9%. However, there is a median bias of b = −0.3 toward higher photo-z solutions. This spec-z sample was not used in calibrating the SED fitting with LePHARE (described in Sect. 6.1) which can be one of the reasons for the bias. Importantly, this high-z spectroscopic sample, once the COSMOS-3D survey is complete, will be highly valuable in recalibrating and improving the photo-z performance for future versions of the catalog, which we plan to build.

To demonstrate the improvement with respect to the previous COSMOS catalogs, COSMOS2020, Fig. 19 shows the σ_MAD and outlier fraction, η, as a function of magnitude for the COSMOS-Web and COSMOS2020 catalogs, computed using the same spec-z sample. Since we are matching the spec-z sample to both catalogs, we use the F444W magnitude to select in magnitude ranges for both catalogs. COSMOS-Web photo-z show better performance in both statistics, especially for fainter magnitudes, where the improvement is about a factor of two over COSMOS2020. This is unsurprising, given the much deeper JWST NIR imagaing.

In Fig. 20 we show the redshift distributions for the total, F277W magnitude-selected and spec-z samples. As expected, the median redshift increases for fainter samples. Thanks to the unique combination of survey area and depth at 1–5 μm we detect numerous candidate z > 8 galaxies. These are studied in detail in the corresponding papers (Casey et al. 2024; Franco et al. 2024; Shuntov et al. 2025; Paquereau et al. 2025, also Franco et al., in prep.).

To evaluate the quality of the photo-z uncertainties, we analyze the cumulative distribution of the ratio between |z_spec − z_phot| and the photo-z 1σ uncertainty. For well calibrated and unbiased uncertainties, the fraction of spec-z that are within the 1σ uncertainty interval $[z_{\text{phot}}^{\text{low68}},z_{\text{phot}}^{\text{upp68}}]$ should be ≈0.68. Figure 21 shows the cumulative distribution of this ratio for several magnitude selected samples, all of which reach values of ≈0.68 at a ratio of unity. This means that the flux uncertainties and consequently photo-z uncertainties are well calibrated and unbiased for all magnitudes.

Finally, we evaluate the quality of the PDF_z using the probability integral transform (PIT) statistic (Bordoloi et al. 2010). The PIT is measured as the Cumulative Distribution Function (CDF), evaluated at the z_spec of the source. The histogram of the PIT for all sources is informative on how accurate and well calibrated the PDF_z are. A flat PIT distribution indicates accurate and unbiased PDF_z; U-shaped (concave) indicates under-dispersed PDF_z; convex corresponds to overdispersed PDF_z; distributions with a slope indicate biased PDF_z. Figure 22 shows the PIT distribution for several magnitude selected samples. All magnitude-selected samples, except the one at 25 < m_F444W < 28 show positive-slope and U-shaped distributions indicating PDF_z that are slightly underdispersed and skewed toward the low-z end. This means that the PDF_z does not adequately capture the z_spec solution, and the peaks at zero and one indicate outliers that can be due to wrong SED templates or problems with the photometry. Additionally, for bright galaxies part of this can be due to the relatively coarse sampling of the redshift range with a step of 0.01, which in some cases can be larger than the photo-z uncertainty. However, since ~75% of the galaxies have 25 < m_F444W < 28, this means that most of our sources have reasonably well-calibrated PDF_z (although we note that the PIT histogram is noisier because of the smaller spec-z sample of faint galaxies)

Fig. 16 Photometric versus spectroscopic redshift comparison for mF444W < 28 (top left) and for different F444W magnitude-selected samples. The dashed teal line marks the one-to-one relation, while the dotted lines correspond to zphot > ±0.15 (1 + zspec). The number of sources in the magnitude selection (N), median absolute deviation (σ), overall bias (b), and the outlier fraction (η) are noted in each panel.

Fig. 17 Photometric versus spectroscopic redshift comparison for different color-selected samples. The top row shows F115W–F150W, while the middle row shows F277W–F444W color-selected samples at < −0.25 (left), −0.25 < color < 0.25 (middle), and > 0.25 (right column). The bottom row shows the performance for quiescent and star-forming galaxies selected using the NUVrJ diagram, i.e., (NUV − r) > 3 × (r − J) + 1 and (NUV − r) > 3.1. The dashed teal line marks the one-to-one relation, while the dotted lines correspond to zphot > ±0.15 (1 + zspec). The number of sources in the magnitude selection (N), median absolute deviation (σ), overall bias (b) and the outlier fraction (η) are noted in each panel.

Fig. 18 Photometric versus spectroscopic redshift comparison using a sample of [OIII] emitters from the COSMOS-3D NIRCam/WFSS survey, and a compilation of PRISM spectra available on DJA from CAPERS and transient programs. The dotted lines correspond to zphot > ±0.15 (1 + zspec). Points are color coded by their F444W magnitude.

Fig. 19 photo-z performance comparison with COSMOS2020. The top and bottom panels show the σMAD and outlier fraction (η) as a function of magnitude bin for COSMOS2020 (dashed lines and square symbols) and COSMOS-Web (solid lines and circle symbols). COSMOS-Web photo-z show better performance in both statistics, especially for fainter magnitudes, where the performance is better by about a factor of 2.

Fig. 20 Redshift distributions for the total (orange), spec-z (black dashed), and different magnitude-selected (blue filled) samples. The legend indicates the median redshift ⟨z⟩ for each sample.

Fig. 21 Cumulative distribution of the ratio between |zspec − zphot| and the photo-z 1σ uncertainty as a function of F444W magnitude. The 1σ uncertainty is taken as the maximum between (zphot − $z_{\text{phot}}^{\text{low68}}$) and ($z_{\text{phot}}^{\text{upp68}}$ − zphot). For unbiased photo-z uncertainties, the cumulative distributions should reach values of ~0.68 for a ratio of unity (marked by dashed lines).

6.3 Separation of stars, AGN, and galaxies

We combine several criteria to classify galaxies, stars, and AGN in the catalog. These criteria are established to maintain a good balance between purity and completeness in the galaxy sample.

First, we describe the quantities necessary to establish these criteria:

• minimum χ² is derived by LePHARE for the galaxy, stellar and AGN template libraries. These three values are noted ${\chi }_{\mathrm{g}}^2,{\chi }_{\mathrm{s}}^2,{\chi }_{\mathrm{q}}^2$, respectively.

• BzK color criteria defined by Daddi et al. (2004) is very efficient to isolate a stellar sequence in a color-color plane. We adopt colors close to Daddi et al. (2004) given our set of filters. We note Bz the magnitude difference between the HSC g and z bands, and zK the magnitude difference between the HSC z and NIRCam F277W bands. This criterion is only applied when the magnitude errors are lower than 0.3 mag in the three bands.

• The effective radius R_eff derived from the Sérsic model fitted by SE++ is used to select point-like sources.

• We define as r150 the ratio between the flux F_ν measured in circular apertures of radius $0_{\cdot }^{\mathrm{\prime }\mathrm{\prime }}25$ and $0_{\cdot }^{\mathrm{\prime }\mathrm{\prime }}1$. This criteria is also used to select point-like sources based on the light profile in the F150W filter.

Sources that satisfy any of the following criteria are classified as stars in our catalog:

• stellar best-fit template ${\chi }_{\mathrm{s}}^2<{\chi }_{\mathrm{g}}^2$ and ${\chi }_{\mathrm{s}}^2<{\chi }_{\mathrm{q}}^2$, coupled with compactness criterion R_eff < $0_{\cdot }^{\mathrm{\prime }\mathrm{\prime }}036$.

  Fig. 23 Star and galaxy separation in the g − z and z − F444W color space. The stars, selected using the criteria described in Sect. 6.3, are shown as teal dots. We only show S/N > 3 sources outside of HSC star masks.

• zK < 0 and gz > 2, which selects a region of the BzK color-color plane in which the stellar locus is well separated from the other sources. We also impose R_eff <$0_{\cdot }^{\mathrm{\prime }\mathrm{\prime }}036$.

• r150 <2 and R_eff < $0_{\cdot }^{\mathrm{\prime }\mathrm{\prime }}036$ for magnitude brighter than 27.5 in F277W. We eliminate potential AGN from this selection by imposing ${\chi }_{\mathrm{s}}^2<{\chi }_{\mathrm{q}}^2$.

We flagged sources with an AGN component that could dominate the source emission, and significantly bias the photometric redshift or the physical parameters. Sources that satisfy any of these criteria are classified as AGN in our catalog:

• ${\chi }_{\mathrm{q}}^2<{\chi }_{\mathrm{g}}^2$ and ${\chi }_{\mathrm{q}}^2<{\chi }_{\mathrm{s}}^2$, coupled with compactness criterion R_eff <$0_{\cdot }^{\mathrm{\prime }\mathrm{\prime }}036$.

• r150 <2 and R_eff < $0_{\cdot }^{\mathrm{\prime }\mathrm{\prime }}036$ for magnitude brighter than 27.5 in F277W. We impose ${\chi }_{\mathrm{q}}^2<{\chi }_{\mathrm{s}}^2$ to remove the stars.

• X-ray detection with Chandra.

Sources are classified in the catalog with the keyword type with 0 for galaxies, 1 for stars and 2 for AGN. In Fig. 23 we show the g − z − F444W color space where stars are marked with teal points, while galaxies are color coded by their z_phot. Most of the sources classified as stars fall on the expected stellar locus of the g − z − F444W color space.

We note these criteria are not exhaustive. There may be certain subtypes of QSOs/AGN that are contributing flux in the UV-MIR wavelengths, and thus the AGN Flag captures a fraction of the total AGN population. Though we note, this fraction likely represents the most luminous and relatively un-obscured QSO/AGN populations, which are the most problematic subtypes that can contaminate robust determination of galaxy properties via the multiwavelength photometry provided in this catalog.

Fig. 22 Probability integral transform (PIT) of the PDFz for magnitude-selected samples. The PIT represents the CDF of the PDFz evaluated at zspec.

7 Physical properties

7.1 LePHARE

With the LePHARE configuration adopted for COSMOS-Web, the physical parameters are derived simultaneously with the photometric redshift, as explained in Sect. 6.1. We produced a PDF_M⋆ associated with the stellar mass by summing the probabilities P_i at each stellar mass bin. By construction, these PDF_M⋆ include properly the photo-z uncertainties. However, in the case of multi-peaked PDF_z, the stellar mass point estimate derived from PDF_M⋆ would not necessarily correspond to the selected redshift point estimate.

To deal with this, we recomputed the physical parameters in a second run by setting the redshift to the point estimate established in the first run, i.e., the median of the PDF_z by default. For homogeneity, we adopt this procedure for all sources of the catalog. In this case, the uncertainties associated with these physical parameters do not propagate redshift uncertainties.

In Fig. 24 we show the stellar mass versus redshift distribution for the complete galaxy sample (cf. Sect. 6.3). This shows a smooth distribution without significant striping or clustering in specific regions in the M_⋆ − z space. We also show the stellar mass completeness computed with the Pozzetti et al. (2010) method. Briefly, we derive a rescaled stellar mass (M_resc) by scaling the F444W magnitude to the magnitude limit of the survey (m_lim) $${\log }_{10}(M_{\mathrm{r}\mathrm{e}\mathrm{s}\mathrm{c}})={\log }_{10}(M_{\star })+0.4(m_{\mathrm{F}444\mathrm{W}}-m_{\mathrm{l}\mathrm{i}\mathrm{m}}).$$(6)

Then, we define the limiting stellar mass, $M_{\star }^{\text{lim}}$, as the 90th percentile of the M_resc distribution. The limiting magnitude is obtained by comparing the COSMOS-Web F444W magnitude number counts with those in the deeper catalog in the PRIMER footprint, requiring f_compl. = N_CWeb/N_PRIMER ≈ 80%, resulting in m_lim = 27.5 (Shuntov et al. 2025). We note that this is different from the fixed-aperture completeness limit discussed in Sect. 3.4.2 because it refers to the total, model magnitude. The limiting stellar mass in several redshift bins is shown in the teal circles in Fig. 24 and the best-fit polynomial function in (1 + z) is shown in the solid line. For comparison, we also show the COSMOS2020 mass completeness function in the dashed line. Thanks to the deep 1–5 μm detection in COSMOS-Web, we achieve a mass completeness improvement by about 1 dex compared to COSMOS2020. This makes COSMOS-Web an excellent catalog to carry out statistical population studies of galaxy evolution out z ~ 10, with complete samples down to M_⋆ ~ 10⁹ M_⊙ at the highest redshifts, and down to M_⋆ ~ 10⁷ M_⊙ in the nearby Universe.

Finally, we demonstrate the morphological information available in the catalog by exploring the correlation between star formation activity and morphology. In Fig. 25 we show the NUV − r − J color distribution for 1 < z < 2 and log(M_⋆/M_⊙) > 10 galaxies, color coded by their Sérsic index (top panel) and B/T ratio (bottom panel). This shows that NUV − r − J quiescent galaxies (marked by the dashed line) are predominantly compact, with high Sérsic index (≳3.5) and B/T ratio (≳0.6). The region between the two populations, the so-called “green valley” is populated by intermediate values of n_S and B/T, while the star-forming region in predominantly disk-dominated. This is in good agreement with trends established at similar and lower redshifts by past works (e.g., Daddi et al. 2005; Trujillo et al. 2007; Barro et al. 2013, 2017; Schreiber et al. 2016). Our catalogs allow for a more detailed exploration of these morphological transformations that accompany galaxy evolution across large redshift ranges.

Fig. 24 Galaxy stellar mass versus redshift diagram for the complete sample of the COSMOS-Web catalog. The circles show the computed stellar mass completeness limit derived by rescaling the stellar masses to the limiting magnitude of the survey and taking the 90th percentile of this distribution (following Pozzetti et al. 2010), while the solid curve shows the best-fit (1 + z) polynomial. The dashed line shows the completeness in the COSMOS2020 catalog that is about 1 dex shallower in M⋆.

Fig. 25 Rest-frame NUV − r − J color distribution for 1 < z < 2 and log(M⋆/M⊙) > 10 galaxies, color coded by their Sérsic index (top) and B/T ratio (bottom). The dashed line marks the selection criteria for quiescent and star-forming galaxies (Ilbert et al. 2013). Compact galaxies (high Sérsic index and B/T ratio) predominantly lie in the quiescent region.

7.2 CIGALE

In addition to LePHARE, we derive a second set of physical properties independently through SED modeling using the CIGALE code (Boquien et al. 2019). This tool enables a detailed analysis of galaxies’ stellar populations and SFHs by fitting their observed photometry across a broad range of wavelengths. In particular, we focus on using CIGALE’s non-parametric SFH model, as this approach is highly suited to the investigation of galaxy evolution over cosmic time (Arango-Toro et al. 2023, 2025; Ciesla et al. 2025; Duan et al. 2024; Wan et al. 2024).

CIGALE utilizes a Bayesian-like methodology to model the SED of galaxies, accounting for stellar and nebular emissions, dust attenuation, and, when necessary, contributions from active galactic nuclei (AGN). The core of this approach is reconstructing the galaxy’s SFH through the sfhNlevels module (Ciesla et al. 2023, 2024), which divides the SFH into bins defined by constant SFRs. This non-parametric approach avoids biases related to specific functional forms, such as exponentially declining models, and offers a flexible framework for accurately capturing complex star formation processes. The non-parametric SFHs are structured using time bins where the SFR can vary independently, and the sfhNlevels module incorporates a continuity-burst prior (Leja et al. 2019; Tacchella et al. 2022), suppressing sudden changes in the SFR while allowing for the episodic bursts of star formation typical in high-redshift galaxies. This adaptability is crucial for accurately representing the diverse SFHs observed in the COSMOS-Web sample, which includes galaxies undergoing rapid changes due to interactions, mergers, and dynamic processes (Scoville et al. 2007; Kartaltepe et al. 2015; Davidzon et al. 2017).

CIGALE estimates a comprehensive suite of physical parameters for each galaxy by comparing its observed photometry to model predictions. These parameters, derived through a Bayesian-like analysis, include the median values and associated uncertainties for stellar masses, star formation rates (both observed and averaged over the last 100 Myr), dust attenuation estimation, and stellar population metallicity. Additionally, we have provided detailed insights into the SFHs, such as the ages since the formation of the first stellar particle or the age at which 50% of the total stellar mass was formed, the integrated SFH, and the SFR values in discrete time bins corresponding to different epochs in the galaxy’s past. We provide as well the direction and normalization of the migration vector, which characterizes the displacement of a galaxy in the stellar mass versus SFR plane over the last 250 Myr, further constraining the recent evolution of star formation activity (see Arango-Toro et al. 2025, for more details).

The accuracy of the SED fitting method is assessed using both mock catalogs and simulations. Mock catalogs, based on the actual COSMOS-Web data, are generated by adding noise to the best-fit model fluxes and re-running the fitting procedure. The results from these mock analyses demonstrate that the stellar masses and SFRs derived using CIGALE are consistent with the true values within small biases (typically <0.14 dex for stellar mass and <0.18 dex for SFR).

Further validation is carried out by comparing the SED-derived physical parameters to those obtained from the Horizon-AGN cosmological simulation (Dubois et al. 2014). The comparison shows excellent agreement, particularly for key parameters such as stellar mass and SFR, with small biases and dispersions. This further confirms the robustness of the CIGALE method in recovering physical parameters from SED fitting. These validations as well as a scientific application of the CIGALE-derived physical properties in the context of the migration of galaxies over the SFR-M_⋆ plane, are detailed in Arango-Toro et al. (2025, Sect. 4.2).

Figure 26 presents the LePHARE NUV − r − J color-color diagram for the COSMOS-Web sample, for galaxies at 1 < z < 3. Figure 27 shows the CIGALE non-parametric SFHs for starforming and quiescent galaxies given the selection criteria over the NUV-r-J diagram (Ilbert et al. 2013). We show the normalized SFHs of individual galaxies in thin colored lines and the median SFH for the whole population in solid black lines. As expected, the star-forming population shows rising and elevated SFRs, where the median SFH rises by ~0.7 dex in the first couple of Myr and then mildly decreases by ~0.3 dex since. The quiescent population shows SFHs that decline by ~1 dex in about 3.5 Gyr. There is some scatter about the median SFH for both populations, with a few individual galaxies showing rising/declining SFHs even though they are selected as quiescent/star-forming from NUV − r − J. This could be due to the limitations of NUV − r − J in efficiently distinguishing dust-reddened and star-forming from truly quiescent galaxies. Nonetheless, Fig. 27 demonstrates the consistency of the SFHs derived from CIGALE for star-forming and quiescent galaxies selected using NUV − r − J colors derived from LePHARE and the synergy of the physical information added by both SED fitting methods.

Finally, we compared the stellar mass measurements from LePHARE and CIGALE to investigate how the different assumptions in SFH and dust modeling would affect the resulting stellar masses. This comparison is shown in more detail in Shuntov et al. (2025, Appendix F). In general,they agree well, albeit with a persistent bias toward higher M_⋆ inferred by CIGALE that is dependent on both redshift and mass. The difference is around 0.1 dex at z < 0.5 and ~0.2–0.3 dex at higher redshifts. Concerning the mass, the difference is higher at the low-mass end and decreases with mass out to z ~ 5.5. At z > 5.5 this trend reverses and CIGALE results in higher masses with increasing LePHARE mass by about 0.1–0.3 dex. This trend is consistently observed when comparing stellar masses derived from non-parametric and parametric SFH (e.g., Leja et al. 2020), likely due to the ability of the former to “hide” significant stellar masses at old ages with little impact on the overall SED.

Fig. 26 Rest-frame NUV − r − J color-color diagram showing the distribution of galaxies at redshifts between 1 < z < 3 of the COSMOS-Web sample. Quiescent galaxies are highlighted in red, while star-forming galaxies are represented in blue. The dashed black line shows the selection criteria (Ilbert et al. 2013) for both populations.

Fig. 27 Star formation histories for star-forming (left) and quiescent (right) galaxies selected with the NUVrJ diagram. The colored, thin lines show the SFHs of individual galaxies, while the solid black lines show the median SFH for the whole population.

8 Conclusions

In this paper we present the COSMOS2025 catalog, the COSMOS-Web galaxy catalog of photometry, morphology, and physical parameters for over 700 000 galaxies in the central area of COSMOS. The foundation of this catalog is deep COSMOS-Web JWST imaging that provides a ~ 28 mag NIR source detection from a positive-truncated $\sqrt{{\chi }^2}$ detection image made from a combination of all four NIRCam bands. We carried out source detection using a hot and cold technique and photometric extraction using two independent approaches: aperture photometry on the space-based ACS/F814W, NIRCam, and MIRI/F770W bands and total photometry from profile-fitting using SE++ on all 37 photometric bands. These include the UltraVISTA DR6 and the HSC PDR3 data, which are both deeper and more homogenous than previous releases.

Our profile-fitting technique provides morphological measurements for every source, both using a single Sérsic together with a combination of bulge and disk light profiles. To maximize scientific applications, we also provide independent morphological measurements based on Sérsic models from Galight as well as types derived from a machine-learning classification. These independent methods show good consistency in the distribution of the morphological properties and in direct comparisons of the same quantities.

We carried out 32-band SED fitting on the total photometry to obtain photo-z and physical parameters from LePHARE. The photo-z are highly accurate with σ_MAD = 0.012 for m_F444W < 28 and stay within σ_MAD ≲ 0.03 as a function of magnitude, color, and galaxy type. We also compared our results with a specz sample of high-redshift 6.5 ≲ z ≲ 10 galaxies identified in COSMOS-3D, CAPERS, and transient programs. Photo-z performance remains excellent with σ_MAD = 0.038 and η = 7.9%. There is a bias of about −0.3 toward higher photo-z, which is likely due to the lack of these objects in our calibration sample. Thanks to the deep NIRCam detection, our catalog improves the stellar mass completeness over COSMOS2020 by about 1 dex, and is ~80% complete for log(M_⋆/M_⊙) ~ 9 at z ~ 9. We demonstrated the correlation between star-formation activity and morphology by showing that compact galaxies (high Sérsic index and B/T ratio) predominately occupy the quiescent region of the NUV − r − J rest-frame color diagram. We obtained additional physical parameters from non-parametric star formation history modeling and SED fitting from CIGALE by fixing the photo-z to the solution from LePHARE.

In our COSMOS2020 paper, we imagined what form COSMOS2025 might take: a new catalog comprising ultra-deep imaging from JWST combined with the final data releases from ground-based telescopes including UltraVISTA and HSC-SSP surveys, together with the rich heritage of existing multiwavelength COSMOS data. The catalog described here is precisely that. This paper builds on the methodological advances of COSMOS2020 (robust light profile fitting for hundreds of thousands of galaxies) and combines them with a rigorous implementation of a two-threshold detection strategy to deliver a catalog of over 700 000 galaxies with photometric measurements in 37 bands. Many papers have used these catalogs in a range of impactful studies (e.g., Casey et al. 2024; Franco et al. 2024; Akins et al. 2025; Gentile et al. 2024; Shuntov et al. 2025; Arango-Toro et al. 2025; Paquereau et al. 2025; Huertas-Company et al. 2025; Kaminsky et al. 2025; Toni et al. 2025; Yang et al. 2025; Nightingale et al. 2025). With these catalogs and images now publicly available, many unexplored scientific avenues are now open to the community.

COSMOS remains a foundation stone of extragalactic astronomy, and many more observations are either underway or planned that today give us a glimpse of what COSMOS2030 might look like. COSMOS-3D (ID #5893, PI: K. Kakiichi) is currently carrying out the largest NIRCam/WFSS survey covering 0.33 deg², including 500 arcmin² in parallel with MIRI. This will also add imaging in two additional NIRCam (F200W and F356W) and MIRI (F1000W and F2100W) bands, as well as WFSS spectra for hundreds of sources at z > 7. The HST Multi-Cycle Treasury program CLUTCH (ID #17802, PI: J. Kartaltepe) will image 0.5 deg², overlapping the areas observed by COSMOS-Web and COSMOS-3D. CLUTCH will observe the rest-frame UV to optical emission for galaxies at z < 7.5 utilizing WFC3-UVIS F225+F275W, ACS F435W+F606W+F814W, and WFC3-IR F098M. Over a wider area, the Euclid mission has already observed COSMOS in both optical and infrared bands, and these data will be delivered as part of the DR1 release in 2026. On slightly longer timescales, both the Roman Space telescope and the Rubin telescope will observe COSMOS. The methods described in this paper provide a clear path to create future COSMOS catalogs that will optimally exploit these unique data.

Data availability

We release the complete COSMOS-Web catalog, as well as the JWST/NIRCam+MIRI imaging described in Franco et al. (in prep.) and Harish et al. (in prep.) at https://cosmos2025.iap.fr/. A README is provided with detailed descriptions of the catalog columns as well as example of simple use cases. We provide an interactive map view of the COSMOS-Web field, constructed using FitsMap (Hausen & Robertson 2022), at https://cosmos2025.iap.fr/fitsmap/. This includes our COSMOS-Web catalog as an overlay, and SED plots for all 784k sources in the catalog. A copy of the catalog is available at the CDS via https://cdsarc.cds.unistra.fr/viz-bin/cat/J/A+A/704/A339.

Acknowledgements

We thank the anonymous referee for the highly constructive feedback that improved the manuscript. Support for this work was provided by NASA through grant JWST-GO-01727 awarded by the Space Telescope Science Institute, which is operated by the Association of Universities for Research in Astronomy, Inc., under NASA contract NAS 5-26555. The Cosmic Dawn Center (DAWN) is funded by the Danish National Research Foundation under grant DNRF140. This work was made possible by utilizing the CANDIDE cluster at the Institut d’Astrophysique de Paris, which was funded through grants from the PNCG, CNES, DIM-ACAV, and the Cosmic Dawn Center and maintained by S. Rouberol. French COSMOS team members are partly supported by the Centre National d’Etudes Spatiales (CNES). We acknowledge the funding of the French Agence Nationale de la Recherche for the project iMAGE (grant ANR-22-CE31-0007). This work has received funding from the Swiss State Secretariat for Education, Research and Innovation (SERI) under contract number MB22.00072. This project has received funding from the European Union’s Horizon 2020 research and innovation programme under the Marie Skłodowska-Curie grant agreement No 101148925. The data products presented herein were retrieved from the Dawn JWST Archive (DJA). D.J.A. is an initiative of the Cosmic Dawn Center (DAWN), which is funded by the Danish National Research Foundation under grant DNRF140 B.T. acknowledges support from the European Research Council (ERC) under the European Union’s Horizon 2020 research and innovation program (grant agreement number 950533). This research was supported by the Excellence Cluster ORIGINS which is funded by the Deutsche Forschungsgemeinschaft (DFG, German Research Foundation) under Germany’s Excellence Strategy - EXC 2094 - 390783311. The Hyper Suprime-Cam (HSC) collaboration includes the astronomical communities of Japan and Taiwan, and Princeton University. The HSC instrumentation and software were developed by the National Astronomical Observatory of Japan (NAOJ), the Kavli Institute for the Physics and Mathematics of the Universe (Kavli IPMU), the University of Tokyo, the High Energy Accelerator Research Organization (KEK), the Academia Sinica Institute for Astronomy and Astrophysics in Taiwan (ASIAA), and Princeton University. Funding was contributed by the FIRST program from the Japanese Cabinet Office, the Ministry of Education, Culture, Sports, Science and Technology (MEXT), the Japan Society for the Promotion of Science (JSPS), Japan Science and Technology Agency (JST), the Toray Science Foundation, NAOJ, Kavli IPMU, KEK, ASIAA, and Princeton University. D.S. carried out this research at the Jet Propulsion Laboratory, California Institute of Technology, under a contract with the National Aeronautics and Space Administration (80NM0018D0004). This paper makes use of software developed for Vera C. Rubin Observatory. We thank the Rubin Observatory for making their code available as free software at http://pipelines.lsst.io/. This paper is based on data collected at the Subaru Telescope and retrieved from the HSC data archive system, which is operated by the Subaru Telescope and Astronomy Data Center (ADC) at NAOJ. Data analysis was in part carried out with the cooperation of the Center for Computational Astrophysics (CfCA), NAOJ. We are honored and grateful for the opportunity of observing the Universe from Maunakea, which has cultural, historical and natural significance in Hawaii.

References

Appendix A Point-spread functions

Fig. A.1 Left: Median of the input star catalog vignettes. Middle: Median of the PSF models. Right: Average relative error between star and PSF model.

Appendix B Source model priors in SE++

In Fig. B.1 we show the priors adopted for the SE++ Sérsic and Bulge+Disk models (Sect. 3.5.4). 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 bulge and disk fits.

Fig. B.1 Adopted priors in SE++.

All Tables

All Figures

	Fig. 1 Survey footprint in the COSMOS field. The regions show the footprints of different imaging instruments color coded accordingly. The JWST footprint from NIRCam and MIRI is shown in the colored areas. The gray regions mark the area affected by HSC bright star masks.
In the text

	Fig. 2 Depths in all bands corresponding to 5σ. These are measured from the variance of the background measured in empty apertures averaged over the field of view. These apertures are 0.15″ in diameter in the NIRCam and HST bands, 0.5″ in MIRI, and 1″ for the ground-based data. No aperture correction is applied. The segment length corresponds to the filter width.
In the text

Fig. 3 Summary of the key steps in preparing the photometric catalog. First, source lists are generated using the “hot and cold” detection technique and used to perform aperture photometry on PSF-homogenized space-based images. Next, sources are grouped and a joint two-dimensional model light profile is extracted for all sources. Using the morphological parameters extracted from these fits, forced photometry is then performed on all bands. Finally, in the post-processing step, the photometric errors for each source are adjusted based on variance measurements in empty apertures and the source lists are cleaned and flagged. The resulting data products include morphological measurements, model-based photometry in 37 bands, and simple aperture photometry in 6 space-based bands, for more than 700 000 galaxies.

In the text

	Fig. 4 Cutout of the NIRCam ${\chi }_{+}^2$ detection image in tile B10, with the custom star mask overlaid. The mask covers bright stars and their large diffraction spikes. At the top right is a large “Dragon’s Breath” artifact (type II), which arises from the bright (K ~ 6.5 Vega mag) star to the north. This feature, and other artifacts, are manually masked.
In the text

Fig. 5 Distribution of pixels in the PSF-homogenized, noise-equalized NIRCam maps (top panels) and χ2 distribution for the quadrature sum of the four NIRCam S/N maps (bottom panel). We show the distributions without truncation (black), with positive truncation ${\chi }_{+}^2$ (blue), and with negative truncation ${\chi }_{-}^2$ (red). The χ2 distribution for degrees of freedom N = 4 and N = 2 are shown in dark and light dashed gray lines. The distribution of the sources’ maximum ${\chi }_{+}^2$ values that are recovered by the hot+cold catalog are shown in the green histogram. Vertical lines mark the equivalent of 3, 4, and 5σ thresholds for χ2 distributions with N = 2.

In the text

Fig. 6 Illustration of our “hot and cold” detection technique. Left: a 90″ × 90″ cutout of the COSMOS-Web ${\chi }_{+}^2$ detection image, showing Kron ellipses of all detected sources found from cold mode (blue) and hot mode (orange) SEP runs. Note the star in the lower-right corner that has been manually masked. Right: catalog completeness and contamination measured using an independent catalog from PRIMER (which reaches one magnitude deeper). Completeness is measured as a fraction of PRIMER sources recovered as a function of F277W magnitude (the deepest COSMOS-Web filter) for cold-mode (blue) and hot-mode (orange) sources. Contaminants are measured as sources found in COSMOS-Web but not detected in PRIMER imaging. We split these measurements into “2 exposure” (bottom row) and “4 exposure” (top row) areas where the imaging depth varies significantly. Completeness at the 8σ (5σ) survey depth exceeds 95% (80%) and contamination at these depths does not exceed 2% (4%).

In the text

	Fig. 7 Cutout of a relatively crowded region in the high-resolution NIRCam (left), the low-resolution ground-based UltraVISTA KS-band image (middle), and the partition image (right). Groups are defined using the partition image such that sources falling in the blobs with connecting pixels (here shown in the same color) are assigned the same group_id.
In the text

Fig. 8 Demonstration of the performance of our cataloging methodology with SourceXtractor++ over a relatively crowded and randomly selected 20″ × 20″ area. The top, middle, and bottom rows show the science, model, and residual images, respectively, in five wide filters: HSC g, HST/ACS F814W, UltraVISTA KS, NIRCam F277W, and MIRI F770W. The science and model images are log scaled, while the residual is linearly scaled between ±6 standard deviations of the residual cutout image. The RGB science and model images are made using F444W, F277W, and F150W.

In the text

	Fig. 9 Number counts for the SE++catalog. Left: counts in NIRCam and MIRI bands. Right: counts in HST/ACS and for ground-based broadband filters. In both cases, counts are computed in 0.35 magnitude bins and normalized by the effective area after masking, namely 0.43 deg2 in all bands except for MIRI, which covers 0.15 deg2.
In the text

Fig. 10 Photometric comparison between SE++ model magnitudes and auto magnitudes from the hot+cold catalog, as a function of auto magnitude. The auto photometry is corrected to total as described in Sect. 3.4.3. The panels show the following bands: F115W (top left), F150W (top center), F277W (top right), F444W (bottom left), HST F814W (bottom center), and F770W (bottom right). In each panel, the solid teal line indicates the running median offset, while the envelope denotes the 3σ-clipped standard deviation. The dotted white curves show the ±1σ and ±3σ envelopes of the combined photometric uncertainties. The vertical dashed teal line marks the 5σ depth in each band. Summary statistics report the median offset, Δ, and standard deviation, σ, for sources brighter than this limit. Sources with S/N < 3 and warn_flag > 2 are removed from the comparison.

In the text

Fig. 11 Photometric comparison with the COSMOS2020 THE FARMER catalog. The density histogram shows the difference between SE++ and THE FARMER as a function of THE FARMER magnitude. The solid yellow line shows the median offset, while the envelope marks the 3σ-clipped standard deviation. The dashed line shows the median offset computed down to the 5σ depth magnitude, shown by the vertical, dashed teal line. The legend shows the average offset and standard deviation for magnitudes brighter than the 5σ depth.

In the text

Fig. 12 Color comparison in the ground-based bands with the COSMOS2020 Farmer catalog. Similarly to Fig. 11, the density histogram shows the color difference between SE++ and the FARMER as a function of magnitude. The solid yellow lines and envelopes show the median offset and the 3σ-clipped standard deviation. The legend shows the averaged offset (also marked in a dashed yellow line) and the standard deviation for magnitudes brighter than the 5σ depth.

In the text

Fig. 13 Distributions and comparisons of several morphological quantities measured by SE++ and Galight for sources with with S/N > 5. Top row: density histogram of the effective radius from the Sérsic fit as a function of F277W magnitude, comparison of Reff (middle) and nS (right) between SE++ and Galight. Bottom row: histograms of the Sérsic total, bulge, and disk effective radii (left), Sérsic index (middle), and axis ratio (right) from SE++ and Galight.

In the text

	Fig. 14 Distributions of the Sérsic index (top row) and B/T in F277W (bottom row) measured from SE++ for spheroid, disk-dominated, irregular, and bulge-dominated galaxy types classified with machine learning. The dashed vertical line marks the median. These include sources with S/N > 5 classified as galaxies.
In the text

	Fig. 15 Magnitude offsets as a function of band wavelength derived as the median difference between the model and observed magnitudes with LePHARE. The shaded regions mark ±0.1, ±0.05 and ±0.01 mag. These are used as an absolute calibration of the photometry when fitting with LePHARE.
In the text

Fig. 16 Photometric versus spectroscopic redshift comparison for mF444W < 28 (top left) and for different F444W magnitude-selected samples. The dashed teal line marks the one-to-one relation, while the dotted lines correspond to zphot > ±0.15 (1 + zspec). The number of sources in the magnitude selection (N), median absolute deviation (σ), overall bias (b), and the outlier fraction (η) are noted in each panel.

In the text

Fig. 17 Photometric versus spectroscopic redshift comparison for different color-selected samples. The top row shows F115W–F150W, while the middle row shows F277W–F444W color-selected samples at < −0.25 (left), −0.25 < color < 0.25 (middle), and > 0.25 (right column). The bottom row shows the performance for quiescent and star-forming galaxies selected using the NUVrJ diagram, i.e., (NUV − r) > 3 × (r − J) + 1 and (NUV − r) > 3.1. The dashed teal line marks the one-to-one relation, while the dotted lines correspond to zphot > ±0.15 (1 + zspec). The number of sources in the magnitude selection (N), median absolute deviation (σ), overall bias (b) and the outlier fraction (η) are noted in each panel.

In the text

	Fig. 18 Photometric versus spectroscopic redshift comparison using a sample of [OIII] emitters from the COSMOS-3D NIRCam/WFSS survey, and a compilation of PRISM spectra available on DJA from CAPERS and transient programs. The dotted lines correspond to zphot > ±0.15 (1 + zspec). Points are color coded by their F444W magnitude.
In the text

Fig. 19 photo-z performance comparison with COSMOS2020. The top and bottom panels show the σMAD and outlier fraction (η) as a function of magnitude bin for COSMOS2020 (dashed lines and square symbols) and COSMOS-Web (solid lines and circle symbols). COSMOS-Web photo-z show better performance in both statistics, especially for fainter magnitudes, where the performance is better by about a factor of 2.

In the text

	Fig. 20 Redshift distributions for the total (orange), spec-z (black dashed), and different magnitude-selected (blue filled) samples. The legend indicates the median redshift ⟨z⟩ for each sample.
In the text

Fig. 21 Cumulative distribution of the ratio between |zspec − zphot| and the photo-z 1σ uncertainty as a function of F444W magnitude. The 1σ uncertainty is taken as the maximum between (zphot − $z_{\text{phot}}^{\text{low68}}$) and ($z_{\text{phot}}^{\text{upp68}}$ − zphot). For unbiased photo-z uncertainties, the cumulative distributions should reach values of ~0.68 for a ratio of unity (marked by dashed lines).

In the text

	Fig. 22 Probability integral transform (PIT) of the PDFz for magnitude-selected samples. The PIT represents the CDF of the PDFz evaluated at zspec.
In the text

	Fig. 23 Star and galaxy separation in the g − z and z − F444W color space. The stars, selected using the criteria described in Sect. 6.3, are shown as teal dots. We only show S/N > 3 sources outside of HSC star masks.
In the text

Fig. 24 Galaxy stellar mass versus redshift diagram for the complete sample of the COSMOS-Web catalog. The circles show the computed stellar mass completeness limit derived by rescaling the stellar masses to the limiting magnitude of the survey and taking the 90th percentile of this distribution (following Pozzetti et al. 2010), while the solid curve shows the best-fit (1 + z) polynomial. The dashed line shows the completeness in the COSMOS2020 catalog that is about 1 dex shallower in M⋆.

In the text

	Fig. 25 Rest-frame NUV − r − J color distribution for 1 < z < 2 and log(M⋆/M⊙) > 10 galaxies, color coded by their Sérsic index (top) and B/T ratio (bottom). The dashed line marks the selection criteria for quiescent and star-forming galaxies (Ilbert et al. 2013). Compact galaxies (high Sérsic index and B/T ratio) predominantly lie in the quiescent region.
In the text

	Fig. 26 Rest-frame NUV − r − J color-color diagram showing the distribution of galaxies at redshifts between 1 < z < 3 of the COSMOS-Web sample. Quiescent galaxies are highlighted in red, while star-forming galaxies are represented in blue. The dashed black line shows the selection criteria (Ilbert et al. 2013) for both populations.
In the text

	Fig. 27 Star formation histories for star-forming (left) and quiescent (right) galaxies selected with the NUVrJ diagram. The colored, thin lines show the SFHs of individual galaxies, while the solid black lines show the median SFH for the whole population.
In the text

	Fig. A.1 Left: Median of the input star catalog vignettes. Middle: Median of the PSF models. Right: Average relative error between star and PSF model.
In the text

	Fig. B.1 Adopted priors in SE++.
In the text