© The Authors 2025

1. Introduction

In the framework of the Λ cold dark matter (ΛCDM) model, galaxies form and evolve within dark matter (DM) halos. Their mass and size grow over cosmic time through both in situ star formation, whereby gas within the galaxy is converted into stars, and ex situ mechanisms, such as the accretion of stars and gas through mergers or interactions with neighboring systems (Toomre & Toomre 1972; White & Rees 1978; Efstathiou & Silk 1983). While mergers play a significant role in the growth of massive galaxies (Bundy et al. 2009), for most spirals and dwarf galaxies (see e.g., Martin et al. 2021) the primary driver of stellar mass assembly is the transformation of gas into stars, regulated by a critical gas density threshold (Quirk 1972; Kennicutt 1989; Schaye 2004; Roškar et al. 2008).

The existence of a gas density threshold should leave an imprint on a galaxy’s stellar mass radial density profile. Once the gas reservoir falls below a critical level, star formation ceases, leading to a characteristic “drop” in the stellar mass density. The sharpness of this edge should provide an insight into the recent star formation history of a galaxy. In fact, over time, stellar migration should smooth out these edges (Roškar et al. 2008; Sellwood 2014; Debattista et al. 2017), gradually erasing signatures of past star formation episodes. The prominence of a galaxy’s outer edge would eventually serve as a timer, with sharper edges indicating more recent star formation episodes (Martínez-Lombilla et al. 2018).

The critical gas density for star formation is theoretically estimated in the range of 3–10 M_⊙/pc² at z = 0 (Schaye 2004). Accounting for a star formation efficiency of 10–30%, one would expect the stellar mass drop to be found at around 1–3 M_⊙/pc². This is in agreement with observational findings for present-day Milky Way-like galaxies (Martínez-Lombilla et al. 2018; Díaz-García et al. 2022).

Detecting systematically such low stellar mass densities require deep imaging with surface brightness limiting depths fainter than μ_g ∼ 26–27 mag/arcsec², which is at the detection limit of the Sloan Digital Sky Survey (SDSS; Abolfathi et al. 2018). For that reason, a significant advancement in understanding the gas density threshold for star formation in disk galaxies was made using a low-surface-brightness (LSB) friendly reduction of Stripe 82 data (Trujillo & Fliri 2016; Román & Trujillo 2018), which is 2 magnitudes deeper than SDSS. Chamba et al. (2022) revealed a moderate mass dependence, with the drop in the stellar mass profile ranging from 0.6 M_⊙/pc² in dwarf galaxies to 3 M_⊙/pc² in ellipticals. Using the location of the 1 M_⊙/pc² isomass contour (R₁; Trujillo et al. 2020) as a proxy for the star formation threshold, or the location of the stellar mass density drop as a measurement of the size of the galaxy (R_edge; Chamba et al. 2022), Trujillo et al. (2020) and Chamba et al. (2022) found that the observed scatter in the stellar mass–size relation decreases by a factor larger than 2 with respect to using the effective radius, R_e. For galaxies in the mass range 10⁷ − 10¹¹ M_⊙, the stellar mass–size relation follows a power law with a slope close to 1/3. In addition, using this feature as a size measurement, Buitrago & Trujillo (2024) found an increase in the size of the disk by a factor of 2, since z ∼1 for Milky Way-like disk galaxies (M_stellar ∼ 5 × 10¹⁰ M_⊙).

The vast majority of photometric studies analyzing breaks and/or galaxy edges rely on one-dimensional (1D) surface brightness or stellar mass density profiles (Pohlen & Trujillo 2006; Bakos et al. 2008; Gutiérrez et al. 2011; Comerón et al. 2012; Martín-Navarro et al. 2012; Peters et al. 2017a,b; Trujillo et al. 2020; Chamba et al. 2022; Chamba & Hayes 2024; Buitrago & Trujillo 2024). While effective, these methods average over azimuthal variations and substructures, potentially missing key morphological features. In contrast, bidimensional (2D) analyses preserve spatial information, allowing to capture asymmetries and localized angular substructures. However, such analyses require deeper data (μ_g > 28 mag/arcsec²) to be feasible and reliable.

Optimized observational strategies and data reduction methods (Jablonka et al. 2010; Trujillo & Fliri 2016; Trujillo et al. 2021; Montes et al. 2021; Román et al. 2023; Montes et al. 2024; Golini et al. 2024) now allow for surface brightness limits of μ_g > 30 mag/arcsec² (3σ; 10″ × 10″ boxes) and similar depths are expected from the upcoming 10-year Legacy Survey of Space and Time (LSST; Ivezić et al. 2019). These deeper data are crucial for systematically investigating the edges of galaxies with an unprecedented accuracy.

Therefore, with the increasing volume of deep observational data, a robust, automated, and unbiased method of identifying galaxy edges is essential, as previous studies have primarily relied on visual detection. In computer vision, boundaries are identified by detecting rapid transitions in brightness, color, or texture (see Jing et al. 2022 for a review). However, defining a galaxy’s edge requires a physically motivated approach rather than purely visual cues. While some methods, such as fitting intersecting exponential functions (Martínez-Lombilla et al. 2023), have been proposed to objectively locate edges, a fully automated approach is crucial for ensuring consistency across large datasets. A recent study by Fernández-Iglesias et al. (2024) applied artificial intelligence to automate the identification of galaxy edges using Hubble images. By analyzing different combinations of color images, their algorithm successfully recovered previously visually identified edges (Buitrago & Trujillo 2024). Similarly, Vega-Ferrero et al. (2025) applied the segment anything model to recover the radial positions of the stellar disk edges from Buitrago & Trujillo (2024). However, machine learning models are often difficult to interpret. The goal of our work is to present a robust method of identifying galaxy edges by locating the minimum of the second derivative of the stellar mass density distribution. This technique provides a reproducible approach to measuring galaxy sizes across large datasets, facilitating a more systematic study of their outermost regions. Additionally, it will enable the investigation of the stellar mass density threshold at the edge across different morphologies, environments (Chamba & Hayes 2024), and will allow for the exploration of the cosmic evolution of galaxy sizes (Buitrago & Trujillo 2024) with the upcoming next generation of datasets.

This paper is structured as follows. Section 2 describes the methodology to locate the edge of galaxies. Section 3 presents the LBT Imaging of Galaxy Haloes and Tidal Structures (LIGHTS; Trujillo et al. 2021; Zaritsky et al. 2024) data used in this work for the galaxy NGC 3486. Section 4 details data processing and mitigation of scattered light effects. Sections 5 and 6 discuss our findings and summarize conclusions. All magnitudes are reported in the AB system (Oke & Gunn 1983).

2. A robust procedure to locate galaxy edges

We propose a method that uses the second derivative of the stellar mass density distribution, Σ_*(r), to determine a galaxy’s edge. By finding the radial position at which the second derivative reaches the outermost¹ minimum, we identify the point of curvature corresponding to the outer sharp change in slope in Σ_*(r). This point marks a natural transition in the galaxy’s structure and provides a robust identification of its edge.

2.1. Motivation for using the stellar mass density profile

