© The Authors 2026

1 Introduction

The blue sky on Earth indicates how far light can scatter from the position of an astronomical light source (the Sun). In optical imaging systems, the light from a point source is concentrated into a “point” that is strongly “spread” (or scattered) across the detector by the atmosphere, telescope, and camera as a 2D “function.” This is referred to as the point spread function (PSF).

However, the standard pipelines of wide-field surveys only provide the PSF to very small sizes (for example 10 arcsec). While this is useful for the study of barely resolved sources (for example high redshift galaxies), it does not allow for the removal of scattered light that has been produced in the image by bright sources. Estimating the extended PSF (up to several arcminutes) is crucial for an unbiased study of many aspects of galaxy evolution, including the drop in completeness of barely resolved objects near bright stars, colors of galaxies, radial profiles and truncations, measurement of the sky, halo shape, tidal features, dwarf galaxies, the circumgalactic medium, and intracluster light. (See Knapen & Trujillo 2017, for a review.)

Bright stars are a limiting factor for large imaging surveys because their light scatters to large distances and faint surface brightness levels, and their impact becomes more significant as the survey gets deeper. The influence of scattered light is so prominent that many projects are forced to mask large portions of their field that are affected. For instance, in the Hyper-Suprime-Cam Strategic Subaru Proposal (HSC-SSP), up to 20% of the total area has been masked (see Fig. 15 of Coupon et al. 2018) due to brighter stars. In the Javalambre Photometric Local Universe Survey (J-PLUS) DR3¹ (Cenarro et al. 2019), approximately 10% of the survey area is lost due to bright stars (less than in HSC-SSP because J-PLUS is shallower).

The accurate characterization of and accounting for the extended PSF is important for many aspects of imaging-based astrophysical research and has been shown to lead to better and more reliable results in areas of study ranging from point sources (Anderson & King 2000; Mighell 2005; Jiménez-Teja et al. 2015; Nardiello et al. 2022) to gravitational lensing (Gillis et al. 2020; Nardiello et al. 2022). It is also critical in studies of extended emission around astrophysical objects (e.g., Sandin 2015). For example, scattered light can create an artificial structure in galaxies that can closely mimic a real structure when it is azimuthally smooth, such as a galaxy outer disk or halo (Michard 2002; Idiart et al. 2002; Trujillo & Fliri 2016; Peters et al. 2017). The imaging of tidal structures around galaxies can also be severely affected by scattered light (van Dokkum et al. 2019; Martínez-Delgado et al. 2008; Lanzetta et al. 2024), as can the study of thick disks in galaxies (Martínez-Lombilla & Knapen 2019; Comerón et al. 2018). Furthermore, the outer parts of planetary nebulae (PNe) are affected by scattered light (Middlemass et al. 1989).

However, many surveys just characterize the central few arcseconds of the PSF, which does not allow for the appropriate removal of the scattered light produced by bright sources in or even outside the imaged area (Trujillo et al. 2001a,b). Quantifying the extended PSF (up to several arcminutes) is therefore crucial for the correct detection of the faint features that highlight important aspects of galaxy evolution and many other galactic sources.

As one pushes the detection limit to fainter surface brightness levels or observes closer to the disk of the Milky Way, the extended PSF becomes one of the most significant limiting factors. Ultimately, their increased presence is one of the factors in a survey’s design. The problem of the PSF is not limited to ground-based surveys. For example, Fig. 8 of Borlaff et al. (2022) shows the expected surface brightness of the PSF on average in the Euclid survey, and Fig. 3 of that paper illustrates how much of the Euclid survey area will be affected by stars.

The common solution of masking the bright stars is not the optimal way to deal with this issue because it leads to the loss of significant sky coverage and contiguity in a survey. Furthermore, the extended PSF is usually much larger than the masked area (which is primarily designed for point-source photometry). Historically, Moffat (1969) and Trujillo et al. (2001b) noticed that a simple Gaussian is not a good parametric fit and defined a new function for this purpose. King (1971) used the Sun to construct an extended PSF that reaches out to six degrees. de Vaucouleurs (1958) used Jupiter to measure the PSF out of five degrees. While the Sun or Jupiter can be used to construct a degree-level PSF, we cannot use such a PSF.

A more robust solution for this problem is to accurately identify the extended PSF and subtract it from bright stars (see, e.g., Mihos et al. 2013; Sandin 2014, 2015; Román et al. 2020; Infante-Sainz et al. 2020; Liu et al. 2022; Garate-Nuñez et al. 2024; Sedighi et al. 2025). There are generally two approaches to constructing the PSF. Infante-Sainz et al. (2020, hereafter I20) uses the coadd of thousands of star images (where other sources are masked) to construct a non-parametric PSF out to 8 arcmin in radius. On the other hand, Liu et al. (2022) uses a parametric Moffat and multipower-law fit to each star, which is only good when the PSF is perfectly circularly symmetric and without rings, spikes, or ghosts (artifacts that are common in large reflective telescopes).

The structure of the PSF is not fixed, and it changes with various parameters. For example, the outer part of the PSF, also known as the aureole, is affected by diffusion and reflection within the instrument (Hasan & Burrows 1995). Michard (2002) selected stars and measured the variation of the outer wings of the PSF with color, concluding that the atmospheric seeing did not change, while the PSF wings depend on the filter and on the time elapsed since the coating of the mirror. Xin et al. (2018) modeled the SDSS PSF and studied the seeing profile, and they found that the PSF changes with wavelength and time. Their PSF was constant at a timescale ranging from 5 to 30 minutes. Liu et al. (2022) showed how accumulated dust on lenses affected the extended PSF.

In this paper we aim to construct, subtract, and study the extended PSF of the J-PLUS DR3 survey. Given that J-PLUS has spikes and rings in its PSF, parametric methods are not useful. Therefore, we used the concepts introduced by I20 and greatly improved upon them for the much more complex data from J-PLUS. Furthermore, as I20 produced PSFs using the final coadded size of each T80Cam image of 2 deg² with a pixel scale of 0.55 arcsec pix⁻¹, the images are therefore 9200 × 9200 pixels (Marin-Franch et al. 2015). The survey design and the scientific aims of J-PLUS are described by Cenarro et al. (2019). In particular, the J-PLUS filters are optimized for sampling the spectral energy distribution of Milky Way stars and nearby galaxies. Also, the photometric calibration and the third data release is described by López-Sanjuan et al. (2024).

In this work, we use 935 single-exposure images from the third J-PLUS data release. This is 10% of the total amount of imaging released so far. The focus of this paper is to introduce the pipeline, analyze the J-PLUS PSF, and review some of the advantages in a subset of the tiles, not the whole survey. We hope to extend the application of this pipeline to future J-PLUS data releases, integrating it across all exposures during the data reduction process. Data reduction generally involves several steps, including dark subtraction, flat-fielding, astrometric calibration, and sky subtraction. Our aim is to incorporate PSF subtraction into the data reduction sequence, placing it before the sky subtraction stage.

The construction of the PSF requires prior knowledge of the location and some properties of the stars. For those, we used Gaia Data release 3 (Gaia Collaboration 2023).

2 Building the extended J-PLUS PSF

