J-PLUS: Turning Off the Bright Stars

Astronomy &Astrophysics manuscript no. paper ©ESO 2025 October 18, 2025 J-PLUS: Turning Off the Bright Stars Sepideh Eskandarlou1, Mohammad Akhlaghi1,2, Johan H. Knapen3,4, Carlos López-Sanjuan1,2, Raúl Infante-Sainz1, Helena Domínguez Sánchez1,5, Zahra Sharbaf3,4, Héctor Vázquez Ramió1,2, Juan Antonio Fernández Ontiveros1,2, César Iñiguez García1, Tamara Civera Lorenzo1, David José Muniesa Gallardo1, Paula R.T. Coelho8, Alessandro Ederoclite1,2, Jesus Varela1, Fran Jiménez-Esteban9, A. Javier Cenarro1,2, Antonio Marín-Franch1,2, Renato A. Dupke6,7, Mariano Moles1, Carlos Hernández-Monteagudo3,4, Rahna P.T.1, David Cristóbal-Hornillos1, Jailson Alcaniz6, Laerte Sodré Jr.8, and Raul E. Angulo10,11 1Centro de Estudios de Física del Cosmos de Aragón (CEFCA), Plaza San Juan 1, 44001 Teruel, Spain 2Unidad Asociada CEFCA-IAA, CEFCA, Unidad Asociada al CSIC por el IAA y el IFCA, Plaza San Juan 1, 44001 Teruel, Spain 3Instituto de Astrofísica de Canarias, La Laguna, 38205, Tenerife, Spain 4Departamento de Astrofísica, Universidad de La Laguna, 38206, Tenerife, Spain 5Instituto de Física de Cantabria (IFCA), CSIC-Univ. de Cantabria, Avda. los Castros, s/n, E-39005 Santander, Spain 6Observatório Nacional - MCTI (ON), Rua Gal. José Cristino 77, São Cristóvão, 20921-400 Rio de Janeiro, Brazil 7University of Michigan, Department of Astronomy, 1085 South University Ave., Ann Arbor, MI 48109, USA 8Instituto de Astronomia, Geofísica e Ciências Atmosféricas, Universidade de São Paulo, 05508-090 São Paulo, Brazil 9Centro de Astrobiología, CSIC-INTA, Camino bajo del castillo s/n, E-28692, Villanueva de la Canãda, Madrid, Spain 10 Donostia International Physics Centre (DIPC), Paseo Manuel de Lardizabal 4, 20018 Donostia-San Sebastián, Spain 11 IKERBASQUE, Basque Foundation for Science, 48013, Bilbao, Spain Received July 22, 2025; accepted October 12, 2025 ABSTRACT Context. Photometric surveys require precise point spread function (PSF) characterization, as it varies across filters and is crucial for accurate photometry and low surface brightness (LSB) studies, such as galaxy halos, tidal features, extended emission from the circum-galactic medium, and intra-cluster light. However, the small PSF size provided by default pipelines suits only barely resolved objects, making it difficult to analyze regions near bright stars (rendering those regions unusable). Aims. We aim to demonstrate the feasibility of subtracting the extended PSF from each J-PLUS DR3 exposure prior to sky subtraction, an approach that has not yet been explored in wide surveys, and to enable a comprehensive analysis of its behavior across detector position and time. Methods. To build an extended, non-parametric PSF, three different ranges of stars are selected to create the central, middle, and outer regions in exposures within ∼2.5 hours of the target. These components are then combined to generate a final PSF for each exposure and filter, spanning 15 mag arcsec−2in surface brightness and 4 arcmin in radius in the broad bands. Results. In narrow-band filters, the J-PLUS PSF exhibits two rings, whereas in broad-band filters, only one ring is observed. Additionally, the position of the ring shifts with filter wavelength: as the filters become redder, the ring radius increases. We find that the precision of sky subtraction can be greatly improved with a PSF-subracted image and that out to a radius of 4 arcmin, there is no significant variation in the extended PSF observed as a function of time or position in the field of view. The radial profile of NGC 4212 (which is close to a star) is also studied before/after PSF-subtraction as a demonstration of the effect. We developed a novel method to determine the central coordinates of saturated stars, and classify stars without using Gaia magnitudes. Additionally, mirror reflections are automatically detected and masked. Furthermore, in combining different stars and various components of the PSF, we avoided the use of a fixed radius by introducing a new method that does not depend on radial measurements. Conclusions. Accurate characterization of the extended PSF and its subtraction improves sky subtraction, increases the effective area of the survey by about 10%, and enables the study of extended large LSB features in wide area surveys like J-PLUS. Our pipeline is published as free software (GNU GPLv3) an can be customized to other surveys such as J-PAS, where its impact will be even greater due to its depth. This paper is fully reproducible and produced from Commit 4860c70. Key words. Surveys – Galaxies: halos – Techniques: image processing – Methods: data analysis – Stars: imaging 1. Introduction The blue sky 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 source (for example high redshift galaxies), it does not allow the removal of the scattered light which has been produced in the image by the bright sources. Estimating the extended PSF (up to several arcminutes) is crucial for an un-biased study of many aspects of galaxy evolution. For instance, the drop in completeness of barely resolved objects near bright stars, colors of galaxies, radial profiles and truncation’s, measurement of the sky, halo shape, tidal features, dwarf galaxies, circum-galactic medium Article number, page 1 of 14 A&A proofs: manuscript no. paper (CGM), and intra-cluster light (ICL). See Knapen & Trujillo (2017) for a review. Bright stars are a limiting factor for large imaging survey because their light scatters to relatively large distances and faint surface brightness levels; becoming 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 (JPLUS) DR31(Cenarro et al. 2019), approximately 10% of the survey area is lost due to bright stars (less than in HSC-SSP because the 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). But it is critical in studies of extended emission around astrophysical objects (e.g. Sandin 2015). For example, scattered light can create artificial structure in galaxies that can closely mimic real structure, when it is azimuthally smooth, 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ínezDelgado 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). Also, the outer parts of planetary nebulae (PNe) are affected by scattered light (Middlemass et al. 1989). However, many surveys just characterize the central few arcsec region of the PSF, which does not allow 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 arcmin) is therefore crucial for the correct detection of the faint features which 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. Such that their increased presence is one of the factors in a survey’s design. The problem of the PSF is not limited to groundbased surveys. For example fig. 8 of Borlaffet 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 can not 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 1https://www.j-plus.es/datareleases/data_release_dr3 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 multi-power-law fit to each star, which is only good when the PSF is perfectly circular symmetric, without rings, spikes or ghosts (artifacts that are common in large reflective telescopes). The structure of the PSF is not fixed and 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, studied the seeing profile, and found that the PSF changes with wavelength and time. Their PSF was constant at a time scale 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 use the concepts introduced by I20 and greatly improve upon them for the much more complex data from JPLUS. Furthermore, whereas I20 produced PSFs using the final coadded size of each T80Cam image is 2 deg2with a pixel scale of 0.55 arcsec pix−1; therefore the images are 9200×9200 pixels (Marin-Franch et al. 2015). The design of the J-PLUS filters is optimized for sampling the spectral energy distribution (SED) of Milky Way (MW) stars and nearby galaxies (Cenarro et al. 2019). The camera is described in detail by Marin-Franch et al. (2015), the survey design and the scientific aims by Cenarro et al. (2019). The photometric calibration and the third data release is described by LópezSanjuan 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 a prior knowledge of the location and some properties of the stars. For those, we use Gaia Data release 3 (Gaia Collaboration et al. 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 1summarizes all the steps for J-PLUS DR3 Article number, page 2 of 14 Sepideh Eskandarlou et al.: J-PLUS: Turning Offthe Bright Stars 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 arcmin. Intermediate-brightness stars are used to connect the middle and center part of the PSF, thus creating a profile which spans up to several arcmin in radius and covering (from brightest to faintest part) up to 15 mag/arcsec2in 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 broad-band 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 pre-processing was necessary: Based on the saturation and non-linearity limit in J-PLUS, saturated pixels were masked in all the exposures used. Based on the saturation and non-linearity limits in J-PLUS, saturated pixels were masked in all exposures. I20 explored the possibility of using Gaia broad-band 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 narrow/medium-band image. The current work is a proof of concept for the Javalambre Physics of the Accelerating Universe Astrophysical Survey (JPAS; Benitez et al. 2014;Bonoli et al. 2021) survey. J-PAS uses 56 narrow-band filters and will thus suffer this narrow band sorting issue much strongly. Due to the varied properties of stars across 56 narrow-band filters, it is inappropriate to use Gaia broad-band magnitudes for sorting stars. Given all the issues above, we chose not to use Gaia magnitudes to select and sort the stars. Stack Exp-1 Exp-2 Exp-3 All exposures within 5h Single exposure Stamps of outerStamps of middleStamps of center Center Middle Outer United PSF Projected PSF 10 15 20 25 Surface Brightness (mag/arcsec2) 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/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 magarcsec−2(See Sect. 2.5). Article number, page 3 of 14 A&A proofs: manuscript no. paper 10−210−1100 10−1 100 101 102 103 104 105 106 107 Distance form center [arcmin] Relative flux Center part radial profile Middle part radial profile Outer part radial profile 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.2. Constructing the central part of the PSF Non-masked (due to saturation and non-linearity) 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. 1displays 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 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 I20 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. Fig. 3: Finding the center of saturated stars from the bleeding pixels (black in bottom-left panel). After collapsing the pixels along X and Y axis, the peak of the distribution shows the center coordinate of the saturated star. 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. 2as the blue curve (smallest radii). Article number, page 4 of 14 Sepideh Eskandarlou et al.: J-PLUS: Turning Offthe Bright Stars 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 non-linearity or saturation2. 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. 2displays 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 is a night with fewer valid observation, we use 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 Section 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 2Saturated pixels occur when a sensor reaches its maximum charge capacity, causing photons to leak into neighboring pixels. Non-linear pixels occur at fainter pixel values where the electron count does not increase at the same rate as the input flux. Therefore, these are not considered part of the PSF. 101102 101 102 103 104 Distance form center [pix] Count Blank (saturation+nonlinearity) level 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). asymmetry has no significant effect) as shown in Fig. 3. We crop each bleeding region, setting all saturated/bleeding pixel values to 1 and all other pixel values to NaN. The image is then collapsed along the horizontal/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 broad-band magnitudes and due to the heavy saturation/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 non-linear 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 Article number, page 5 of 14 A&A proofs: manuscript no. paper 0 100 200 10 15 20 25 Distance form center [arcsec] SB [mag/arcsec2] Top right Top left Bottem left Bottem right (Top right)(Top left) (Bottom right) (Bottom left) a b c fe d 20 21 22 23 24 25 26 27 28 Surface Brightness (mag/arcsec2) 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 are found (in blue, red, green, and purple) and shown in panel b. The quadrant which has less signal-to-noise ratio (in this case the violet one) is 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. 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 will be 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 non-linearity level of J-PLUS’s reduced exposure with a horizontal dashed black line. The brightest star is the one for which the smallest nonblank (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/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. Fig. 5illustrates 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 −90degrees), 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. Afterwards, 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 Article number, page 6 of 14 Sepideh Eskandarlou et al.: J-PLUS: Turning Offthe Bright Stars a b 3 arcmin c 20 25 30 Surface Brightness (mag/arcsec2) 20 25 30 15 20 25 30 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, note the brightest star which 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. 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 Band 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 signalto-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 may happen that a bright star falls on a large-scale diffuse emission, such as cirrus, galaxy halos, and intra-cluster light (ICL), 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 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). 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. 1is illustrated on the left-hand side of the last row. In the given example, the united PSF spans a surface brightness range of 15magarcsec−2, 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 MakeProfiles (Sec. 8.1 of Akhlaghi 2024). This result is displayed on the right-hand side of the last row in Fig. 1as “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. 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 Article number, page 7 of 14 A&A proofs: manuscript no. paper 0 10 20 uJAVA J0378 J0395 J0410 0 10 20 J0430 gSDSS J0515 rSDSS 10−210−1100 0 10 20 J0660 10−210−1100 iSDSS 10−210−1100 J0861 10−210−1100 zSDSS Distance from center [arcmin] Surface brightness [mag/arcsec2] Fig. 7: Radial profile of three different tiles in different exposures in 12 bands. The broad-band PSFs have one peak while the narrow-bands 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. 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. 6displays 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. 4. Results In this section, we investigate how the PSF changes with filter, time and CCD position. These parameters are known to affect the shape of the PSF but our automated method on a wide area survey like J-PAS allows us to study the relevant characteristics on a large number of different scenarios. 4.1. Dependence on filter J-PLUS employs 12 filters, including both broadand the narrow-band. 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 larger/shared 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. Article number, page 8 of 14 Sepideh Eskandarlou et al.: J-PLUS: Turning Offthe Bright Stars 16 18 20 22 24 26 02468 −2 −1 0 1 2 Distance from center (arcmin) Surface Brightness (mag/arcsec2) Fig. 8: Effect of time on the PSF: different stars in same position on detector, imaged two hours apart in rSDSS band. In both of them 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 are 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 arcminutes. The error bars in this panel were calculated using the standard deviation after applying median absolute deviation (MAD) 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. 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. Fig. 7shows 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. 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 use 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 10 arcmin 1234567 8 −2 −1 0 1 2 Distance from center (arcmin) Fig. 9: Dependence of PSF on position in 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 arcminutes is less than 0.5 mag/arcsec2 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 use 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 generate 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 4arcmin. At larger radii, the Article number, page 9 of 14