If a gas density threshold is indeed present (Schaye 2004; Roškar et al. 2008) affecting star formation, both surface brightness and stellar mass density profiles should reflect a structural transition: one in terms of light (μ) and the other in terms of mass (Σ_*). In this work we use the stellar mass density profile. The reason lies in the nature of the stellar populations. In surface brightness profiles, particularly in bluer bands, variations in the luminosity of young star-forming regions, spiral arms, and clumps introduce localized fluctuations that can create multiple substructures in the profile. These can blur the detection of the outermost drop in the averaged lights profile, as they introduce local minima of the second derivative. In contrast, the stellar mass density profile, derived from the full spectral energy distribution or from color combinations, smooths out variations in stellar populations across the disk, reducing the influence of substructures that are more prominent in light-based profiles (see Fig. A.1). Working with Σ_*(r) profiles also has some drawbacks. The most important one is the addition of larger uncertainties as it relies on the combination of multiple surface brightness profiles plus the adoption of some stellar population models. For these reasons, the ideal scenario would be to work with ultra-deep near-infrared data, for which the underlying mass distributions can be traced more directly. In this work we rely on g Sloan and r Sloan data to get the stellar mass density profile of galaxies, as these are the available filters in the LIGHTS survey.

2.2. Mathematical derivation

As was mentioned in the introduction, the end of in situ star formation should lead to a drop in the stellar mass distribution of a galaxy disk. Therefore, in the majority of galaxy disks and dwarfs, the mass distribution should be well represented by a broken exponential function (a function described in Erwin et al. 2008) that has two defined exponential declining zones with different scale lengths. We define the transition point as the edge of the galaxy, R_edge.

Within the galaxy edge (r < R_edge), the slope of the disk’s stellar mass surface density profile, Σ_*(r), is parametrized by the scale length, h₁. Beyond R_edge, the profile declines with a second scale length, h₂. The sharpness of the edge is parametrized by α. For an infinitely thin disk, the broken exponential isdescribed as

$$\begin{array}{c}\hfill {\mathrm{\Sigma }}_{\ast }(r)\propto e^{-\frac{r}{h_1}}{[1+e^{\alpha (r-R_{\mathrm{edge}})}]}^{\frac{1}{\alpha }(\frac{1}{h_1}-\frac{1}{h_2})}.\end{array}$$(1)

The first derivative, $\frac{\partial }{\partial r}{\mathrm{\Sigma }}_{\ast }(r)$, with respect to the radial position indicates how the intensity changes along the radius. For a broken exponential function on a logarithmic scale, Log(Σ_*(r)), the first derivative, $\frac{\partial }{\partial r}\mathrm{Log}({\mathrm{\Sigma }}_{\ast }(r))$, shows a sharp change at the point where the slope changes, R_edge, preceded and followed by two constants (equal to – 1/h₁ before R_edge and to – 1/h₂ beyond). The second derivative, $\frac{{\partial }^2}{\partial r^2}\mathrm{Log}({\mathrm{\Sigma }}_{\ast }(r))$, gives information about the curvature or the rate of change in the slope. At the point where the slope changes abruptly, the second derivative has a minimum. Therefore, $\frac{{\partial }^2}{\partial r^2}\mathrm{Log}({\mathrm{\Sigma }}_{\ast }(r))$ presents a minimum corresponding to R_edge. In other words, the search for the minimum of $\frac{{\partial }^2}{\partial r^2}\mathrm{Log}({\mathrm{\Sigma }}_{\ast }(r))$ corresponds to the search for the location of R_edge.

2.3. Application to a simulated disk

To illustrate the implementation of our second-derivative method of edge detection, we generated a simulated disk using IMFIT (Erwin 2015), an image-fitting tool designed for modeling galaxy surface brightness profiles. The disk follows a BrokenExponential light distribution, with an inner scale length of h₁ = 45 pixels transitioning to a steeper declining outer slope of h₂ = 11 pixels at R_edge = 200 pixels. We assume a sharp edge with α = 10, a position angle (PA) of 30 deg (counterclockwise from the X axis), and an axis ratio (AR) of 0.7.

The top left panel of Fig. 1 shows the stellar mass density distribution of the simulated disk in Cartesian coordinates, where the dashed ellipse marks the point of change in the slope (R_edge). The top central panel presents the same distribution in polar coordinates and the sinusoidal wave corresponding to the edge location, which allows one to examine possible asymmetries in the edge structure.

Fig. 1. Edge detection method. The figure presents the stellar mass density map in Cartesian and polar coordinates (top left and top central panels), and the averaged surface stellar mass density profile (bottom left panel) of a simulated disk, computed in three different cases as is explained in the main text. The first derivative of the stellar mass density profile (bottom central panel) and the second derivative in both 2D (upper right) and 1D (bottom right) are also shown. Gray lines represent the surface stellar mass density and its derivatives using 10° wedges along azimuthal angles. Dashed red lines indicate Redge, corresponding to the minimum of the second derivative.

2.3.1. 1D edge detection

The bottom row of Fig. 1 summarizes the 1D analysis, showing (from left to right) the stellar mass density profile, its first derivative, and its second derivative in three different cases: the analytical BrokenExponential function used to generate the simulated disk with IMFIT, the averaged stellar mass surface density profile of the disk, and the 30-degree² polynomial that fits the averaged stellar mass surface density profile of the disk. The averaged stellar mass surface density profile of the disk was extracted using elliptical apertures with a fixed PA and inclination, and a radial step size of one pixel with the astscript-radial-profile in Gnuastro (Infante-Sainz et al. 2024).

The derivatives of the BrokenExponential function and the polynomial fit were calculated analytically. The derivatives of the averaged stellar mass density profile of the simulated disk were computed using a custom implementation of the five-point-stencil derivative³ method in Python (Sauer 2011). This numerical method leverages a higher-order Taylor series expansion, yielding a more precise derivative by drawing on multiple neighboring points. The approach reduces local errors and smooths out irregularities, though it does cause the loss of the first and last data points. However, this is inconsequential here, as our primary focus is on the outer regions of the galaxy, which are well covered within our data range. The edge location remains consistent in the three cases, demonstrating the robustness of our method using simulated data. The minimum of the second derivative corresponds to R_edge (dashed red lines in Fig. 1).

2.3.2. 2D edge detection

While fixed elliptical apertures are a standard method in surface photometry, they average over the entire azimuthal range, potentially smoothing out localized variations. In contrast, wedge⁴-based profiles allow for a more localized analysis, capturing directional differences and small-scale features, even though they still require assumptions about the geometry in each sector.

The elliptical symmetry of the simulated disk ensures that the minimum of the second derivative along each radial direction, R_edge, θ, corresponds to R_edge by design. The top right panel of Fig. 1 presents the second derivative map in polar coordinates⁵, where a clear sinusoidal pattern emerges (in pink), corresponding to the edge. This structure marks the edge in different azimuthal directions. The dispersion between gray lines in the bottom panels of Fig. 1 represents the uncertainty in computing the derivatives using 10° wedges along different azimuthal angles. No error bars are shown, as we assume an idealized scenario with infinite signal to noise (S/N). Despite the disk model being noise-free, minor inaccuracies in locating R_edge, θ and artificial oscillations of the derivatives arise from numerical artifacts introduced during the discrete generation of the 2D image, and radial or angular binning of the profile. The choice of width of the wedge, depending on the S/N of real data, is discussed in the next sections.

3. Target selection and data

We applied the methodology described above to the LIGHTS data of NGC 3486, an M33-like galaxy (i.e., with a similar stellar and total mass to M33) located at 13.6 Mpc (Kourkchi et al. 2020) and in a region of low Galactic extinction (A_g = 0.08, A_r = 0.06; Schlafly & Finkbeiner 2011). At this distance, the full width at half maximum (FWHM) of point-like sources of LIGHTS data is 1.6″ corresponding to 100 pc. NGC 3486 constitutes a case study; we aim to perform the analysis to the entire LIGHTS dataset. NGC 3486 was observed with the Large Binocular Cameras (LBCs) Blue and Red mounted at the Large Binocular Telescope (LBT) in Arizona as part of the program DDT2019B12 (PI: M. D’Onofrio). Each camera comprises four CCDs with a pixel scale of 0.224″/pix. The LBC setup includes four CCDs, with each CCD covering 7.8′ × 17.6′ section of the sky and a total camera field of view (FoV) of ∼ 23′ × 25′. For these observations, the g Sloan filter was employed in the LBC Blue camera, and the r Sloan filter was utilized in the LBC Red camera. The data collection took place on April 3, 2021, under dark sky conditions. By using dithered exposures of 180 seconds each, with four arcminutes offsets, we obtained a total on-source integration time of 1.5 hours (see Zaritsky et al. 2024 for more details on observations).