As summarized in the introduction, we adopt a non-parametric approach to constructing the extended PSF, based on the previous work by I20. But while I20 constructed a single PSF for the whole SDSS, this work constructs it for every exposure; many other improvements have been introduced that are further elaborated in the text. Fig. 1 summarizes all the steps for J-PLUS DR3 tile 85369 in the rSDSS filter, which contains the Coma cluster (just under the label on the top image).

The overall approach (as a summary of the subsections below) is that non-saturated stars in the image are used to construct the very central few pixels of the PSF, while the outskirts of the brightest stars are used for the outer part of the PSF, out to distances of several arcminutes. Intermediate-brightness stars are used to connect the middle and center part of the PSF, thus creating a profile which spans up to several arcminutes in radius and covering (from brightest to faintest part) up to 15 mag/arcsec² in surface brightness. These extended PSFs are created without any assumptions on the functional form of the PSF and are thus robust to any peculiarities that are present in many optical systems. The subsections below we describe the intricacies of constructing each part of the PSF.

For PSF construction, an internal sky estimate is derived and subtracted using the NoiseChisel program (Akhlaghi & Ichikawa 2015). This intermediate step should not be confused with the formal sky subtraction performed later in the data reduction pipeline, as the NoiseChisel parameters may differ between the two stages.

As shown in Fig. 1, the different components of the PSF are constructed by generating stamps around known stars. These stamp images were produced with the astscript-psf-stamp tool, as described in Sect. 10.8.3 of Akhlaghi (2024). Each stamp is a cropped image centered – sub-pixel aligned when necessary – on the given coordinates, with all foreground and background sources masked. The stamp size is consistent across filters and is manually chosen based on broadband filters, as detailed in the following subsections.

2.1 Selection of stars to use

Individual stars are necessary at very different magnitudes to construct the various parts of the extended PSF. In particular, we choose stars with a good parallax from Gaia DR3: their parallax should be greater than three times the parallax error. This is done to avoid confusion with quasars for the central and middle parts. The selected stars must be sufficiently isolated from nearby bright sources to minimize contamination. Isolation is defined using the stamp size as the reference radius: a star is considered isolated if no brighter source is detected within this radius. Such stars are then consistently used to construct the different components of the PSF. Before using the images of each star, some preprocessing was necessary: Based on the saturation and nonlinearity limits in J-PLUS, saturated pixels were masked in all exposures.

I20 explored the possibility of using Gaia broadband magnitudes for their selection and sorting of the brightest stars. They found that Gaia is incomplete for objects brighter than 7th magnitude, and no stars with magnitude below 1.7 mag are present. Furthermore, stars may be brighter or fainter in specific narrow bands, and the Gaia magnitudes of such stars may not match their actual brightness within a narrowband or medium-band image.

The current work is a proof of concept for the Javalambre Physics of the Accelerating Universe Astrophysical Survey (J-PAS; Benitez et al. 2014; Bonoli et al. 2021) survey. J-PAS uses 56 narrowband filters and will thus suffer this narrow band sorting issue much strongly. Due to the varied properties of stars across 56 narrowband filters, it is inappropriate to use Gaia broadband magnitudes for sorting stars. Given all the issues above, we chose not to use Gaia magnitudes to select and sort the stars.

Fig. 1 Flow chart of extended the PSF construction for a J-PLUS tile. All images are scaled to the color bar at the bottom. The top row shows the coadded image of three exposures which are shown in second row and the second row shows the three exposures that were used to build it. The third row shows the exposure(s) used for the various parts of the PSF in the second exposure of the second row: the single exposure itself for the center and middle parts and all exposures ±2.5 hours before and after the exposure for the outer part. The marked green and red stars are used for creating central (Sect. 2.2) and middle parts (Sect. 2.3), respectively. The stars under pink crosses (on the right) are used to create the outer PSF (Sect. 2.4). The fourth row presents a selection of the star stamps, with other sources masked, while the fifth row displays the coadded stamps. These are subsequently combined to construct the “Unified PSF,” shown in the final row. Its 1D radial profile up to 9.17 arcmin (in radius) is then projected into the circular “Projected PSF.” The surface brightness of the united PSF (from the central pixel to outer) changes about 15 mag arcsec−2 (See Sect. 2.5).

2.2 Constructing the central part of the PSF

Nonmasked (due to saturation and nonlinearity) stars in each exposure are selected to construct the central part of its PSF. Gnuastro’s Segment (Akhlaghi 2019a) and MakeCatalog (Akhlaghi 2019b) programs are used to mask other sources and find these stars. Additionally, all those with any bright stars or galaxy in the neighborhood were rejected. The stars were then sorted based on the mean value of the brightest three pixels and the brightest 100 stars were selected. As mentioned in Sect. 2.1 this allows us to not use Gaia’s magnitudes. The number of stars is high enough to reach an adequate signal-to-noise ratio (larger than 3) in the central 0.46 arcmin stamp of the PSF, the stamp size is constant in the pipeline but can be adjusted by the user. The positions of the stars thus selected are marked with a green cross in the “single exposure” box of Fig. 1.

The “stamps of center” box in Fig. 1 displays the image stamps of several non-saturated stars used to construct the central part of the composite PSF. All pixels with a signal-to-noise ratio of less than three are set to “NaN” (not a number). Only the four brightest stars out of the 100 selected stars are not masked and fully considered in the stacking procedure to construct the PSF. The number four is chosen arbitrarily only to ensuring sufficient signal in the connecting region with the next part of the PSF, which will cover the majority of these pixels. This procedure ensures that fainter stars are excluded from the outer regions, thereby preventing the addition of unwanted noise.

The next step in creating the central part of the PSF is to coadd the stamps of all individual stars. But this can only be done after they are all normalized to the same flux. Most previous approaches such as I 20 used a specific radius to obtain the normalization value in the specific radius range. Similarly, Bazkiaei et al. (2024) demonstrated that selecting an appropriate normalization radius is particularly challenging in large pipelines and plays a crucial role in creating the extended PSF for the Vera C. Rubin Observatory’s Legacy Survey of Space and Time (LSST). J-PLUS has even more complex issues because of its number and diversity of filters and the goal of this work (to have an extended PSF for each exposure).

For the normalization, we consider the brightest star as a reference. We use Gnuastro’s astscript-psf-stamp, which applies a consistent masking strategy across all steps described in this paper. This tool crops both the reference and target stars, masks extraneous objects, and centers the images with sub-pixel precision. Then, all pixels fainter than a signal-to-noise ratio of 5 are masked and we find the normalization value by dividing the sum of the remaining pixels of each star by that of its reference star. As all stars are centered in the central pixel of each stamp, this method automatically leads to the correct normalization value for each star but without having to manually specify a certain radius.

After normalization we coadd all the stamps by taking the mean after 3σ-clipping to construct the central part of the PSF (shown in the left panel of the fifth row of Fig. 1). The radial profile (using Gnuastro’s astscript-radial-profile; see Infante-Sainz et al. 2024) of the central part of the PSF is shown in Fig. 2 as the blue curve (smallest radii).

Fig. 2 Radial profiles of the coadds used for different PSF components (fifth row of Fig. 1). The blue, green, and red profiles correspond to the central, middle, and outer regions, covering radii of 0.46, 0.92 and 9.17 arcmin (in radius), respectively. Vertical lines indicate the points where these regions are joined.

2.3 Constructing the middle part of the PSF