The data reduction process followed a similar strategy to the one outlined in Trujillo et al. (2021) and developed in Zaritsky et al. (2024). Our pipeline’s main steps involved deriving flat field images, background removal, astrometric and photometric calibrations, and mosaic co-addition. The final stacked data cover a region of 40.2′ × 40.2′ around NGC 3486. With this observational strategy and optimized data reduction we achieved a surface brightness limiting depth of 31.4 mag/arcsec² and 30.5 mag/arcsec² in the g and r Sloan bands, respectively, corresponding to a 3σ fluctuation of the background in areas equivalent to 10″ × 10″ (Trujillo & Fliri 2016; Román et al. 2020). Our data allow us to reach surface brightness limits about 4 mag/arcsec² deeper than the previous imaging of this particular source with SDSS. We show in Appendix B the full FoV of the LIGHTS data of NGC 3486.

4. Methods

The edge detection process for NGC 3486 follows a structured approach. We began by masking foreground and background sources to minimize contamination. Next, although LIGHTS data are background-subtracted, in some cases the subtraction was not precise enough to reliably represent the local surrounding background of the galaxy, leading to slight over- or underestimation. Therefore, we corrected for the local background. Then, we applied the point spread function (PSF) correction to account for the scattering of the galaxy’s light due to the effect of the PSF. After this, we proceeded to get the surface stellar mass density map of NGC 3486. Finally, we derived the location of the edge of NGC 3486 using the methodology described in Sect. 2. Figure 2 provides a schematic overview that serves as a reference for the methodology presented in the following sections.

Fig. 2. Flowchart illustrating the steps used to account for the effect of the PSF on the outermost regions of NGC 3486 (see main text for details). The panels display LIGHTS data of NGC 3486. Each panel covers a FoV of 15.5′ × 15.5′.

4.1. Mask and local background treatment

To achieve accurate background subtraction, we first masked foreground and background sources. The mask was generated by combining the g and r band masks. The first step in creating the masks involves running NoiseChisel and Segment (Akhlaghi & Ichikawa 2015; Akhlaghi 2019). To improve the detection of contaminating sources, we applied unsharp masking, which removes the light from the galaxy and facilitates the identification of individual contamination sources. Unsharp masking involves smoothing the original image with a Gaussian filter corresponding to the FWHM of point-like sources and subtracting this smoothed image from the original. Finally, we manually inspected the deep color image of NGC 3486 (see Fig. B.1) to mask any remaining undetected objects. The final mask applied to NGC 3486 is shown in Fig. C.1 in Appendix C.

To compute the local background correction, once the data were masked, we manually placed circular apertures (each with a 20″ radius) in specific regions approximately 6.5′ from the galaxy center. These locations were selected to be relatively free of contamination, avoiding bright sources (V < 17 mag) and diffuse light brighter than μ_g ∼ 28 mag/arcsec². We adopted this conservative approach to ensure robust background estimation. We then computed the resistant (3σ-clipped) median of the unmasked background pixels within these regions. This process was performed separately for each filter.

4.2. Removing the effect of the PSF in the outermost regions of the galaxy

Light scattering from the PSF can significantly alter the observed properties of a galaxy’s outer disk by redistributing mostly light from the central regions, leading to an artificial brightening of the outskirts (Trujillo et al. 2001; Sandin 2014, 2015). The effect depends on the instrument and wavelength used (Duc et al. 2015; Croston et al. 2006; Mihos et al. 2015; Abraham & van Dokkum 2014). In this work, we mitigated PSF effects in the galaxy’s outermost regions by combining model-based convolution techniques (Trujillo & Fliri 2016; Rich et al. 2019; Peters et al. 2017b; Martínez-Lombilla & Knapen 2019; Infante-Sainz et al. 2020) with direct image deconvolution methods, such as PyOperators (Karabal et al. 2017). The following sections detail the steps taken in this process.

4.2.1. Breaking the PSF-Stellar halo degeneracy

The problem: It can be challenging to distinguish a faint stellar halo present around galaxies from the light-broadening effects caused by the PSF. This confusion, which we term PSF-Stellar halo degeneracy, makes it difficult to tell if a galaxy’s outer light distribution is due to an extended stellar halo or merely an artifact of the PSF. To illustrate the problem, we generated two BrokenExponential models using the makeimage task in IMFIT, with the same PA, AR, and outer edge radius, R_edge. Model1 has an outer slope (h₂) that is 20% shallower than the one of Model2. In addition, Model2 is composed of a faint Exponential stellar halo component. We applied the extended g band PSF of the LBC (Sedighi et al. 2025) to both models. The averaged stellar mass surface densities of Model1 and Model2 are illustrated in the left and central panels of Fig. 3, respectively, along with their components and the profile of the PSF. The parameters used to create the models are reported in Appendix D.

Fig. 3. Surface stellar mass density profiles of two simulated disks. Left: model1, consisting of a broken exponential disk only. Center: model2, similar to Model1, but with a steeper outer disk slope (h2) and a faint stellar halo. Right: comparison of the PSF-convolved profiles for both models.

The outer profile (R > 200 pix) in both models can be interpreted as either only a PSF effect on the declining disk light or a faint stellar halo combined with the PSF. This uncertainty highlights the PSF-Stellar halo degeneracy, making it difficult to distinguish between PSF broadening and a different outer profile, h₂, plus a faint stellar halo. This test can be extended to higher inclinations, which amplify PSF effects, emphasizing the need for accurate PSF deconvolution when analyzing faint outer structures in galaxies.

In summary, the PSF redistributes photons, artificially broadening the outer regions of galaxies and smoothing the slope of the outer disk (h₂). In general this depends on the shape of the PSF. This redistribution does not alter the location of the edge in the stellar mass density distribution. However, it is key to correct for the PSF to ensure accurate recovery of the intrinsic shape of the profile, particularly the outer slope (h₂).

Wiener-Hunt deconvolution to break the degeneracy: We used the extended PSF of the LIGHTS survey (Sedighi et al. 2025) to deconvolve the background-subtracted data by applying the Wiener-Hunt (Orieux et al. 2010) deconvolution algorithm. To mathematically frame the effect of the PSF, the observed data g(x, y) from the telescope can be expressed as a convolution of the “true” image, f(x, y), with the PSF, h(x, y), plus an additive noise term, η(x, y):

$$\begin{array}{c}\hfill g(x,y)=f(x,y)\ast h(x,y)+\eta (x,y).\end{array}$$(2)

In an ideal, noise-free scenario, we can use Fourier transforms to work with these data components in frequency space. Specifically, the Fourier transforms, G(u, v) of g(x, y) and H(u, v) of h(x, y), would allow us to isolate F(u, v), the Fourier transform of the true image. From this, we could reconstruct f(x, y) by transforming back into spatial space. However, because H(u, v) acts as a low-pass filter, it tends to amplify high-frequency noise. Thus, a noise suppressor is needed. The deconvolution process using Wiener-Hunt filter was performed in the Fourier domain:

$$\begin{array}{c}\hfill F(u,v)=\frac{G(u,v)H^{\ast }(u,v)}{{|H(u,v)|}^2+K},\end{array}$$(3)

where H^*(u, v) is the complex conjugate of H(u, v) and K is the regularization parameter that controls noise suppression. In this work we used Wiener-Hunt deconvolution implemented in the Python package skimage⁶, in which the regularization parameter is automatically estimated using a stochastic iterative process (Orieux et al. 2010).

To visualize the effect of performing the Wiener-Hunt deconvolution on the data, we extracted the surface brightness profiles for NGC 3486 in both g and r bands with astscript-radial-profile and a bin size of 10″⁷. We extracted the light profiles of the galaxy using AR = 0.765 ± 0.020 and a PA of 13.0 ± 0.7° (clockwise from the X axis). This choice of parameters is justified in the Appendix E and compared with previous estimates. To maintain consistency, we applied the same mask used for sky background estimation (see Fig. C.1). The errors on these profiles were calculated by combining the uncertainties of each elliptical annulus with the one of the local sky estimate, as is outlined in Golini et al. (2024).

In Fig. 4, we show in green the resulting g and r surface brightness profiles prior to deconvolution, which exhibit an inner component extending to about 40″, followed by a exponential decline. At around 200″, a slope change occurs, with the profile steepening in both bands, marking the edge of the visible disk seen in Fig. B.1.

Fig. 4. Surface brightness profiles of NGC 3486 in the g band (upper panel) and r band (central panel) using LIGHTS data. The surface brightness profiles shown correspond to the background-subtracted data (green), the Wiener-Hunt deconvolved data (blue), and the data for which the PSF effect has been removed (orange). These profiles were obtained using a fixed ellipticity and PA representative of the outer parts of the galaxy (see the main text for details). Only values above 3σ the noise of the data are shown. The pink lines correspond to the IMFIT models we used to model the Wiener-Hunt deconvolved light distribution of the galaxy. The vertical dashed black lines highlight the location of Redge with the corresponding shaded gray uncertainty, dashed vertical purple lines highlight the location of RPSFeff (the radius where the PSF starts to be noticeable), and horizontal blue lines represent the data surface brightness limiting depth. The g − r average radial profiles prior and after deconvolution are shown in the bottom panel.

The Wiener-Hunt deconvolved surface brightness profiles in the g and r bands are also shown in Fig. 4. Compared to the original (background-subtracted) profiles, the deconvolved profiles exhibit a steeper outer slope beyond R = 190″. We visually define this transition point as R_PSFeff. This effect is due to the redistribution of light from the galaxy’s outskirts to the inner regions, compensating for the PSF-induced broadening in the original data. For instance, post-deconvolution, the inner region (R < 20″) appears 0.1 mag brighter, while outer regions (R > 190″) are up to 1 mag fainter. The edge location is preserved.

The Wiener-Hunt filter is a widely used deconvolution technique, but like other deconvolution methods such as Richardson-Lucy (Richardson 1972; Lucy 1974) or ones based on spatial regularizations such as the assumption of sparsity in a wavelet domain (Starck & Murtagh 1994; Starck et al. 2006), it introduces image smoothness, correlates pixel values, and alters background statistics, especially in regions of LSB (Murli et al. 1999; Orieux et al. 2010; Ramos et al. 2020; Yoo & Ahn 2014).

Despite this drawback, the Wiener-Hunt filter remains robust to saturation effects, PSF model inaccuracies, and varying background noise levels when averaging surface brightness distributions. These tests are illustrated in Appendix F along with a comparison to wavelet-regularized deconvolution. However, due to the effects of the deconvolution methods on the LSB structures, the treatment of the outermost region of the galaxy requires a more sophisticated procedure, which we explain below.

4.2.2. Creation of a model of the galaxy not affected by PSF

To remove the PSF effect from the NGC 3486 data, while preserving the noise structure of the LIGHTS image, we proceeded as follows. First, we modeled the averaged surface brightness profile of NGC 3486 in both bands after Wiener-Hunt deconvolution (blue profiles in Fig. 4) using a custom implementation of the IMFIT Exponential function in Python.

We modeled the Wiener-Hunt deconvolved profile in the range between R_PSFeff (vertical dashed purple lines in Fig. 4), which is the point at which the PSF begins to significantly affect the outermost regions of the light profiles in the g and r filters (i.e., > 25 mag/arcsec²), to the radius where the uncertainty reaches 0.1 mag. This results in an outer exponential scale length (h₂) of 18.8″ in both filters. Second, we generated a 2D image of this Exponential model using makeimage. When constructing the 2D outer disk component, we adopted the PA and inclination parameters obtained in Appendix E. We did not include a bulge or stellar halo component to avoid imposing a priori assumptions. Stellar halos often exhibit asymmetries, spatial offsets, or orientations misaligned with the disk, making them difficult to model reliably. Our final model of the galaxy, not affected by the PSF, consists of two parts: an inner region (R < R_PSFeff) represented by the Wiener-Hunt deconvolved image, which enhances spiral structures that would otherwise be difficult to model, and an outer region (R > R_PSFeff) that transitions to an exponential declining profile (shown in pink in Fig. 4).

4.2.3. Subtraction of the model and residuals addition: PSF deconvolved data

After generating the two deconvolved models (one for g and one for the r band) of NGC 3486, we employed the astconvolve tool in Gnuastro to convolve these models with the LIGHTS’ PSF. To isolate structural components beyond the PSF effect, we subtracted these convolved images from the original galaxy data using the astarithmetic tool in Gnuastro, resulting in a “residuals” image. These residuals highlight the galaxy’s asymmetric faint features in the outer region of the galaxy that are not accounted for in the symmetric model.

Following this, we added the residuals image back into the model created in Sect. 4.2.2, resulting in a refined image with PSF effects effectively removed. This stage of analysis ensures that any faint residual light detected in the galaxy’s outer regions corresponds to real, low-luminosity stellar structures rather than artifacts introduced by the PSF.

The surface brightness profiles, after PSF correction, extracted using fixed ellipses, are shown in orange in Fig. 4. In both bands, the galaxy light distribution after PSF removal shows a steeper decline, up to ∼ one magnitude fainter, starting at R ∼ 190″. We remark that, to avoid contamination, all surface brightness profiles shown in Fig. 4 were computed after applying a mask (see Fig. C.1) to the images.

4.3. g − r color and surface stellar mass density

With PSF effects removed, the galaxy’s outer light distribution now reliably represents the intrinsic stellar light of the galaxy and we can proceed to the analysis of the g − r color and stellar mass density profiles. The averaged g − r color profile (shown in the lower panel of Fig. 4) of NGC 3486 shows a typical U-shaped profile (Bakos et al. 2008; Azzollini et al. 2008) with an inner red bulge (g − r ∼ 0.7) transitioning to a bluer disk with an average g − r ∼ 0.3 at R = 100″ and then reddening to g − r ∼ 0.45 beyond R = 230″.

We derived the 2D surface mass density map of NGC 3486 using the mass-to-light ratio (M/L)_λ versus color relation prescribed in Roediger & Courteau (2015). Explicitly we used the g Sloan data as reference as it is our deepest imaging dataset and the g − r color. We assumed a Chabrier initial mass function (IMF; Chabrier 2003). The color image was created by subtracting the g and r surface brightness 2D maps and correcting for Galactic extinction using the following coefficients: A_g = 0.08 and A_r = 0.06 (Schlafly & Finkbeiner 2011). Inclination correction was included following the prescriptions in Trujillo et al. (2020). The resulting Cartesian and polar maps of the stellar surface mass distribution of NGC 3486 are shown in the top left and top central panels of Fig. 5, respectively.