The method for constructing the middle part of the PSF is quite similar to that used for the central part. The middle part of the PSF for each tile is also composed of the stars from that tile, but here we use the faintest stars affected by nonlinearity or saturation². Information on how these pixels are identified and masked can be found in Sect. 2.3.2 (“Saturated pixels and Segment’s clumps”) of Akhlaghi (2024).

We then match the resulting catalog with Gaia coordinates and select those stars whose parallaxes are greater than three times of the parallax error. The stars are then sorted by their signal-to-noise ratio and the 50 faintest saturated stars are selected for stacking. These stars are marked with red crosses in the “Single exposure” box of Fig. 1. The choice of 50 is arbitrary and defined by the user, and the actual number used may vary across different tiles.

The remaining procedure follows the same steps as in the construction of the central part of the PSF. Adequate faint, saturated stars are selected, the sky is subtracted, and a stamp image is generated for each (Fig. 1, middle panel of the fourth row). The stamp size is fixed at 0.92 arcmin across the pipeline, although it can be modified by the user. During the stamp creation, all extraneous objects are masked. Pixels with a signal-to-noise ratio lower than three are then removed, with the exception of the four brightest stars. As mentioned in the previous section, retaining only high signal-to-noise pixels ensures that faint sources do not contribute additional noise, particularly in the outskirts. The resulting stamps are normalized following the method outlined earlier and then co-added to form the composite middle part. The green curve in Fig. 2 displays the radial profile of the middle part of the PSF. It covers 0.92 arcmin of the total PSF.

2.4 Constructing the outer part of the PSF

The construction of the outer part of the PSF is the most novel aspect of this work, and also the most complicated. This is because the outer part of the PSF can only be observed around the brightest stars. The “All exposures in 5 hours” box in Fig. 1 outlines the construction steps of the outer part of the PSF, using the brightest stars observed. As the number of bright stars in each tile is limited and we need to increase the signal-to-noise ratio in the outer part of the PSF, we use exposures taken with the same filter within a 5-hour interval of each exposure. If there was a night with fewer valid observations, we used only the images taken during that time. Within these five hours, the telescope pointing can change, we have studied the behavior of the outer part of the PSF in different times in Sect. 4.2. In all stages of constructing the outer PSF, as in the central and middle regions, the sky is subtracted and extraneous objects are masked, following the methodology described earlier.

The outer part of the PSF, spanning a radial range of 9.17 arcminutes, is constructed using the brightest saturated stars. A crucial, but by no means trivial, step is to determine their central position in the image. We could not use Gaia for the central positions of these stars because it is not complete for the brightest stars.

To find their centers we use the bleeding area of these stars (which is nearly symmetric for the T80Cam detector; the slight asymmetry has no significant effect) as shown in Fig. 3. We crop each bleeding region, setting all saturated and bleeding pixel values to 1 and all other pixel values to NaN. The image is then collapsed along the horizontal and vertical axes producing a 1D column that shows how many masked pixels were in each. The distribution along the axes can be seen in the top and right side of the bleeding region of Fig. 3. The maximum values in X and Y axes are shown with a gray dashed line, defining the coordinates of the center of the saturated star. Note that this method is specifically designed for use with the J-PLUS detector and is only generalizable to detectors which have a symmetric bleeding pattern.

After defining the central coordinates of the saturated stars, the next critical problem with very bright stars is the selection and sorting of these bright stars. If this is not carefully accounted for and we stack stars with a wide range of brightness, an artificial truncation will occur in the outer part of the PSF caused by the noise in the outer part of images of the fainter stars. Like before, due to the diversity in the J-PLUS filters, we cannot simply use the Gaia broadband magnitudes and due to the heavy saturation and bleeding in their centers, a simple read-out of the clumps from Segment (as in the middle part) is not enough. Aperture photometry is also not reliable because of two factors: the bleeding area can be different and we need stars that are also located on an edge of an image.

To sort the bright stars by flux, the stars with a saturated area larger than 50 pixels are first selected. For clarification, saturated and nonlinear pixels are recorded as NaN, so their area is determined from the contiguous NaN regions, and their original pixel values are excluded from the analysis. In contrast, pixels identified as spikes are included in the analysis, because saturated pixels are not part of the PSF, whereas spikes are. This procedure is applied to all exposures taken within the 5-hour interval of each target. We crop the images of these stars to a size of 250 pixels, as this dimension nearly includes the stellar rings (which is discussed later but are not considered in this analysis) and remains sufficiently distant from the mirror reflection region. Finally, their radial profiles are generated. The cropped images and the radial profiles of some of these stars are shown in Fig. 4 (note that in all images, darker pixels indicate higher brightness), where we show the saturation and nonlinearity level of J-PLUS’s reduced exposure with a horizontal dashed black line. The brightest star is the one for which the smallest non-blank (saturation+nonlinearity level) radius is the most distant from the center, the purple one in the figure (note that darker pixels are brighter in this paper’s images). After identifying this radius among all the stars, and in order to avoid saturated pixels that might appear around the bleeding area, we selected a point 5 pixels beyond the smallest non-blank radius as the minimum boundary of the light gray area. This value was then multiplied by 2 to determine the maximum extent of the light gray area in Fig. 4. The flux within the interval is summed and used for sorting stars by brightness. The interval does not terminate when the purple profile concludes, as we excluded the star’s ring to avoid confusion due to the varying brightness of the ring among different stars. Among all the exposures shown in Fig. 1, we sorted and selected stars based on this methodology. All the stars used to construct the outer part of the PSF are marked with pink crosses in Fig. 1.

Another important obstacle in the construction of the outer part of the extended PSF are the internal reflections (also known as ghosts). Such reflections occur at faint levels in the vicinity of bright stars due to light bouncing around in the instrument and telescope system and are very hard to model (Slater et al. 2009). They do not originate from diffraction (and are therefore not part of the PSF). Therefore, we need to mask them for each star when constructing the outer part of the PSF.

Figure 5 illustrates the automated procedure we use for the definition and removal of internal reflections. Panel a of Fig. 5 shows one of the stars which is contaminated by the internal reflection. We first divide each image stamp into four quadrants: top right (0–90 degrees), top left (90–180 degrees), lower left (180–270 degrees) and lower right (270–360 degrees). Radial profiles are created for the flux within each quadrant (see panel b). We find the radial profile with the lowest σ-clipped mean value as the reference, as it is the cleanest from reflections. For the example star of Fig. 5, the bottom right region (270–360 degrees, colored in violet) is chosen as best and less contaminated part of the star by the reflection. Afterward, an azimuthally symmetric 2D image is created from the selected radial profile.

Subtracting this 2D from the stamp image of the star greatly enhances the ghost (see panel d of Fig. 5). We then warp the subtracted image to a scale of 1/2 (so each output pixel covers 2 × 2 input pixels) and run NoiseChisel (Akhlaghi & Ichikawa 2015; Akhlaghi 2019a,b) to detect the diffuse signal. Panel e shows that the three separate reflections of this star are now clearly detected after the use of NoiseChisel. The final result of this part of the analysis is shown in panel f, where the detected pixels are converted to the previous resolution and all the detected areas are masked in the original stamp image. While the internal reflections and wings of fainter stars have been nicely masked the spikes are kept in panel f. The technical details and implementation steps of finding the spikes and un-masking them are provided in Appendix B and Fig. B.1.