Fig. 5. Edge detection in NGC 3486. Upper row: stellar mass surface density map in Cartesian and polar coordinates (left and central panels) and polar map of its second derivative (right). Bottom row: radial profile of averaged stellar mass surface density as a function of radius (central panel) along with its second derivative (right panel). A polynomial fit to the data and its second derivative are also shown (blue). In all panels shaded regions indicate measurement uncertainties. The outer boundary, Redge is indicated in all panels by a dashed black line. Σ*(Redge) is shown as a horizontal red line in the bottom central panel. Yellow dots in the upper panels mark Redge, θ.

After constructing the stellar surface mass density map, we applied the mask generated in Sect. 4.1 to exclude sources not belonging to the galaxy. Then, to extract the radial profile, we used fixed elliptical apertures as in Sect. 4.2.1. The resulting stellar mass density profile is presented in orange in the bottom central panel of Fig. 5. The profile reveals an inner luminous component, followed by a gradual exponential decline out to R ∼ 200″. Beyond this radius, the mass density profile steepens.

4.4. One-dimensional edge detection: R_edge

We computed the numerical derivatives of the averaged stellar mass density profile of NGC 3486. The minimum of the second numerical derivative lies at R_edge, num = 213″ ± 5″⁸. To test the result, we fit a polynomial function to the average stellar mass density profile and we computed its analytical derivatives (its second derivative is shown as a dashed blue line in the bottom right panel of Fig. 5). The minimum obtained from the polynomial fit occurs at R_edge, fit = 198″. Considering the uncertainty between these two approaches, we adopted a final R_edge = 205″ ± 5″, indicated by the dashed black line. This radius corresponds to a mean surface stellar mass density of 1.3 ± 0.1 M_⊙/pc², a μ_r(R_edge) = 25.60 ± 0.01 mag/arcsec² and a μ_g(R_edge) = 25.95 ± 0.01 mag/arcsec². We discuss the implications of these values in the following sections.

We illustrate R_edge in two ways for a more intuitive understanding. First, in the upper left panel of Fig. 5, R_edge isrepresented as an ellipse. Second, R_edge is represented as a sinusoidal wave on the polar map of the galaxy (central upper panel of Fig. 5).

In the top right panel of Fig. 5, we show the 2D map of the second derivative, representing changes in the angular direction and highlighting the natural rise of the edge of the galaxy. The location of the edge is also shown as a vertical dashed black line in Fig. 4. Visual inspection tends to place the edge at a smaller radius, especially in profiles with smoother transitions between the inner and outer exponential regions (i.e., a lower α parameter in Eq. (1)). This is a perceptual bias: smoother transitions can mislead the eye into identifying a drop earlier than it actually occurs. In contrast, our edge detection method uses the second derivative to pinpoint the inflection point (see Appendix G for a visual example of this effect).

4.5. Bidimensional edge detection: R_edge, θ

In order to characterize the edge of the galaxy in two dimensions, we considered the following. To estimate the uncertainty in the intensity profile at each radial location (E_tot, r), we computed the standard deviation of the counts within anelliptical aperture σ_tot(r). Under Gaussian approximation, the uncertainty in the total intensity measurement is given by $E_{\text{tot},r}={\sigma }_{\text{tot}}(r)/\sqrt{N_{\text{pix, aperture}}(r)}$. Therefore, dividing the elliptical aperture into 60-degree wedges reduces the number of pixels per wedge by a factor of 6 and increases the uncertainty by a factor of 2.5.

This division into 60 degree wedges still provides a sufficient S/N for reliable edge detection. Therefore, to map the edge at different angular positions (R_edge, θ), we divided the galaxy into 60-degree wedges. This means that at the semimajor axis, we defined a wedge spanning PA = −30 deg to 30 deg, then shifted by 30 degrees to define another wedge from PA = 0 deg to 60 deg, and so on. This means 12 wedges in total.

For each wedge, we extracted the radial profile using masked data and sampled it using the task astscript-radial-profile. As in the case of the 1D approach, the final value of R_edge, θ is the mean between the minimum obtained with the numerical derivation and the analytical derivation of the polynomial fit. R_edge, θ are represented as yellow dots in the upper panels of Fig. 5. The errors associated with these estimates correspond to ±5″ (the half width of the step). Fig. 6 shows the distribution of the R_edge, θ and corresponding Σ_⋆(R_edge, θ) compared to the value obtained with the 1D analysis. The 1D and 2D methods yield consistent edge values, with angular positions of R_edge, θ clustering around R_edge (few kiloparsec), capturing small asymmetries in the disk’s structure. We calculated the average absolute deviation of all R_edge, θ from R_edge, and express this value as a percentage of R_edge. This gives a mean asymmetry of ∼5% for NGC 3486. The stellar mass densities at the edge Σ_⋆(R_edge, θ) cluster around the mean edge density Σ_⋆(R_edge), ranging between 0.75 and 4.5 M_⊙/pc². In Fig. 6, R_edge, θ values are color-coded by their corresponding angular coordinate, θ, revealing that neighboring angular sectors exhibit similar R_edge, θ, suggesting local coherence in edgestructure.

Fig. 6. Distribution of stellar surface mass density at Redge, θ. The green box corresponds to 1.3 ± 0.1 M⊙/pc2 at Redge. The orange line corresponds to the median of the Σ⋆(Redge, θ). The colored region corresponds to a 1σ dispersion. Coincidence of locations at Redge, θ occurs because of the limited resolution of the method. Points are color-coded by θ.

5. Discussion

Our results indicate that it is possible to accurately detect the edge of galaxies through an unbiased mathematical definition, using the stellar surface mass density distribution. NGC 3486 shows an almost symmetric disk (asymmetry < 5%) with a sharp decline beyond R_edge = 205″ ± 5″ equal to 13.8 ± 0.6 kpc at a distance of 13.6 Mpc. This corresponds to μ_r(R_edge) = 25.60 ± 0.01 mag/arcsec², μ_g(R_edge) = 25.95 ± 0.01 mag/arcsec², a mean surface stellar mass density of 1.3 ± 0.1 M_⊙/pc², and Σ_⋆(R_edge, θ) values between 0.75 and 4.5 M_⊙/pc².

We note that, assuming the IMF does not vary with radius, switching to a different IMF (e.g., Salpeter) would uniformly scale the stellar mass density profile up or down. While this would affect the Σ_⋆(R_edge) estimate, it would not alter dramatically the overall shape of the profile. Consequently, the absolute values of Σ_⋆(R_edge) and Σ_⋆(R_edge, θ) reported here areIMF-dependent and would vary accordingly.

The results of this work therefore align with the ones of Chamba et al. (2022) for galaxies of a similar mass and redshift (M_* ∼ 10¹⁰ M_⊙). They report R_edge values ranging from 6 to 20 kpc (∼16 kpc on average) and a stellar mass density at the edge for S0/a-Irr galaxies between 0.2 and 4 M_⊙/pc², with an average of ∼1.15 M_⊙/pc² (see their Fig. 5). The galaxies in Chamba et al. (2022) that have both a similar stellar mass and an average g − r color to that of NGC 3486 (∼0.3) have average edge radii greater than 10 kpc and tend to be larger than their redder counterparts.

5.1. The formation of a sharp truncation in NGC 3486

NGC 3486 shows a broken exponential stellar surface mass density profile, characterized by a change in slope that separates the decline of the inner region from the steeper outer regions. This transition appears sharp, with a ratio of h₂/h₁ ∼ 5 after removal of the effect of the PSF. This sharp transition is similar to that observed in the edge-on galaxies NGC 5907 and NGC 4565 (Martínez-Lombilla et al. 2018). The near-ultraviolet (NUV) and far-ultraviolet (FUV) profiles (see Fig. A.1) also show a sharp drop, probably indicating negligible star formation beyond the edge. A similar break is also suggested in the infrared profile from Richards et al. (2016), where a change in slope is observed at ∼200″ (see their Fig. A.1).

The presence of this clear edge suggests an ongoing star formation process. Over time, this sharp boundary is expected to fade as stellar migration redistributes stars from the inner regions inside R_edge into the outer disk (Sellwood 2014; Debattista et al. 2017), gradually erasing signatures of past star formation. Further evidence for this star formation link can be found in the g − r color profile. Within R_edge the galaxy shows active star formation, while beyond this radius the g − r color becomes redder, possibly due to stellar migration. This is consistent with the results of Martínez-Lombilla et al. (2018), who interpret the gradual reddening beyond the edge as passive aging of the stellar population.

Under the hypothesis that R_edge marks the limit of in situ star formation, the formation and evolution of the edge is intimately linked to the abundance of HI gas. The absence of star formation in the outer disk is likely tied to the availability of cold gas and feedback effects. Observations of HI and molecular gas distributions in galaxies with similar profiles show that gas densities drop below the critical threshold for star formation beyond the edge. Using Very Large Array (VLA) data, Richards et al. (2016) found that NGC 3486 has an asymmetric HI gas distribution that extends from the northwest and southeast sides of the HI disk. They also infer a mean inclination from the HI data at 200″ of ∼40 deg that is consistent with the AR of 0.765 found in this work. In Fig. H.1 (see Appendix H) we have tentatively superimposed the HI contours from Richards et al. (2016) on the LIGHTS data. The asymmetric HI gas distribution is further analyzed by Smith et al. (2021), who reports two HI arcs extending from the southeastern and northwestern parts of NGC 3486 (see their Fig. A.9) and rotating with the rest of the galaxy. We hypothesize that the observed HI asymmetry could be due to an interaction with a satellite galaxy (PGC 033184; see Fig. B.1) located to the north of NGC 3486 at z = 0.008. The gas content of NGC 3486 was also studied by Broeils & Van Woerden (1994) using short 21-cm line observations from the Westerbork Synthesis Radio Telescope. Fig. 7 shows a comparison between their derived surface gas mass density profile (blue) and the surface stellar mass density of this work (orange). We have adapted their profile to the distance assumed in this paper. At R_edge, the surface HI gas mass density reaches ∼5 M_⊙/pc², consistent with observations of nearby spiral galaxies (Martin & Kennicutt 2001; Leroy et al. 2008). In NGC 3486, the star formation efficiency at R_edge is around 25%, as the surface gas mass density at this radius is ∼3.8 times higher than the stellar mass density. This result is therefore consistent with theoretical models suggesting that the star formation threshold in galaxies is regulated by thermal instabilities arising from the transition between warm and cold gas phases (Schaye 2004).

Fig. 7. Surface stellar (orange; this work) and gas (blue; Broeils & Van Woerden 1994) mass density profiles of NGC 3486. The stellar mass density profile shows an edge feature at 205″ ± 5″ (gray band), which coincides with the radius at which the gas profile also decays more rapidly. We assume a distance of 13.6 Mpc (Kourkchi et al. 2020). Σ⋆(Redge) is marked with a horizontal red line.

5.2. Method limitations and future directions

We evaluated the performance and implications of our new edge detection method, which mitigates the effects of scattered light and helps to identify the edge location (and therefore the characterization of the size using this feature as a proxy) in galaxies. This method is crucial in deep imaging, providing a robust approach to identifying galaxy edges based on second derivatives. This technique also provides the first wedge-based characterization of the border of galaxies. The code is modular and automated, including profile extraction, derivative calculation, and edge detection in both 1D and 2D. We acknowledge that masking LSB features, deciding on the sampling of the surface brightness and mass profiles, or adapting wedge widths still involves manual decisions and may carry subjective elements. Computing the second derivative can introduce numerical errors, which may stem from point spacing or limited floating-point precision. Our approach minimizes smoothing to preserve data authenticity, but this remains a balancing challenge. While computationally efficient overall, the most resource-intensive step remains the PSF correction of the data. Future work (Golini et al., in prep.) will explore the method’s reproducibility across datasets with varying depths, noise levels, resolution (i.e., different sampling of the mass profile), and PSF models, and expand its application to a broader range of galaxy properties (e.g., inclination, size, mass).

6. Conclusions

We have introduced a robust method of measuring the edges of galaxies by leveraging the second derivative of the stellar surface mass density profile. This approach serves as a complementary quantitative alternative to traditional visual-based methods (Martínez-Lombilla et al. 2018; Trujillo et al. 2020; Díaz-García et al. 2022; Chamba et al. 2022; Buitrago & Trujillo 2024) of identifying galaxy boundaries. Our technique enables both 1D and 2D analysis, allowing for a comprehensive characterization of galaxy outskirts, asymmetries, and structural evolution. By applying our method to deep imaging observations of NGC 3486 from the LIGHTS (Trujillo et al. 2021; Zaritsky et al. 2024) survey (μ_r, lim = 30.5 mag/arcsec² and μ_g, lim = 31.4 mag/arcsec²; 3σ in 10″  × 10″ boxes), we identified an outer edge at 205″ ± 5″, in strong agreement with the visual boundary. Our wedge-based analysis quantified the galaxy’s asymmetry at 5%, confirming the consistency between the isophotal ellipse parameters and the edge contours at different angular directions. The sharp transition at the edge occurs at a stellar surface mass density of ∼1 M_⊙/pc², reinforcing the possible link between galaxy edges and the threshold of in situ star formation (Quirk 1972; Kennicutt 1989; Schaye 2004; Roškar et al. 2008).

The computational efficiency and adaptability of this method make it well suited for integration into survey pipelines, ensuring consistent edge detection across different datasets, such as LSST (Ivezić et al. 2019), Euclid (Scaramella et al. 2022), Roman (Bailey et al. 2023), and ARRAKIHS (Guzmán 2024). In summary, this method offers a promising advance in galaxy edge detection, while aligning well with the needs of deep imaging surveys (such as background and scattered light removal). Continued development will enable its application across varied datasets, paving the way for an improved understanding of galaxy structures and the processes thatshape them.

Acknowledgments