After selecting, centering, ranking the brightest stars, and masking the mirror reflections, we proceed to construct the outer part of the PSF. This process is carried out by normalizing the image stamps (see Sect. 2.2). It should be noted that during the normalization step, mirror reflections appear mostly in the outer regions of the PSF. Even after masking, pixels with a high signal-to-noise ratio remain, most of which do not include mirror reflections. These pixels are then used to determine the normalization value. Finally, normalization is followed by co-adding the 22 brightest stars (see Fig. 1). It is important to note that this number applies only to the current tile, and the number of bright stars may vary across different tiles. The resulting radial profile of the outer part of the Coma cluster PSF, derived in the rSDSS filter, is illustrated by the red line in Fig. 2.

It could happen that a bright star falls on a large-scale diffuse emission, such as cirrus, galaxy halos, or intracluster light (see Román et al. 2020 and Liu et al. 2025). If that star is the only bright star used for the outer PSF, these sources of diffuse emission will affect the extended PSF reconstruction. But we usually have more than one star in a very different part of the equatorial sky (another exposure: similar atmosphere conditions, but different RA, Dec). Since they are not affected by the same type if diffuse emission, the effect on the final coadd of the outer part will be removed (or dramatically decreased).

Fig. 3 Finding the center of saturated stars from the bleeding pixels (black in the bottom-left panel). After collapsing the pixels along the X and Y axes, the peak of the distribution shows the center coordinate of the saturated star.

Fig. 4 Sorting the brightest stars in an image by magnitude using flux at a certain radius. Colored radial profiles correspond to stars in the top of the image (same same color border). Stars are sorted based on flux within the vertical gray region (which is just beyond the saturated radius of the brightest star).

2.5 Constructing the final PSF

The united PSF for each exposure is constructed by combining the individually generated central, middle, and outer components, as described earlier. The approach used to merge these components is another key innovation of this work: instead of applying fixed radii to determine the normalization factor, we adopt a dynamic method based on overlapping pixels between the different parts of the PSF, similar to the technique outlined Sect. 2.2. Specifically, the middle component is first merged with the central component, followed by the integration of the outer component. The “United PSF” box in Fig. 1 is illustrated on the left-hand side of the last row. In the given example, the united PSF spans a surface brightness range of 15 mag arcsec⁻², extending from the center to the outer regions, and covers a radial range of 9.17 arcminutes.

Ultimately, to obtain the “Projected PSF”, the 1D radial profile of the united PSF is projected into a two-dimensional circular image using the --customtable feature of Gnuastro’s Make-Profiles (Sect. 8.1 of Akhlaghi 2024). This result is displayed on the right-hand side of the last row in Fig. 1 as “Projected PSF”. The projected PSF is necessary in this work because of the relatively small number of stars in the 5 hour interval chosen here: not fully covering the area of the PSF (white regions in the “united PSF” of Fig. 1). If there are a sufficient number of bright stars to cover the whole PSF area, the “Projected PSF” will not be necessary.

To clarify, the spikes are included in the united PSF, meaning that their pixel values are also considered in projected PSF (for more details, see Appendix B). Since the spikes are a part of the PSF and their pixel values contribute to the scaling, it is preferable to include them in both the normalization and subtraction steps. This process enables reconstruction of the PSF, particularly in regions affected by mirror reflections and saturated pixels that have been masked.

Fig. 5 Flowchart of identifying and masking internal reflections (ghosts). The star’s image (panel a) is divided into quadrants by azimuthal angle (0–90, 90–180, 180–270, 270–360 degrees). The radial profile of different angles is found (in blue, red, green, and purple) and shown in panel b. The quadrant that has a smaller signal-to-noise ratio (in this case the violet one) was chosen to create a 2D circular projection of that profile (panel c). The 2D projection is subtracted from the input, highlighting all problematic signal (including ghosts and wings of fainter stars) in panel d. We detect all of those low signal-to-noise with NoiseChisel in panel e and then mask them in panel f. The pixels that belong to spikes are re-inserted, as described in the text.

3 Subtracting the PSF

To remove the final PSF in each exposure, an extended PSF is first constructed in each exposure. The scale factor is then determined using the “United PSF”. In the next step, the “Projected PSF” is scaled for each star and subsequently subtracted. Before explaining the procedure for determining the scale factor for subtraction, it is important to note that both the background and the foreground objects surrounding each star are masked using the astscript-psf-stamp tool. This tool is particularly relevant to non-isolated stars whose outer wings overlap, as it allows other sources to be masked. In this way, accurate normalization can be performed.

To calculate this scale, we first determine the median value of the sky along with its standard deviation. The median of the sky’s standard deviation is then multiplied by 3.5 and added to the median sky value, which is defined as the threshold value. In the next step, for both the united PSF and the star image, all pixels below this threshold are masked (hereafter referred to as the threshold-PSF and threshold-star). It is important to emphasize that the normalization is not based on radius; rather, only the pixels below this threshold are retained, and the corresponding area varies according to the magnitude of the stars. To obtain the scale factor from the thresholded images we do not use the raw pixels of the united PSF and star (which can be very different), but the difference of each pixel with its 4-connected neighbors (see Fig. 6 of Akhlaghi & Ichikawa 2015) for the united PSF and star (which is more precise, since it is relative not absolute). This produces four images for each neighbor for the PSF and star (eight images in total). We then divide the PSF and star image of each neighbor to find the normalization value for that neighbor. This value is approximately the same in all pixels (which had a high signal-to-noise ratio), so, the sigma-clipped median of all the pixels is used to obtain a single value for that neighbor. Ultimately, the normalization value for a certain star is found by taking the sigma-clipped median of the four neighbor normalization values.

To subtract the stars in the field, the selected stars in a certain arbitrary magnitude range will be subtracted based on their brightness: the brightest star is first subtracted, followed by the next brightest star and so on. An example of this PSF subtraction is demonstrated in Fig. 6. NGC 4212 is shown before and after PSF subtraction, also including the scattered light field from the stars (that has been subtracted).

Panel b in Fig. 6 displays scattered light from stars located at the left edge of the image, partially outside the field of view. Despite not being fully visible, their light significantly affects the NGC 4212 galaxy. This highlights the importance of accurately constructing and subtracting the extended PSF, as well as the considerable extent to which stellar light can scatter.

Fig. 6 NGC 4212 in the J-PLUS rSDSS filter before (panel a) and after (panel c) PSF subtraction. Panel b displays the PSF that was subtracted from panel a to create panel c, which also displays the residuals. In particular, we note that the brightest star is outside the left edge of the displayed regions. Since panel b does not contain noise, its color bar does not match those of the other panels.

Fig. 7 Radial profile of three different tiles in different exposures in 12 bands. The broadband PSFs have one peak while the narrowbands one have two peaks. It is due to the ring around star. The position of the ring is changed based on the filter, becoming larger from blue to red filters.

4 Results

In this section we investigate how the PSF changes with filter, time, and charge-coupled device (CCD) position. These parameters are known to affect the shape of the PSF, but our automated method on a wide area survey such as J-PAS allows for the study of relevant characteristics in several of different scenarios.

4.1 Dependence on filter

J-PLUS employs 12 filters, including broadband and the narrowband filters. Images of stars have different properties in different filters depending on the optics of the telescope. In Fig. 7, we illustrate how the PSF varies across different filters by presenting the radial profiles of the PSFs obtained from three distinct exposures in three separate tiles.

In particular, we see a secondary peak/ring in the narrowbands, whereas in the broadband filters, either only one thicker ring is visible or the second ring is much smaller than the first ring (for example, in gSDSS and uJAVA). This peak is located within a radius of approximately 50 arcseconds in the broad uJAVA filter, whereas in the redder filters, the peaks shift to a radius of 55 to 80 arcseconds. The large peak is narrow in the narrow bands but becomes wider in the broad bands. This is because in the broad bands the wavelength range is wide enough that the two rings, whose diameter changes with wavelength, are washed out and are observed as one single ring but broader ring. The rings of the light found around bright stars are related to the diffraction pattern and the PSF shape; the diameter is proportional to the wavelength. These rings, of which there can be one, two, or more, have their origin in the obscuration caused by the secondary mirror and the intermediate baffle situated between the primary and secondary mirrors. In other telescopes we do not see these rings so significantly, because in JAST80 the secondary mirror is very large in relation to the primary mirror. The increased size of the secondary mirror is due to our optics being ‘faster’, resulting in a ‘shorter’ telescope. Also, a baffle is installed to prevent direct light from entering the camera.

Another feature is the variation in the shape of the outer PSF wing across filters and tiles. At the outermost radii of some profiles in Fig. 7, such as gSDSS and rSDSS, a drop is observed. This occurs due to over-subtraction of the sky, meaning that part of the PSF flux has been mistakenly removed as background. However, this effect is not consistent across all cases, even within the same filter; for example, the light-green profile of gSDSS does not show such a drop. This issue will be addressed in a future version of the pipeline.

In addition, the number of stars used to construct different parts of the PSF varies across tiles and filters. Figure 7 shows that, although the number of stars differs for each part, these variations do not lead to substantial differences between the resulting PSFs; the only notable difference appears in the outer region, which, as mentioned, is caused by sky over-subtraction.

Fig. 8 Effect of time on the PSF. Different stars in same position on the detector are shown imaged two hours apart in the rSDSS band. In the two stars, extra objects and internal reflection are masked, and the brightest star (blue) is normalized to the faint star (red). The size of the black circles is 2, 4, and 6 arcmin. The surface brightness radial profile of each star is shown in its own color. The surface brightness differences are shown in the lower panel, where the difference between the two stars remains nearly zero within a radius of 4 arcmin. The error bars in this panel were calculated using the standard deviation after applying median absolute deviation clipping to the surface brightness differences across every 50 radii. The 4 arcmin limit is imposed by the brightness of the fainter star (if it was brighter, we could verify at a larger distance). This shows that the PSF does not change significantly in this time interval until this radius.

4.2 Dependence of PSF on time

To generate the extended PSF, we require bright stars. However, the availability of stars that are bright enough is restricted within each exposure. We thus used stars imaged in diverse fields in a certain time range around the exposure to be corrected, to enhance the signal-to-noise ratio. We assume that the PSF is stable enough to do this, and ascertain whether the PSF undergoes changes across different time frames. In this section we analyze the variations in the PSF to acquire a comprehensive understanding of its dynamic behavior over time.

We used two stars, highlighted in red and blue in Fig. 8, situated at approximately the same CCD position but observed at different times: 2 hours apart in the rSDSS filter. We generated image stamps of both stars and mask out additional objects and mirror reflections. We depict the resulting images of the stars (indicated by red and blue) in Fig. 8, where the black circles illustrate radii of 2, 4, and 6 arcmin. After normalization, we show the radial surface brightness profiles of both stars in Fig. 8, again in red and blue, along with the variation in surface brightness at each radius.

We find that the difference in the PSF between the two stars observed two hours apart, shown in the lower panel of Fig. 8, is negligible at radii smaller than 4 arcmin. At larger radii, the discrepancy is dominated by noise from the fainter star, which may result from sky over-subtraction, as discussed in Sect. 4.1 regarding Fig. 7. Further tests will be conducted on this in future versions of this pipeline. The error bars in this panel were calculated using the standard deviation after applying median absolute deviation clipping to the surface brightness differences across every 50 radii. This indicates that stars captured at two-hours intervals do not exhibit differences in surface brightness. The PSF of J-PLUS exposures does not differ much based on time, at least in this test: the conditions for doing observations have successfully removed outliers. We note that this is a single test and the PSF may not be this stable in other exposures. In next versions, we will study the expansion of the time interval for selecting stars to derive the outer part of the PSF to be able to incorporate more stars to enhance the contiguity (no blank regions) of the united PSF and its signal-to-noise ratio.

4.3 Dependence on position

To generate the PSF, we use stars at various positions across the field of view. It is well known that reflections of the mirror vary with field of view position, but they are not part of the PSF and as shown in Fig. 5, they are completely masked in all observed stars. In this section, we investigate whether any discernible change in the PSF occurs relative to the field of view position in the outer parts of the PSF.

For this investigation, we use the tile displayed in Fig. 9, as observed in the rSDSS filter. This tile was chosen primarily for its bright stars situated at various positions across the field of view. The three stars marked with colored crosses are used for this test. The yellow star is used as the reference due to its proximity to the center of the field of view. We normalized the stars marked in blue and red after masking the mirror reflection and additional objects, to match the yellow star. The distance of the red object from yellow is 20 arcminutes, while the blue one is 30 arcminutes away from yellow. We show images of the three stars, before normalization, in the middle panel of Fig. 9, where the different black circles illustrate radii of 2, 4, and 6 arcmin.

The lower panel of Fig. 9 displays the variation in surface brightness between the blue and red stars compared to the yellow one. A caveat here is that the reference star, shown in yellow, is not the brightest one, which makes the analysis more challenging. The blue profile residual lies almost directly above the red one, which is because the blue star is shallower than the red star. In the future, this will be studied in more detail by constructing dedicated PSFs in different parts of the image from more exposures. Again, beyond a radius of about 4 arcminutes, the difference is dominated by noise and is not reliable. But within 4 arcminutes, we see that the outer PSF does not exhibit significant variation with position on the field of view. This test verifies that our usage of all stars positioned across the various field of view positions to construct the outer PSF. A crucial aspect of this work is that we entirely masked the mirror reflection in this study. Without this accurate and robust masking, the PSF would indeed undergo changes with position.

Fig. 9 Dependence of PSF on position in the field of view. This tile is chosen because it has many stars which are located in different CCD position. The red, yellow and blue stars are chosen to check how PSF will change based on the CCD position. In all of them extra object and internal reflection are masked and showed with same color in the middle. Black circles depict to 2, 4 and 6 arcminutes radius. The yellow star is considered as reference star and the other ones normalize to it and the difference of surface brightness at different radii are shown at the bottom panel. Difference at shorter radius than 5 arcmin is less than 0.5 mag/arcsec2.

4.4 Sky estimation and stacks

The extended PSF of bright stars is a considerable factor in robust sky measurements during the reduction phase of any survey (Borlaff et al. 2022; Liu et al. 2023; Watkins et al. 2024). To study this effect we subtracted the sky image provided by NoiseChisel on each exposure and coadded the three exposures. In one round there was no PSF subtraction and in the second round we subtracted the extended PSF before estimating the sky.