We thank the referee for the constructive comments that helped improve the quality of the manuscript. GG acknowledges support from the PID2022-140869NB-I00 grant from the Spanish Ministry of Science and Innovation. IT acknowledges support from the ACIISI, Consejería de Economía, Conocimiento y Empleo del Gobierno de Canarias and the European Regional Development Fund (ERDF) under a grant with reference PROID2021010044 and from the State Research Agency (AEI-MCINN) of the Spanish Ministry of Science and Innovation under the grant PID2022-140869NB-I00 and IAC project P/302302, financed by the Ministry of Science and Innovation, through the State Budget and by the Canary Islands Department of Economy, Knowledge, and Employment, through the Regional Budget of the Autonomous Community. This research also acknowledge support from the European Union through: “UNDARK” and “Excellence in Galaxies – Twinning the IAC” of the EU Horizon Europe Widening Actions programmes (project numbers 101159929 and 101158446) and (MSCA EDUCADO, GA 101119830). Views and opinions expressed are however those of the author(s) only and do not necessarily reflect those of the European Union or European Research Executive Agency (REA). Neither the European Union nor the granting authority can be held responsible for them. RIS acknowledges financial support from the Spanish Ministry of Science and Innovation through the project PID2022-138896NA-C54. MM acknowledges support from grant RYC2022-036949-I financed by the MICIU/AEI/10.13039/501100011033 and by ESF+, and program Unidad de Excelencia María de Maeztu CEX2020-001058-M. SR acknowledges support from the GEELSBE2 project with reference PID2023-150393NB-I00 funded by MCIU/AEI/10.13039/501100011033 and the FSE+ and also the Consolidación Investigadora IGADLE project with reference CNS2024-154572. SR acknowledges financial support of the Department of Education, Junta de Castilla y León and FEDER Funds (Reference: CLU-2023-1-05). AAR acknowledges funding from the Agencia Estatal de Investigación del Ministerio de Ciencia, Innovación y Universidades (MCIU/AEI) under grant “Polarimetric Inference of Magnetic Fields” and the European Regional Development Fund (ERDF) with reference PID2022-136563NB-I00/10.13039/501100011033. J.R. acknowledges financial support from the Spanish Ministry of Science and Innovation through the project PID2022-138896NB-C55 Facilities: LBC-Red, LBC-Blue. Software: Astropy (Price-Whelan et al. 2018), Astrometry.net (Lang et al. 2010), Gnuastro (Akhlaghi & Ichikawa 2015), IMFIT (Erwin 2015), Matplotlib (Hunter 2007), NumPy (van der Walt et al. 2011), SciPy (Oliphant 2007; Millman & Aivazis 2011), source-extractor (Bertin & Arnouts 1996), SCAMP (Bertin 2006), SWarp (Bertin 2010).

References

Appendix A: Profiles of NGC 3486 at different wavelengths

We show in Fig. A.1 the NUV and FUV average profiles of NGC 3486 derived in Cejudo et al. (2025), and the averaged profiles using LIGHTS data in the g and r bands. In the same plot, we report the stellar mass density profile (Σ_*) of the galaxy obtained in this work (see main text for details). Σ_*, on average, is smoother than the lights profiles.

Fig. A.1. Light profiles of NGC 3486 using different filters. Optical and stellar mass density profiles are derived using LIGHTS data in this work. NUV and FUV profiles are taken from Cejudo et al. (2025).

Appendix B: Color image of LIGHTS data

Fig. B.1 shows the LIGHTS data of NGC 3486. The color image with gray-scaled background was constructed using astscript-color-faint-gray (Infante-Sainz & Akhlaghi 2024), by combining the g and r filters. The black and white background is from the g band.

Fig. B.1. Color composite image of NGC 3486 and its surroundings, constructed from Sloan g and r band data using astscript-color-faint-gray. Regions with low S/N appear white. The black and white background is from the g band. North is up and east is to the left. At a distance of 13.6 Mpc, 10′ correspond to 39.4 kpc.

Appendix C: Mask

Fig. C.1 shows the mask used in this work for the LIGHTS data of NGC 3486 created using a combination of g and r masks.

Fig. C.1. Mask of NGC 3486 used in this work. The mask is created by combining Noisechisel, Segment, unsharp masking and manual masking based on g and r band data of the LIGHTS Survey.

Appendix D: Structural parameters of Model1 and Model2

The IMFIT parameters of the simulated Model1 and Model2 share similarities but differ in key aspects. To create both models we used a BrokenExponential profile with a PA of −20° (clockwise from the X axis), an inclination of 39.6° and an inner scale length h₁ of 75 pixels. The location of the transition between inner and outer exponential declines was set to 200 pixels, and α is equal to 10 in both models. Differences arise in the outer scale length h₂, where in Model1 it is set to 30 pixels, while in Model2 to 25 pixels. Additionally, in the halo Exponential component, Model2 has central intensity value of 0.0004 counts/pix and scale length of 270 pixels.

Appendix E: Characterising the position angle and axis ratio of NGC3486

To extract the surface brightness profile of NGC 3486, we first used the ellipse routine in photutils (Jedrzejewski 1987; Bradley et al. 2019) to fit elliptical isophotes to the combined g and r band data. This method provides a representative PA and ellipticity of the outer galaxy regions. In the initial run of ellipse, all parameters varied freely. In the subsequent run, we fixed the center to the median coordinates derived in the first pass (RA (J2000) = 165.0994 deg, Dec (J2000) = 28.9751 deg)⁹. The fitting began at 2′ with linearly spaced elliptical annuli. We found an average AR of 0.765 ± 0.020 and PA of 13.0 ± 0.7° (clockwise from the X axis) between radii of 170″ and 220″. These values align with prior measurements, such as AR = 0.792 and PA = 8° (clockwise from the X axis; Smith et al. 2021) and PA = 5° (clockwise from the X axis) with an inclination of 36.8 from gas kinematics (Richards et al. 2016).

Appendix F: Tests on Wiener deconvolution

PSF, saturation and noise. To evaluate the accuracy of the Wiener-Hunt deconvolution method, we applied it to the simulated disk Model1. By comparing the deconvolved images to their original (not convolved) counterparts, we assessed how well the method retrieves the intrinsic averaged surface brightness profiles and we visually compared the background variations. The profiles were extracted using fixed ellipses for each model. Our findings are discussed below.

Fig. F.1 presents a comparative analysis of disk models, illustrating the effects of the PSF and Wiener-Hunt deconvolution under different conditions. The first column represents the base model, while subsequent columns introduce noise and saturation effects. Each row shows a different processing stage. The first row displays the original model (assuming infinite S/N), the same model with added Gaussian (σ = 0.003) noise using the task mknoise-sigma-from-mean function in Gnuastro (Noisy model), the Noisy model with inner 20″ replaced with Not a Number (NaNs) simulating saturation (Sat. Noisy model), and the original model adding higher Gaussian (σ = 0.03) noise (Noisy2 model). The Noisy model simulates an image with a limiting surface brightness of μ_lim ∼ 29.1 mag/arcsec² (3σ; 10″ × 10″), while the Noisy2 model corresponds to μ_lim∼ 26.6 mag/arcsec² (3σ; 10″ × 10″). These values are representative of typical surveys such as Stripe 82 and SDSS.

In the second row, the same models were convolved with the LIGHTS’ PSF (Sedighi et al. 2025) using the –psf option in IMFIT. The PSF declines as a power law (r^β) with an exponent β = -2.1. On the third and forth row, the panels display Wiener-Hunt deconvolution applied with the correct PSF and with a not optimal PSF model (PSF-2) with β = −2.2. In the fifth row, we show the averaged surface brightness profiles corresponding to each column, comparing the different processing steps. The colorbar indicates the surface brightness scale. Fig. F.1 demonstrates how, after Wiener-Hunt deconvolution, the images introduce correlated noise in the outskirts while preserving the edge location and, how Wiener-Hunt deconvolution, even in case of noise and saturation, can partially recover the original structure, down to the limiting depth of the image.

Comparison with wavelet-regularized deconvolution. Wiener-Hunt deconvolution is not the only method we can use to address the PSF-Stellar halo degeneracy. We compared Wiener-Hunt deconvolution with wavelet-regularized deconvolution techniques. Before delving into this comparison, we introduce the concept of wavelets as follows. The quality of the deconvolved images can be improved by performing a regularized deconvolution. If the noise term in the image formation assumption of Eq. 2 is Gaussian with zero mean and diagonal covariance matrix, the deconvolution can be performed in a Bayesian framework by optimizing the following metric:

$$\begin{array}{c}\hfill L={\displaystyle \sum _{x,y}}{\parallel g(x,y)-f(x,y)\ast h(x,y)\parallel }^2+\lambda R(f),\end{array}$$(F.1)

where ∥f∥₂ is the ℓ₂ norm¹⁰.