To show the impact of the PSF subtraction, three different exposures for one tile are coadded. Figure 10 highlights the influence of scattered light from stars on the sky map. The left panel displays the sky map from NoiseChisel before subtraction of the PSF, while the middle panel depicts the map after PSF subtraction. The right panel corresponds to the difference between the two sky maps: one before and one after PSF subtraction. The contours in the last panel represent varying levels of surface brightness (24, 23, 22 mag arcsec⁻²).

In the left-hand map, the presence of a bright star is evident, yet after PSF subtraction the sky map appears more uniform. As illustrated in Fig. 6, PSFs outside the field can also contribute, particularly when bright stars are located near the field edges. Consequently, the potential impact of stars situated beyond the field of view on sky estimation values should not be overlooked. This observation underscores the importance of PSF subtraction in accurately assessing sky values.

An important point to highlight here is that we used NoiseChisel to measure the sky. This improvement will be much more prominent for surveys that use SExtractor (Bertin & Arnouts 1996) in their sky subtraction because NoiseChisel is known to be much more precise in this (see Akhlaghi & Ichikawa 2015).

4.5 Effect on nearby galaxy surface brightness profiles

In Fig. 11, the surface brightness radial profile of the NGC 4212 (which was shown in Fig. 6) in a single exposure and in the rSDSS filter is shown before and after PSF subtraction. The red profile corresponds to the surface brightness of NGC 4212 before PSF subtraction, whereas the red profile shows it after PSF subtraction. Zooming into the central part of the profile reveals that the PSF of nearby stars also affects the central part (although indeed at a smaller scale). From the center outward, the difference becomes more significant.

It is important to note that scattered light in a galaxy profile originates not only from the stars but also from their bulges or prominent central regions. In Fig. 11, we considered only the scattered light from the stars and not from the galaxy’s center. This is because our goal here is only to remove the stars, not full scatter light. Including the scattered light from the bulge is necessary in a dedicated study of nearby galaxies, as it affects all radii at all azimuthal angles (de Jong 2008; Sandin 2014), whereas the scattered light from stars that are not exactly in the galaxy center will not have this symmetricity.

Until now, galaxies such as this (that have a bright star nearby) were usually discarded in most statistical studies during sample selection. But as shown here, through the robust estimation and subtraction of the extended PSF, we are able to include them in our sample selection.

5 Discussion and summary

The vast majority of surveys tend to focus solely on the central few arcseconds of the PSF in their standard data products. This overlooks the critical need for comprehensive characterization and removal of scattered light (PSF) to much larger distances across the field, especially light originating from bright sources within the image. Additionally, all filters exhibit a distinct PSF, adding another layer of complexity to the issue. The diverse colors exhibited by bright sources also introduce a complex, multicolor interplay within the scattered light field, potentially introducing systematic biases in both photometry and color-based measurements of background astrophysical sources.

Large-scale multicolor surveys such as J-PLUS present a unique opportunity with their wide coverage across the sky and numerous filters. The extended PSF of bright stars is a problem in the J-PLUS survey (as in other wide-field surveys, and more significant as they go deeper), affecting 10% of its total area. In this paper we have demonstrated the unique features of our pipeline to estimate and subtract the extended PSF by investigating its changes from one exposure to another, depending on the CCD position and time. We also examined the effects of PSF on the sky background and the radial profile of NGC 4212. Our pipeline is based on the earlier work of I20 but has the following novelties:

• The extended PSF is estimated and subtracted on a perexposure level;

• We do not use radii to normalize stars. Normalization happens in three different parts of the process: when coadding the stars to create each part, when uniting the various parts, and when normalizing the PSF to each star that was to be subtracted. The manually inserted radii were a major issue in J-PLUS because of the diversity of the filters (some broadband, some medium band, and some narrowband);

• A robust solution to the problem of internal reflections (ghosts) has been implemented;

• To determine the central coordinates of saturated stars, we utilized the pixel distribution along the x and y axes within the bleeding region;

• Due to the limitations caused by stellar emission and absorption lines as well as the constraints of narrowband filters, Gaia is not used for sorting bright stars. Consequently, a new algorithm has been developed to address this issue without ancillary data;

• Low-level scripts have been added to Gnuastro (Sect. 10.8 of Akhlaghi 2024) with an extensive documentation and a hands-on tutorial (Sect. 2.3 of Akhlaghi 2024) to simplify the implementation and customization of this method in other surveys.

• We constructed PSF profiles reaching up to 15 mag arcsec⁻² in surface brightness and extending up to 4 arcmin for each individual J-PLUS DR3 exposure.

The J-PLUS PSF subtraction pipeline that was introduced here is written in a generic and reproducible format using Maneage (Akhlaghi et al. 2021). While it was implemented as a proof of concept in this paper (only on five J-PLUS tiles), it is automatic, and we plan to apply it to the full J-PLUS survey in future releases. Furthermore, there are several improvements that will be implemented in future, for example, the sky subtraction issue mentioned at the end of Sect. 4.1 and optimization in speed and memory usage.

A crucial aspect that warrants detailed examination in a separate study is validation through simulations. Such simulations must account for several relevant conditions, including lunar distance, zodiacal light, atmospheric fluctuations, airglow, and stray light. This task becomes particularly complex in the case of J-PLUS, as it employs narrow-, medium-, and broadband filters, and each produces PSFs of varying shapes. Additionally, internal reflections depend on the CCD position, while diffraction spikes shift accordingly. Another topic for further study is the dependence on time and detector position in all filters (here we just focused on rSDSS).

The most significant part of the pipeline that is custom written for J-PLUS is the interface to the database for data access and downloading. The lower-level components have also been implemented in Gnuastro and are easy for any instrument to use. This will allow the adoption of the pipeline in other large surveys, such as J-PAS and the Euclid Wide Survey, or the LSST and other projects, including the James Webb Space Telescope (JWST) and Analysis of Resolved Remnants of Accreted galaxies as a Key Instrument for Halo Surveys (ARRAKIHS). Apart from how the data are accessed, how the data are named, and the format used to save the data, all the configuration parameters that need to be changed for a different survey are in the reproduce/analysis/config directory.

Through the robust PSF subtraction that we performed, almost 10% of the area of J-PLUS will become usable (sources near bright stars will not be flagged), and LSBs optimized science cases such as PNe, cirrus, and the halos of galaxies in the local Universe (among others) will become possible with J-PLUS. The density of brighter stars increases when moving to lower galactic latitudes, causing the loss of many galaxies in those regions of the sky. Furthermore, with an automated and optimized pipeline that robustly subtracts the PSF (such as the one introduced here), it will be possible to go closer to the Galactic plane. Our findings underscore the significant impact of PSF subtraction on ensuring accurate generation of sky maps and precise determination of an object’s color and thus its importance in enhancing the reliability and robustness of large surveys.

Fig. 10 Effect of PSF subtraction on the estimation of the sky level. The left panel shows the sky map from NoiseChisel before PSF subtraction, while the middle one shows it after subtraction. The right panel shows the difference between these two panels and all the contours corresponding to the surface brightness of 24, 23, and 22 mag arcsec−2.

Fig. 11 Effect of PSF on the radial profile of NGC 4212 before (red) and after (green) PSF subtraction. The 2D effect was shown in Fig. 6. The PSF affects the whole galaxy, from the central to the outer parts (much more significantly on the outer parts).

Data availability

This project was developed within the reproducible framework of Maneage (Managing data lineage, Akhlaghi et al. 2021, latest Maneage commit a575ef8, dated 12 May 2025). The implementation was carried out on an x86_64 architecture with Little Endian byte order. A complete list of the software packages employed, together with their specific versions, is provided in Appendix C. The source code for the project is publicly available on SoftwareHeritage for longevity as swh:1:dir:4382ed451649e77a56b998b9324996a283158a85³ and all the output data products are available on Zenodo:17348653⁴.

Acknowledgements

We would like to extend our sincere gratitude to Ignacio Trujillo and the anonymous referee for their invaluable comments and contributions throughout this project. The authors acknowledge PID2024-162229NB-I00, CNS2023-145339, PID2022-136505NB-I00 and ICTS-MRR-2021-03-CEFCA from the Spanish Ministry of Science and Innovation (MCIN/AEI/10.13039/501100011033). Co-funded by the European Union (EU) MSCA EDUCADO, GA 101119830. Views and opinions expressed are however those of the author(s) only and do not necessarily reflect those of the EU. Neither the EU nor the granting authority can be held responsible for them. Based on observations made with the JAST80 telescope and T80Cam camera for the J-PLUS project at the Observatorio Astrofísico de Javalambre (OAJ), in Teruel, owned, managed, and operated by the Centro de Estudios de Física del Cosmos de Aragón (CEFCA). We acknowledge the OAJ Data Processing and Archiving Unit (UPAD; Cristóbal-Hornillos et al. 2012) for reducing the OAJ data used in this work. Funding for the J-PLUS Project has been provided by the Governments of Spain and Aragón through the Fondo de Inversiones de Teruel; the Aragonese Government through the Research Groups E96, E103, E16_17R, E16_20R, and E16_23R; the MCIN/AEI y FEDER, Una manera de hacer Europa) with grants PID2021-124918NB-C41, PID2021-124918NB-C42, PID2021-124918NA-C43, and PID2021-124918NB-C44; the Spanish Ministry of Science, Innovation and Universities (MCIU/AEI/FEDER, UE) with grants PGC2018-097585-B-C21 and PGC2018-097585-B-C22; the Spanish Ministry of Economy and Competitiveness (MINECO) under AYA2015-66211-C2-1P, AYA2015-66211-C2-2, AYA2012-30789, and ICTS-2009-14; and European FEDER funding (FCDD10-4E-867, FCDD13-4E-2685). The Brazilian agencies FINEP, FAPESP, and the National Observatory of Brazil have also contributed to this Project.

References

Appendix A J-PLUS database

To get the relevant metadata of all J-PLUS exposures, we used the ADQL query below, which can be submitted in the J-PLUS webpage⁵ after login.

SELECT tile.ref_tile_id, rc.tile_id,
       rc_id, filter.name as filter_name,
       rc.RA,rc.DEC, rc.OBSERVATION_DATE,
       rc.FWHM_MEAN, rc.EXPTIME, rc.AIRMASS,
       tile.name as tile_name, tile.zpt
FROM jplus.ReducedIndividualFrame AS rc
     JOIN jplus.Filter AS filter
          ON filter.filter_id = rc.filter_id
     JOIN jplus.TileImage AS tile
          ON tile.tile_id = rc.tile_id

The output of the above query identifies each exposure with an “RC-ID.” The required exposures can be downloaded from a generic link, which must be modified by replacing the RC-ID with the corresponding exposure number. For example, to obtain the image with an RC-ID of 1040978, one would replace the RC-ID with 1040978 in the link provided.⁶

Appendix B Identifying stars spikes

In Fig. 5 we demonstrated the process of masking internal reflections (ghosts). However, the spikes produced by bright stars were also masked in the process. These spikes come from diffraction and are formed by the secondary mirror and baffle holders. Thus, the spikes are a component of the PSF. However, the reflection from the mirror does spread the light and thus is not considered part of the PSF. Therefore, spikes should not be masked out before deriving the PSF, whereas the mirror’s reflection should be masked.

When a structure depends on distance and azimuthal angle from a central point, such as the diffraction spikes of stars, projecting the pixels to polar coordinates significantly simplifies the analysis. For this project, a new option for generating polar plots has been added to Gnuastro (see Eskandarlou & Akhlaghi 2024). To restore the spikes, we create a polar plot of each star stamp after masking the internal reflection. We show in the top panel of Fig. B.1 the polar plot of the stamp image in Fig. 5(f). We then collapse the polar plot in the radius axis with a σ-clipped mean, which is shown as the red curve overlaid on the top panel of Fig. B.1. At the angles where spikes exist, a strong peak is present. Based on the discovered angles, the pixels of the spikes are removed from the mask of Fig. 5(e). The lower part of Fig. B.1 shows how the internal reflection is masked but the spikes are kept.

However, we found that the spikes from bright stars depend on the CCD position and should be incorporated during the PSF subtraction step to accurately derive the final PSF. Since the spikes are a part of the PSF, their pixel values are included throughout the process to improve normalization. Nevertheless, due to the reason mentioned above, they are currently not be subtracted in the final stage of PSF subtraction. Their full implementation in the subtraction process will be addressed in future studies.

Fig. B.1 Unmasking of the spikes in the image of bright stars using polar plots. The top panel shows the output of the polar plot. The x-axis shows the azimuthal angle, and y-axis shows the radius. The white pixels were originally bleeding pixels. The red profile shows the collapsed value at each azimuthal angle. Based on this plot, the angle of the spike was extracted and removed from the masked pixels, shown in Fig. 5(e).

Appendix C Software acknowledgement