This solution is known as the maximum a-posteriori (MAP) estimate. The first term represents the log-likelihood of the observations given the deconvolved image, while the second term encodes prior assumptions about the image structure through the functional R. If no regularization is imposed (R = 0), the solution reduces to the well-known maximum likelihood (ML) estimate, which is highly sensitive to noise. The prior acts as a regularization term, incorporating prior knowledge about the image to stabilize the solution. For example, Tikhonov regularization uses R(f) = ∥f∥². In this application, we build on the ideas of sparsity-based regularization to deconvolve the images. In this case, the regularization term is given by the ℓ₁ norm of a transformed version of the image:

$$\begin{array}{c}\hfill R(f)={\parallel W(f)\parallel }_1,\end{array}$$(F.2)

where W(f) represents the isotropic undecimated wavelet transform (IUWT), also known as the starlet transform (Starck & Murtagh 1994; Starck et al. 2006). The IUWT is a redundant wavelet transform that decomposes an image into multiple frequency components, providing a sparse representation well-suited for astronomical images. For a more detailed explanation of the mathematical foundations behind wavelet-regularized deconvolution methods, we refer the reader to the following works: Starck & Murtagh (1994, 2007), Starck et al. (2010), Carrillo et al. (2012). Fig. F.2 shows a cropped region of the outskirts of NGC 3486 (3′ × 3′) before and after deconvolution and compares the averaged surface brightness profiles obtained with Wiener-Hunt and wavelet-regularized deconvolution methods using fixed elliptical apertures. We find close agreement up to 280″ corresponding to μ_g =,29 mag/arcsec².

Fig. F.1. Effectiveness of Wiener-Hunt deconvolution on a simulated disk in case of inaccurate PSF modeling, saturation and background noise (see main text for details).

Fig. F.2. Effect of different deconvolution methods on the background of LIGHTS data. Gray-scaled images correspond to the g band. The bottom right panel shows the averaged surface brightness profile of NGC 3486 after applying Wiener-Hunt (blue) and wavelet-regularized (red) deconvolution techniques, compared to the original profile (green).

Appendix G: Stability of R_edge with varying α

Figure G.1 shows the BrokenExponential profile in arbitrary units (top panel) and its second derivative (bottom panel) for two disks sharing the same PA, AR, central intensity, inner and outer scale length (h₁, h₂), and R_edge but different α. While the location of the minimum of the second derivative remains unchanged since R_edge does not vary, for fixed h₂ a smoother transition (lower α) leads to visually place R_edge at inner locations.

Fig. G.1. Derivatives of BrokenExponential disks with different α parameters.

Appendix H: HI isocontours

Fig. H.1 shows the HI contours from Richards et al. (2016) on top of the gray-scaled LIGHTS data of NGC 3486.

Fig. H.1. Low spatial resolution HI integrated intensity contours from new VLA data from Richards et al. (2016) overlapped with gray scaled image of LIGHTS data of NGC 3486.

All Figures

Fig. 1. Edge detection method. The figure presents the stellar mass density map in Cartesian and polar coordinates (top left and top central panels), and the averaged surface stellar mass density profile (bottom left panel) of a simulated disk, computed in three different cases as is explained in the main text. The first derivative of the stellar mass density profile (bottom central panel) and the second derivative in both 2D (upper right) and 1D (bottom right) are also shown. Gray lines represent the surface stellar mass density and its derivatives using 10° wedges along azimuthal angles. Dashed red lines indicate Redge, corresponding to the minimum of the second derivative.

In the text

	Fig. 2. Flowchart illustrating the steps used to account for the effect of the PSF on the outermost regions of NGC 3486 (see main text for details). The panels display LIGHTS data of NGC 3486. Each panel covers a FoV of 15.5′ × 15.5′.
In the text

	Fig. 3. Surface stellar mass density profiles of two simulated disks. Left: model1, consisting of a broken exponential disk only. Center: model2, similar to Model1, but with a steeper outer disk slope (h2) and a faint stellar halo. Right: comparison of the PSF-convolved profiles for both models.
In the text

Fig. 4. Surface brightness profiles of NGC 3486 in the g band (upper panel) and r band (central panel) using LIGHTS data. The surface brightness profiles shown correspond to the background-subtracted data (green), the Wiener-Hunt deconvolved data (blue), and the data for which the PSF effect has been removed (orange). These profiles were obtained using a fixed ellipticity and PA representative of the outer parts of the galaxy (see the main text for details). Only values above 3σ the noise of the data are shown. The pink lines correspond to the IMFIT models we used to model the Wiener-Hunt deconvolved light distribution of the galaxy. The vertical dashed black lines highlight the location of Redge with the corresponding shaded gray uncertainty, dashed vertical purple lines highlight the location of RPSFeff (the radius where the PSF starts to be noticeable), and horizontal blue lines represent the data surface brightness limiting depth. The g − r average radial profiles prior and after deconvolution are shown in the bottom panel.

In the text

Fig. 5. Edge detection in NGC 3486. Upper row: stellar mass surface density map in Cartesian and polar coordinates (left and central panels) and polar map of its second derivative (right). Bottom row: radial profile of averaged stellar mass surface density as a function of radius (central panel) along with its second derivative (right panel). A polynomial fit to the data and its second derivative are also shown (blue). In all panels shaded regions indicate measurement uncertainties. The outer boundary, Redge is indicated in all panels by a dashed black line. Σ*(Redge) is shown as a horizontal red line in the bottom central panel. Yellow dots in the upper panels mark Redge, θ.

In the text

	Fig. 6. Distribution of stellar surface mass density at Redge, θ. The green box corresponds to 1.3 ± 0.1 M⊙/pc2 at Redge. The orange line corresponds to the median of the Σ⋆(Redge, θ). The colored region corresponds to a 1σ dispersion. Coincidence of locations at Redge, θ occurs because of the limited resolution of the method. Points are color-coded by θ.
In the text

	Fig. 7. Surface stellar (orange; this work) and gas (blue; Broeils & Van Woerden 1994) mass density profiles of NGC 3486. The stellar mass density profile shows an edge feature at 205″ ± 5″ (gray band), which coincides with the radius at which the gas profile also decays more rapidly. We assume a distance of 13.6 Mpc (Kourkchi et al. 2020). Σ⋆(Redge) is marked with a horizontal red line.
In the text

	Fig. A.1. Light profiles of NGC 3486 using different filters. Optical and stellar mass density profiles are derived using LIGHTS data in this work. NUV and FUV profiles are taken from Cejudo et al. (2025).
In the text

	Fig. B.1. Color composite image of NGC 3486 and its surroundings, constructed from Sloan g and r band data using astscript-color-faint-gray. Regions with low S/N appear white. The black and white background is from the g band. North is up and east is to the left. At a distance of 13.6 Mpc, 10′ correspond to 39.4 kpc.
In the text

	Fig. C.1. Mask of NGC 3486 used in this work. The mask is created by combining Noisechisel, Segment, unsharp masking and manual masking based on g and r band data of the LIGHTS Survey.
In the text

	Fig. F.1. Effectiveness of Wiener-Hunt deconvolution on a simulated disk in case of inaccurate PSF modeling, saturation and background noise (see main text for details).
In the text

	Fig. F.2. Effect of different deconvolution methods on the background of LIGHTS data. Gray-scaled images correspond to the g band. The bottom right panel shows the averaged surface brightness profile of NGC 3486 after applying Wiener-Hunt (blue) and wavelet-regularized (red) deconvolution techniques, compared to the original profile (green).
In the text

	Fig. G.1. Derivatives of BrokenExponential disks with different α parameters.
In the text

	Fig. H.1. Low spatial resolution HI integrated intensity contours from new VLA data from Richards et al. (2016) overlapped with gray scaled image of LIGHTS data of NGC 3486.
In the text