This research was done with the following free software programs and libraries: Bzip2 1.0.8, CFITSIO 4.5.0, CMake 3.31.5, cURL 8.11.1, Dash 0.5.12, Discoteq flock 0.4.0, Expat 2.6.4, File 5.46, Fontconfig 2.16.0, FreeType 2.13.3, Git 2.48.1, GNU Astronomy Utilities 0.22.90-ccb84 (Akhlaghi & Ichikawa 2015), GNU Autoconf 2.72, GNU Automake 1.17, GNU AWK 5.3.1, GNU Bash 5.2.37, GNU Binutils 2.43.1, GNU Bison 3.8.2, GNU Compiler Collection (GCC) 14.2.0, GNU Coreutils 9.6, GNU Diffutils 3.10, GNU Emacs 28.1, GNU Findutils 4.10.0, GNU gettext 0.23.1, GNU gperf 3.1, GNU Grep 3.11, GNU Gzip 1.13, GNU Integer Set Library 0.27, GNU libiconv 1.18, GNU Libtool 2.5.4, GNU libunistring 1.3, GNU M4 1.4.19, GNU Make 4.4.1, GNU Multiple Precision Arithmetic Library 6.3.0, GNU Multiple Precision Complex library, GNU Multiple Precision Floating-Point Reliably 4.2.1, GNU Nano 8.3, GNU NCURSES 6.5, GNU Readline 8.2.13, GNU Scientific Library 2.8, GNU Sed 4.9, GNU Tar 1.35, GNU Texinfo 7.2, GNU Wget 1.25.0, GNU Which 2.23, GPL Ghostscript 10.04.0, Help2man, Less 668, Libffi 3.4.7, Libgit2 1.9.0, libICE 1.1.2, Libidn 1.42, Libjpeg 9f, Libpaper 1.1.29, Libpng 1.6.46, libpthread-stubs (Xorg) 0.5, libSM 1.2.5, Libtiff 4.7.0, libXau (Xorg) 1.0.12, libxcb (Xorg) 1.17.0, libXdmcp (Xorg) 1.1.5, libXext 1.3.6, Libxml2 2.13.5, libXt 1.3.1, Lzip 1.25, OpenSSL 3.4.0, PatchELF 0.13, Perl 5.40.1, pkg-config 0.29.2, podlators 6.0.2, Python 3.13.2, util-Linux 2.40.4, util-macros (Xorg) 1.20.2, WCSLIB 8.4, X11 library 1.8, XCB-proto (Xorg) 1.17.0, xorgproto 2024.1, xtrans (Xorg) 1.5.2, XZ Utils 5.6.3 and Zlib 1.3.1. The L^AT_EX source of the paper was compiled to make the PDF using the following packages: biber 2.20, biblatex 3.20, caption 68425 (revision), courier 61719 (revision), csquotes 5.2o, datetime 2.60, fancyvrb 4.5c, fmtcount 3.10, fontaxes 2.0.1, footmisc 7.0b, fp 2.1d, helvetic 61719 (revision), kastrup 15878 (revision), logreq 1.0, mweights 53520 (revision), newtx 1.756, pgf 3.1.10, pgfplots 1.18.1, preprint 2011, setspace 6.7b, sttools 3.1, texgyre 2.501, times 61719 (revision), titlesec 2.17, trimspaces 1.1, txfonts 15878 (revision), ulem 53365 (revision), xcolor 3.02, xkeyval 2.9, xpatch 0.3 and xstring 1.86. We are very grateful to all their creators for freely providing this necessary infrastructure. This research (and many other projects) would not be possible without them.

All Figures

Fig. 1 Flow chart of extended the PSF construction for a J-PLUS tile. All images are scaled to the color bar at the bottom. The top row shows the coadded image of three exposures which are shown in second row and the second row shows the three exposures that were used to build it. The third row shows the exposure(s) used for the various parts of the PSF in the second exposure of the second row: the single exposure itself for the center and middle parts and all exposures ±2.5 hours before and after the exposure for the outer part. The marked green and red stars are used for creating central (Sect. 2.2) and middle parts (Sect. 2.3), respectively. The stars under pink crosses (on the right) are used to create the outer PSF (Sect. 2.4). The fourth row presents a selection of the star stamps, with other sources masked, while the fifth row displays the coadded stamps. These are subsequently combined to construct the “Unified PSF,” shown in the final row. Its 1D radial profile up to 9.17 arcmin (in radius) is then projected into the circular “Projected PSF.” The surface brightness of the united PSF (from the central pixel to outer) changes about 15 mag arcsec−2 (See Sect. 2.5).

In the text

	Fig. 2 Radial profiles of the coadds used for different PSF components (fifth row of Fig. 1). The blue, green, and red profiles correspond to the central, middle, and outer regions, covering radii of 0.46, 0.92 and 9.17 arcmin (in radius), respectively. Vertical lines indicate the points where these regions are joined.
In the text

	Fig. 3 Finding the center of saturated stars from the bleeding pixels (black in the bottom-left panel). After collapsing the pixels along the X and Y axes, the peak of the distribution shows the center coordinate of the saturated star.
In the text

	Fig. 4 Sorting the brightest stars in an image by magnitude using flux at a certain radius. Colored radial profiles correspond to stars in the top of the image (same same color border). Stars are sorted based on flux within the vertical gray region (which is just beyond the saturated radius of the brightest star).
In the text

Fig. 5 Flowchart of identifying and masking internal reflections (ghosts). The star’s image (panel a) is divided into quadrants by azimuthal angle (0–90, 90–180, 180–270, 270–360 degrees). The radial profile of different angles is found (in blue, red, green, and purple) and shown in panel b. The quadrant that has a smaller signal-to-noise ratio (in this case the violet one) was chosen to create a 2D circular projection of that profile (panel c). The 2D projection is subtracted from the input, highlighting all problematic signal (including ghosts and wings of fainter stars) in panel d. We detect all of those low signal-to-noise with NoiseChisel in panel e and then mask them in panel f. The pixels that belong to spikes are re-inserted, as described in the text.

In the text

Fig. 6 NGC 4212 in the J-PLUS rSDSS filter before (panel a) and after (panel c) PSF subtraction. Panel b displays the PSF that was subtracted from panel a to create panel c, which also displays the residuals. In particular, we note that the brightest star is outside the left edge of the displayed regions. Since panel b does not contain noise, its color bar does not match those of the other panels.

In the text

	Fig. 7 Radial profile of three different tiles in different exposures in 12 bands. The broadband PSFs have one peak while the narrowbands one have two peaks. It is due to the ring around star. The position of the ring is changed based on the filter, becoming larger from blue to red filters.
In the text

Fig. 8 Effect of time on the PSF. Different stars in same position on the detector are shown imaged two hours apart in the rSDSS band. In the two stars, extra objects and internal reflection are masked, and the brightest star (blue) is normalized to the faint star (red). The size of the black circles is 2, 4, and 6 arcmin. The surface brightness radial profile of each star is shown in its own color. The surface brightness differences are shown in the lower panel, where the difference between the two stars remains nearly zero within a radius of 4 arcmin. The error bars in this panel were calculated using the standard deviation after applying median absolute deviation clipping to the surface brightness differences across every 50 radii. The 4 arcmin limit is imposed by the brightness of the fainter star (if it was brighter, we could verify at a larger distance). This shows that the PSF does not change significantly in this time interval until this radius.

In the text

Fig. 9 Dependence of PSF on position in the field of view. This tile is chosen because it has many stars which are located in different CCD position. The red, yellow and blue stars are chosen to check how PSF will change based on the CCD position. In all of them extra object and internal reflection are masked and showed with same color in the middle. Black circles depict to 2, 4 and 6 arcminutes radius. The yellow star is considered as reference star and the other ones normalize to it and the difference of surface brightness at different radii are shown at the bottom panel. Difference at shorter radius than 5 arcmin is less than 0.5 mag/arcsec2.

In the text

	Fig. 10 Effect of PSF subtraction on the estimation of the sky level. The left panel shows the sky map from NoiseChisel before PSF subtraction, while the middle one shows it after subtraction. The right panel shows the difference between these two panels and all the contours corresponding to the surface brightness of 24, 23, and 22 mag arcsec−2.
In the text

	Fig. 11 Effect of PSF on the radial profile of NGC 4212 before (red) and after (green) PSF subtraction. The 2D effect was shown in Fig. 6. The PSF affects the whole galaxy, from the central to the outer parts (much more significantly on the outer parts).
In the text

Fig. B.1 Unmasking of the spikes in the image of bright stars using polar plots. The top panel shows the output of the polar plot. The x-axis shows the azimuthal angle, and y-axis shows the radius. The white pixels were originally bleeding pixels. The red profile shows the collapsed value at each azimuthal angle. Based on this plot, the angle of the spike was extracted and removed from the masked pixels, shown in Fig. 5(e).

In the text