Ultra-deep imaging of NGC1052-DF2 and NGC1052-DF4 to unravel their origins

,


Introduction
NGC 1052-DF2 1 and NGC 1052-DF4 (van Dokkum et al.  2018a, 2019), two low surface brightness galaxies ( µ V eff > 25 mag arcsec −2 ) in the line of sight of the NGC 1052 group 2 , have attracted the attention of the extragalactic community in recent years.Both galaxies have large effective radii (R eff > 1.5 kpc) for their typical stellar masses (∼10 8 M ; van Dokkum et al. 2018aDokkum et al. , 2019)).Interestingly, the very low velocity dispersion of their globular cluster (GC) systems (4−10 km s −1 ) suggests that the amount of dark matter they contain is either very small or negligible.The possibility of a low dark matter content sparked a lively debate, much of which focused on how best to estimate the intrinsic velocity dispersion of NGC 1052-DF2 (see, e.g., Martin et al. 2018;van Dokkum et al. 2018b;Fensch et al. 2019;Emsellem et al. 2019), the distance to these objects (see, e.g., Trujillo et al. 2019;Monelli & Trujillo 2019;Zonoozi et al. 2021), and/or the possibility that these galaxies are rotating (see, e.g., Lewis et al. 2020).While the above concerns can mitigate the problem of the absence of dark matter, especially in the case of NGC 1052-DF4, where the velocity dispersion of the GC system is extremely small, the absence of dark matter remains a strong hypothesis.
Numerous hypotheses have been proposed to explain the peculiar properties of these galaxies, but none has provided a definitive answer (see, e.g., Bennet et al. 2018;Famaey et al. 2018;Kroupa et al. 2018;Fensch et al. 2019;van Dokkum et al. 2022).Among the proposed formation mechanisms, the most compelling explanation is that these galaxies had a typical dark matter content, which they subsequently lost through interactions with neighboring objects (Ogiya 2018;Jackson et al. 2020;Macciò et al. 2021;Ogiya et al. 2022).
Since the distribution of dark matter is more extended than the stars, during a tidal interaction with a more massive neighbor, the galaxy first loses its dark matter and then its stars begin to be stripped along two symmetrical tidal tails emanating from two opposite sides of the galaxy.These tails exhibit an S-shaped structure with a width equal to the diameter of the disrupted object (e.g., Johnston et al. 2002;Moreno et al. 2022).Recent simulations (see, e.g., Moreno et al. 2022;Katayama & Nagamine 2023) show that if the dark matter in NGC 1052-DF2 and NGC 1052-DF4 has been stripped by gravitational interaction with a more massive object, then the tidal tails should be visible in ultra-deep imaging, reaching a surface brightness of µ V ∼ 30 mag arcsec −2 (even if the interaction took place several gigayears ago).
Therefore, the key to testing these tidal stripping scenarios is to have very deep imaging of the galaxies under study.Several papers in the past have examined these two galaxies with deep data to analyze their outskirts and surroundings, looking for (among other things) tidal features (see, e.g., Müller et al. 2019;Trujillo et al. 2019;Román et al. 2021).The tidal stripping scenario for NGC 1052-DF4 has recently been confirmed observationally by Montes et al. (2020) using the IAC80 telescope and by Keim et al. (2022) with deep imaging from the Dragonfly Telephoto Array.Both studies have independently shown, using very deep imaging, that NGC 1052-DF4 has the characteristic S-shaped structure indicative of tidal interaction.While the two teams agree on the tidal distortion of the galaxy, the massive galaxy causing the distortion is a matter of debate.Montes et al. (2020) argue that the observed tidal tails point to NGC 1035, a medium-sized spiral very close to NGC 1052-DF4 in projection, which also shows signatures of being distorted.On the other hand, Keim et al. (2022) suggest that the massive galaxy producing the tidal features of NGC 1052-DF4 is NGC 1052.
For the other low-mass galaxy, NGC 1052-DF2, the situation is far from clear.Montes et al. (2021), using deep images taken with the Isaac Newton Telescope (INT), find no evidence for tidal tails in the outer parts of the galaxy.Instead, they find a surface brightness profile that is cored in the inner region and decreases exponentially in the outer part.They argue this could indicate the presence of a low-inclination disk, supporting the evidence for rotation suggested by Lewis et al. (2020).However, Keim et al. (2022) interpret this core-like shape in the innermost region of the galaxy as a signature of the tidal interaction of NGC 1052-DF2 with NGC 1052.
While the images obtained so far for the two galaxies are very deep (µ V ∼ 29−29.5 mag arcsec −2 ; 3σ in 10 × 10 boxes), they seem insufficient to determine the presence of tidal disruption features.The aim of the present work is to probe the galaxies with images obtained by the Gemini telescopes that are about 1 mag deeper than previous data published.Following the motivation of Moreno et al. (2022) and Katayama & Nagamine (2023), we hope that with these ultra-deep data the presence of tidal distortions, if present, can be unambiguously revealed.
The paper is organized as follows.Section 2 provides an overview of the observational strategy employed and outlines the data reduction pipeline used to obtain our ultra-deep images.Section 3 presents the results of our analysis, including the surface brightness profiles of NGC 1052-DF2 and NGC 1052-DF4.In Sect. 4 we compare our results with previous research.Finally, we summarize our conclusions in Sect. 5.All magnitudes in this paper are given in the AB magnitude system.

Data
The ultra-deep images of NGC 1052-DF4 and NGC 1052-DF2 were obtained with the two Gemini Multi-Object Spectrographs (GMOS) 3 mounted on the Gemini North and South telescopes (Hook et al. 2004).The GMOS cameras cover a 5 × 5 square area with truncated corners.Each camera consists of three adjacent charge-coupled device (CCD) sensors, and each of the three CCDs is read by four amplifiers.The first two amplifiers of CCD1 and the last two of CCD3 are outside the imaging field of view (FoV) and, therefore, they do not contribute to the final dataset.On the new Hamamatsu CCDs for GMOS-S, the gaps between the CCDs are 4.88 wide (61 pixels in unbinned mode), while the gaps between the CCDs for the GMOS-N Hamamatsu CCDs are 5.4 wide (67 pixels in unbinned mode).We worked in 4 × 4 pixels binned mode resulting in a final pixel scale of 0.32 pix −1 .Due to the specific configuration of the telescope, vignetting from the On-Instrument WaveFront Sensor (OIWFS), which provides guidance and tip-tilt correction information to the telescope control system, appears as an artifact on the data, which is properly accounted for as described in the following Sect.2.1.
3 http://www.Gemini.edu/instrumentation/gmosNGC 1052-DF2 was observed on October 28, 2021, with GMOS-S (r filter) and on November 8, 2021, with GMOS-N (g filter).The data of NGC 1052-DF4 (r filter) were obtained on October 31, 2021, with GMOS-S.The mean seeing during the observations was 1.3 in r filter for DF4 and 0.88 and 1.2 in r and g filters, respectively, for DF2.
To perform photometric calibration of Gemini data and check the consistency of our results, we used public Sloan Digital Sky Survey (SDSS) and Dark Energy Camera Legacy Survey (DECaLS) data that are currently available for the objects.We retrieved the SDSS gand r-band imaging data from the DR14 SDSS (Abolfathi et al. 2018) Sky Server.The magnitude zero point for the SDSS dataset is 22.5 mag and the pixel scale is 0.396 pix −1 .We downloaded gand r-band optical data from DECaLS (Proposal ID 2014B-0404; PIs: David Schlegel and Arjun Dey) archive with a zero point of 22.5 mag and a pixel scale of 0.262 pix −1 .We used the DECaLS DR10 data from bricks 0405m085 and 0405m083 for NGC 1052-DF2 and NGC 1052-DF4 fields, respectively.

Observational strategy
To reach the depths we are aiming for in this work (i.e., limiting surface brightness of around µ g ∼ 31.0 mag arcsec −2 , 3σ in 10 × 10 ), we need to deal with several observational biases that can affect deep imaging, such as scattered light, saturation, or ghosts.In addition, we need an accurate estimation of the flat field and a careful treatment of the sky background.To this end, we designed a special observing strategy for the Gemini telescopes that consists of a dithering pattern with offsets of 10 .The offset is needed to fill the gap between the CCDs in the cameras and to minimize scattered light from the telescope optics.This offset also ensures that we minimize the errors introduced by quantum efficiency and differences between the amplifiers, since the galaxy always overlaps on the central CCD (CCD2).The final dithering pattern around both galaxies is shown in Figs.A.1 and A.2. NGC 1052-DF2 is not centered in the Gemini FoV to better characterize the scattered light from the bright star to the east (see Sect. 2.2.7).
The differences between the night sky and the twilight illumination can lead to subtle differences in the flat field.The twilight light gradient can introduce unwanted gradients in our data.Therefore, the twilight flats provided by the observatory are unsuitable for our purpose.A more accurate flat-field correction is obtained using the same scientific exposures obtained with a large displacement dithering pattern (typically of the size of the object under observation; see Trujillo & Fliri 2016).
Unfortunately, the main limitation of using this approach with the Gemini telescopes revolves around the size of the OIWFS patrol field.When making large offsets (>10 ), the risk is to push the guide star beyond the patrol field's boundaries.Consequently, employing the same science field with extensive offsets could lead to loss of guiding, potential vignetting in one or more images, and a reduction in the effective FoV, resulting in a smaller area with a higher signal when stacking images.This is true for NGC 1052-DF4, where the Guide star is located in the external border of the OIWFS patrol field while in the case of NGC 1052-DF2, the Guide star was closer to the GMOS FoV.
To mitigate this problem, every six science observations of the galaxy, the telescope is pointed at a nearby empty region of the sky (∼7 away from the galaxy).In such a region, we took five unguided exposures on sky.With this sky data we created the flat field.Using this technique, we retrieved ∼20 min of sky observations for each night to build the flat field from the data.A99, page 2 of 21 Notes.The table shows the exposure time and the number of frames taken during each night for both the sky and the on-source images.The surface brightness limit of the data (3σ; 10 × 10 boxes) and the camera used are also given.
At the end of the observing run, we got 1 h 30 min for NGC 1052-DF4 (r filter), 1 h 20 min for NGC 1052-DF2 (r filter), and 1 h 30 min for NGC 1052-DF2 (g filter).Table 1 lists the exposure times for each filter and the number of frames, together with the limiting surface brightness depths of the data.Bias frames were taken at the end of the nights as part of the observing programs.They consist of 20 bias frames per CCD.

Frame processing
In this subsection, we explain in detail how the sky images and the science images are processed and assembled to obtain the final co-addition.The main steps of the process are: (i) CCD assemblage (gain correction), (ii) the derivation of the flat field using sky images, (iii) custom processing of galaxy images, and (iv) mosaic co-addition.
The pipeline makes effective use of the Gnuastro Software (Akhlaghi & Ichikawa 2015) and its tasks.

Bias and gain correction
First, for each of the three CCDs and for each individual image (both on-source and on-sky), we generated a bad pixel mask.This mask account for various cosmetic defects present on the CCDs, such as the readout failure of the detector, over-scan regions, and bad columns.In practice this mask implies changing these pixels to not-a-number values.
Also, we masked all pixels with values greater than 5 × 10 4 ADUs to account for the sensor nonlinearity.After this, we produced bias-corrected images.To do so, we created a master bias (for each CCD) by combining the 20 bias frames we have for each of the cameras and nights with a resistant mean (using a 3σ rejection).
After subtracting the master bias, we flattened the response of the CCDs, since they are an array of four different amplifiers, each with a different gain.For this reason, we used the following strategy to obtain the most homogeneous image possible.We started by selecting one of the central amplifiers as a reference for the other amplifiers.We masked all the sources of the images to avoid their light contamination.Also, we assumed continuity of the background between two adjacent amplifiers.This is a reasonable assumption, since we do not expect strong changes in the illumination of the camera, nor background relevant gradients between adjacent pixels on the sky.In our case, the background pixels of the reference amplifier are selected from a column at the edge with a width of 10 pixels.The median of these pixels is our reference background for that amplifier.We did the same estimation for the adjacent amplifier using the adjacent column.We then calculated the ratio between the two quantities to multiply the gain we want to use with respect to the reference gain.Once these two amplifiers are corrected, the same procedure is applied to the other amplifiers of the CCD.While the above procedure nicely flatten the gain levels at the boundaries of the amplifiers for each individual image, it still leave some doubt about how good our initial hypothesis is about the continuity of camera illumination and background light between adjacent amplifiers.For this reason, we tried to minimize any potential deviation from this hypothesis from frame to frame by applying a final gain correction that has the information of all individual frames of a given night.We assumed that the gain of the amplifier remains constant during the few hours of observation of the source and the sky.In this case, the correction that we finally applied to leverage out the amplifiers is the median of all the corrections between any two amplifiers that we find from one frame to another.We find that the ratio between the backgrounds applied to do the flattening correction is extremely stable over the observation run, with a variability around 0.1%.

Flat-field correction
In this section, we explain how we combined our on-sky images to build the flat field in a two-step process, following a strategy similar to that described in Trujillo et al. (2021).Once the raw sky images are bias-and gain-corrected (as described in Sect.2.2.1), we computed a crude flat-field image using a median of the normalized sky images for each CCD, without any masking.In order to follow the illumination pattern of the entire GMOS camera, the normalization of each individual CCD that makes up the camera does not use the pixels of the entire CCD, but only the portion of the CCD that is expected to be illuminated with the same efficiency.In practice, this means that each CCD image is normalized using the area inside a circular annulus of 100 pixels width centered on the central pixel of CCD2 (i.e., the central CCD of the camera) and radius of 420 pixels.With this radius, we make sure that the annulus covers the entire FoV of the camera (CCD1 + CCD2 + CCD3).Then, for each CCD, we calculated the resistant median of the pixels within the corresponding annulus section.This value was used to normalize the flux in each frame before they were added together to create the flat field.
The rough flat fields are the combination using the sigmaclipped (3σ) median of the normalized sky images.Then, to get an improved master flat field, we used the corrected (by the first rough flat field) sky CCDs images to build a mask of the astronomical sources using NoiseChisel (Akhlaghi 2019).This was done because even with a first rough flat-field correction the images show the astronomical sources very well.With these pixels masked, we next built a much more accurate flat field, following the previous steps again.
Having the master flat field generated, each bias-and gaincorrected science frame was divided by the master flat field.These master flat-field-corrected images were then used to generate better masks of the astronomical sources of the science images.These masks were therefore used to improve the A99, page 3 of 21 previous step related to the gain correction.We, consequently, repeated the previous steps to improve our reduction.With the very final version of the master flat field, we performed the flatfield correction of our science frames.
Finally, to avoid vignetting problems toward the corner of the detectors, we removed from the images all the pixels where the master flat field have an illumination lower than 90% of the central pixels.

Astrometry
In the following, we use the entire assembled frame of the camera (i.e., CCD1 + gap + CCD2 + gap + CCD3) as the unit of work.To construct the complete pointing of the camera, we took the header of the preprocessed images with THELI (Schirmer 2013), a data processing pipeline for astronomical images, to have a first astrometric solution.We then used the makenew operator in Gnuastro to create the gap image between the CCDs (with different widths depending on the camera) and we stitched the CCDs and gaps together in a single fits file with a single extension.On the images we ran SExtractor (v.2.25.2;Bertin & Arnouts 1996) to generate the input catalogs for SCAMP (v.2.10.0;Bertin 2006) to improve the astrometry and the distortion of the camera.SCAMP was run twice, once with MOSAIC TYPE4 set to LOOSE and once in UNCHANGED mode.In the initial iteration, the CCD images are considered as separate and mechanically unconnected devices.Pattern matching is conducted independently for each extension.In the subsequent iteration, no preprocessing is assumed and the new distortion headers are computed on the basis of the first iteration.This approach allows for a more accurate characterization of distortions, achieving a precision of 0.01 .These new headers are our final world coordinate system (WCS) solutions for each pointing.

Sky subtraction
Sky subtraction is a crucial step in reducing low surface brightness data; if not done correctly, it can introduce artificial gradients and affect the faint objects under study.In this work, since the FoV of the camera is not very large, we worked under the hypothesis that the sky in each frame can be well described by a single constant.
To compute the sky, we first masked the brighter sources and the diffuse light with Noisechisel.To account for potential non-detections by Noisechisel of diffuse light surrounding bright objects, a manual mask is applied after the initial Noisechisel detection process.
In the case of the NGC 1052-DF2 field, we added a manual circular mask around the bright star (i.e., HD 16873 with m V = 8.34 mag and G8III, located at RA = 40.5395deg and Dec = −8.3768deg; see Fig. A.1) close to our FoV.The manual mask has a radius of 4 .This radius was chosen so that at this radial distance the expected contribution of scattered light from this bright star is ∼28.5 mag arcsec −2 .This is roughly the limiting surface brightness depth of each of our individual exposures.Consequently, we used such masks to try to minimize the effect of the scattered light on our determination of the sky background.For the NGC 1052-DF4 field instead, we built an extra elliptical mask around the galaxy NGC 1035 with a position angle (PA) of 20 • (clockwise from the north axis), an axis ratio of 0.28, and a maximum major axis of 4 .In the same field, we also added a circular mask with a radius of 1 over the nearest bright star (RA = 39.7778;Dec = −8.1121)toward the west side of the galaxy.In both fields, we computed the sky values of each frame as the 3σ clipped median of the non-masked pixels and subtract these values from the data.Before removing the sky value, we visually inspected Noisechisel's output on the distribution of background pixels within the image to ensure that our assumption of a constant value is accurate.We performed a visual inspection to identify any remaining gradients in the image.No gradients were observed.

Photometry
Once the sky value has been removed from each individual exposure, we converted the Gemini pixel values from analog digital units (ADUs) into into linear physical units of nanomaggies5 .To do this, we used the SDSS data covering the same FoV around our target.For both the Gemini images and the SDSS mosaic, we stored a catalog of the flux of the sources in the field using circular apertures with a diameter of 2 with Gnuastro astmkcatalog.These catalogs were then checked against the Gaia DR3 star catalog to make sure that we were using verified point-like sources for the photometry.We stored only the parameters of the stars that are not saturated in the Gemini data (i.e., we used only stars with magnitudes g > 18 and r > 17.8 mag) in order to minimize errors in the estimation of the flux of the sources.In addition, we only used stars with high enough signal-to-noise in SDSS images, that is, g and r magnitudes of less than 21 mag.For all unsaturated stars, we computed the ratio of the flux in the SDSS circular aperture (F SDSS,ap ) to the flux in the Gemini aperture (F Gem,ap ) for each science exposure.Then we multiplied the Gemini image (im G ) by the resistant median of all these ratios.At the end of the process we have the pixel values of each pointing in units of nanomaggies (im G,n ).

Image co-addition
In the final step of data reduction, the individual exposures are added together to form the final mosaic.The following steps are motivated by the fact that the conditions of our observations change during the night.For example, the airmass increases or decreases depending on the trajectory of the target in the sky as it moves further or closer to zenith (see Fig. A.3).In addition, other meteorological phenomena such as the passage of a cloud can change the brightness of the sky value.These local variations manifest themselves as a degradation of sky quality: the higher the standard deviation of the (photometrically calibrated) sky background pixel values, the worse the image quality.To account for these effects and improve the accuracy of the low surface brightness features, the images taken under better conditions should weigh more than the images taken under worse conditions.
To carry out the co-addition of the individual frames, we used a weighted average to stack the images, where the weight assigned to the image is the ratio of the standard deviation of the sky of the best exposure to the standard deviation of the sky of the ith frame.However, before stacking the data in this way, it is necessary to mask those pixels that have been affected by unwanted signals such as cosmic rays.These pixels are identified using the following procedure.First we used Swarp (Bertin 2010) to put the images into the same astrometric solution.Then

Original
Model Residual NGC1052-DF2 we created a mask to reject the pixels that have a value that deviates from the resistant mean (3σ rejection is sufficient for our purpose) of all pixels at the same WCS location, in all frames.
We then masked the outliers and stacked the data using this sigclip-weighted mean.
The co-added image is significantly deeper than any single image, and therefore a large number of very low surface brightness features, previously invisible, emerge from the noise.These features subtly affect the sky determination of our individual science images, and for this reason it is necessary to mask these regions and repeat the entire sky determination process on the individual exposures with an improved mask.In short, we repeated the sky estimation (and subsequent reduction steps) described above on the individual images using the improved masks generated by this first data co-addition.

Scattered light removal in the NGC 1052-DF2 field
As mentioned above, one of the main difficulties in obtaining a high quality image for low surface brightness science in the field around NGC 1052-DF2 is the presence of a very bright star (HD 16873) close to the galaxy.Therefore, we needed to develop a specific strategy to remove the light gradient produced by the scattered light of such a bright source outside the FoV of our co-added image.Figure 1 illustrates the three stages of subtraction of the scattered light emitted by the star HD 16873 near NGC 1052-DF2.The left panel of the Fig. 1 shows that the light from the star HD 16873 manifests itself as a clear excess of light on the northeast side of the final reduction of NGC 1052-DF2.
To remove the light pollution from the bright star, we did the following.First, we masked all major bright sources in the final Gemini mosaic with Noisechisel.Second, only by using the non masked pixels in the Gemini stacked image, we obtained the radial profile of the scattered light produced by that star using circular apertures centered on the star position (i.e., RA = 40.5395deg, Dec = −8.3768deg).In practice, we get the light profile of the star only in the Gemini FoV region (from ∼2.4 to ∼9 ).Such a profile is fitted with a power law shape βr α .We find that the slope of the light due to the scattered light from the star is well described by a power law with α = −2.6 for the g filter and α = −2.75 for the r filter.We then converted the power-law radial profile into a 2D circular model (middle panel of Fig. 1) centered on the star position and finally subtracted it from the original image.The right panel of Fig. 1 shows the Gemini data after removal of the scattered light model.

Scattered light removal around NGC 1052-DF4
In the case of NGC 1052-DF4, there are no bright stars (m V < 9 mag) in close proximity (d < 5 ) to the galaxy as in the case of NGC 1052-DF2.Nevertheless, there are a few bright objects that could introduce unwanted scattered light and complicate the analysis of the outskirts of this galaxy.One is the galaxy NGC 1035 to the east (see Fig . A.2).The other is the star at RA = 39.7779 and Dec = −8.1126, to the west of NGC 1052-DF4.This star has a brightness of Gmag = 14 mag.
To subtract both the nearby star and the galaxy, we created light models using build-ellipse-model in Photutils.The advantage of this method is the ability to change the PA and ellipticity of the isophotes at different radii to create more realistic models and minimize residuals.To obtain a satisfactory fit of NGC 1035, whose physical center is outside the FoV of the Gemini camera, we used a technique using Gnuastro's astwarp to combine the Gemini data with surrounding INT (Montes et al. 2021) data.To prepare the data for modeling, we also applied a mask based on the Noisechisel detection map.Then, by deriving the optimal models for the star and galaxy, we performed the subtraction of these models (middle panel of Fig. 2) from the original data.The right panel of Fig. 2 shows the Gemini data after removal of the scattered light models.

Surface brightness limits of the final co-adds
The NGC 1052-DF4 and NGC 1052-DF2 mosaics are generated using the data reduction pipeline described above.Each field has a different scattered light distribution that is processed individually.Once this was done, we estimated the surface brightness depth of our images.The surface brightness limits we obtained using the metric 3σ fluctuations in 10 × 10 boxes are 30.9mag arcsec −2 (g band) and ∼30.5 mag arcsec −2 (r band; see Appendix A of Román et al. 2020 for more details on the implementation of the chosen metric).To put our new images in context with previous datasets, in Table 2 we show the limiting surface brightness of earlier images using the same metric as in the work presented here.With the exception of the g-band data A99, page 5 of 21

Analysis
In this section we use our ultra-deep images to check for the presence of tidal features and to analyze the light distribution around NGC 1052-DF2 and NGC 1052-DF4.

NGC 1052-DF2
To put the depth of the current dataset into context, Fig. 3 shows a region of 3 × 3 around NGC 1052-DF2 using different datasets, namely SDSS, DECaLS (DR10), INT, and Gemini.The limiting surface brightness of each dataset is shown in the figure.The color of the images is constructed using astscriptcolor-faint-gray (Infante-Sainz & Akhlaghi 2024), by combining the g and r filters, with the white and black background corresponding to the g filter.The apparent size of the galaxy increases with image depth.

Radial surface brightness profiles
To derive the surface brightness profiles of NGC 1052-DF2 in both bands, we first created a common mask for both filters using an image resulting from the sum of the g and r bands.Once all sources of contamination are masked out, we ran ellipse (Jedrzejewski 1987) in Photutils (Bradley et al. 2019) to obtain the best PA and ellipticity of the isophotes in the outer parts of the galaxy.This code fits elliptical isophotes to the 2D images of galaxies using the method described in Jedrzejewski (1987).We first ran ellipse allowing all parameters to vary freely.In the second run, we fixed the center to the median centers of the isophotes returned by ellipse in the first iteration 6 .ellipse fits the isophotes both outward and inward from an initial semimajor axis (SMA) length.In both cases, we In Fig. 4 we make a comparison of the ellipticity and PA radial profiles derived from our data with those found in the literature.This analysis is relevant because changes in ellipticity and PA with radius could indicate changes in the structure of the galaxy.These changes could be due either to internal structures (triaxiality, spiral arms, etc.) or to external effects (tidal distortions, minor mergers, etc.).To make such a comparison, we used the combined g + r image8 .This was done to provide a consistent comparison with other analyses in previous papers that used this combined information.There is generally good agreement (within the error bars) between the different datasets.between 40 and 80 (green horizontal lines in Fig. 4).We find an axis ratio b/a = 0.83 ± 0.01 and PA = 51.1 ± 2.7 deg.Based on the elliptical and PA radial profiles, NGC 1052-DF2 appears to have two separate structural components, before and after a radius of ∼40 .Since the image of the galaxy presented in Fig. 3 shows no obvious signature of distortion, we favor the hypothesis that this structural change is most likely related to an internal change in the morphology of the galaxy.We discuss this idea in Sect. 4.
Before extracting the surface brightness profile, we needed to correct our images for any residual sky contribution near the galaxy.To do this, we computed the value of the background level of the data in a masked elliptical annulus centered on the galaxy.This elliptical annulus has the shape that is representative of the outer part of the galaxy.That is, it has the mean PA and mean axis ratio values we derived above.The annulus has inner and outer radii of 120 and 140 , respectively.
Once the local sky value is subtracted from the image, we extracted the radial surface brightness profiles of the g and r filters using astscript-radial-profile (Infante-Sainz et al 2024).To get a better representation of the outer part of the galaxy, we set the ellipticity and PA to the average values we found above.This is particularly useful if we want to extract information beyond R = 80 where the uncertainties in measuring the ellipticity and PA of the outer isophotes increase significantly.Using the same ellipticity and PA values, we extracted the surface brightness profiles of the SDSS data and DECaLS data.
For consistency, the mask that we applied to these shallower datasets is the one derived for the Gemini data, but scaled to the respective pixel sizes (0.396 pix −1 for SDSS and 0.262 pix −1 for DECaLS).The results of doing this are shown in Fig. 5.
The errors in the surface brightness profiles are calculated taking into account the two main sources of uncertainty in the estimate of the intensity at a given radial distance.Such an intensity is the sum of the intensity of the source itself and the intensity of the background.In short, we want to calculate the error associated with the sum I tot (r) = I source (r) + I background using the pixels located at a given radial distance within an elliptical annulus aperture.
The intensity of the background in ground-based optical dark observations, as in our case, is mainly produced by the brightness of the sky atmosphere (µ V ∼ 22 mag arcsec −2 ).There are other sources of background light that could be relevant, such as zodiacal light and the scattered light produced by nearby bright sources.In the particular case of our observations, these other background sources are subdominant with respect to the brightness of the night sky emission.For example, the scattered light from the nearby star HD 16873 is contributing a surface brightness of µ V > 28 mag arcsec −2 to the surroundings of NGC 1052-DF2.
For both the source and background intensities, the error is characterized by the shot noise, and therefore the uncertainty per pixel for each component could be written as √ I source and I background .In practice, since we are dealing with many pixels within the elliptical aperture, the characterization of the noise of each component can be done with good accuracy by measuring the standard deviation of the counts of each component within the elliptical aperture using only the pixels that are not masked (i.e., N pix,ellipse (r)).Thus, when measuring the intensity profile of the object, we can estimate σ tot (r) by calculating the standard deviation of the counts within such an elliptical aperture.Under the Gaussian approximation, the uncertainty E tot,r when measuring I tot is therefore To estimate the error associated with the determination of the background intensity, we can use the background pixels of the image to accurately determine σ background , which, assuming the sky background intensity is a constant, is independent of the radial position of the object isophote we are analyzing.However, the uncertainty E background,r in characterizing the contribution of I background to I tot (r) depends on the number of pixels within the elliptical aperture.For this reason, we estimated the E background,r when measuring the contribution of I background at each radial distance as Having measured the uncertainty at E tot,r and E background,r , we can then estimate the uncertainty associated with the source itself, E source,r : In Fig. 5 we compare the surface brightness profiles in both filters (and their combination) obtained from our ultradeep Gemini images with those derived using the SDSS and A99, page 8 of 21 DECaLS datasets.Here, the surface brightness profiles correspond to those derived using fixed ellipticity and PA.There is a notable change in the surface brightness slope around 60 .This sudden decrease (also called truncation) has been observed in other dwarf galaxies with similar stellar masses (Chamba et al. 2022).The fact that the truncation is visible in both bands simultaneously, and also in the DECaLS data, gives us confidence that this feature is likely real.The surface brightness profiles when the ellipticity and PA are allowed to vary are shown in Finally, taking advantage of our very deep dataset, we derived the effective radius using the surface brightness profiles in the two different filters.By integrating the light profiles extracted using fixed ellipses, we obtain R eff,g = 25.8 ± 0.07 and R eff,r = 25.6 ± 0.07 .If the ellipticity and PA of the ellipses from which the surface brightness profiles are obtained are left free, the effective radii grow a bit: R eff,g = 28.3 ± 0.17 and R eff,r = 28.6 ± 0.19 .This is as expected, since there is more light in the outer parts of the galaxy (see Fig. B.2).In this work, using fixed ellipses, for SDSS we get R eff,g = 21.5 ± 1 and R eff,r = 22.1 ± 1 , while for DECaLS we get R eff,g = 23.5 ± 0.2 and R eff,r = 24.1 ± 0.2 .In these cases the effective radii are lower because the data are shallower and therefore the light beyond a given radius is not measured.The Gemini effective radii are also comparable to those published in previous papers.For example, using HST data, van Dokkum et al. (2018a) finds R eff,HST = 22.6 (F606W; similar to the r band).This lower value is reasonable considering that the outer part of the galaxy is brighter than originally thought by fitting a Sérsic law to the central region (see Fig. 5, middle panel).With INT deep data, Montes et al. (2021) measure 28.4 ± 0.3 using fixed ellipses, in very good agreement with our measurements when ellipticity and PA are left free.Finally, using Dragonfly, Keim et al. (2022) derives R eff,DF2 = 24.8 .This last measurement is obtained by a Sérsic fit with n = 0.6, so it may be missing the light in the outer parts of the galaxy.

Color and surface mass density profiles
To better understand the nature of NGC 1052-DF2, we examined its color g − r radial profile and its stellar surface mass density profile.This is shown in Fig. 6.The left and middle panels of Fig. 6 show the g − r color profile of NGC 1052-DF2 as a function of radial distance from the galaxy center for Gemini, INT, and DECaLS data.The color profile has been corrected of Galactic extinction using the following coefficients: A g = 0.08 and A r = 0.06 (Schlafly & Finkbeiner 2011).In the case of Gemini and DECaLS data, the ellipticity and PA are fixed as explained above.The INT profile corresponds to the one published in Montes et al. (2021) and was also obtained using fixed ellipses.The INT color profile agrees on a bluish behavior toward the edge of the galaxy.This is not so clear in the DECaLS data.The Gemini profile seems to lie between these two shallower datasets.The small offsets (<0.02 mag) between the different datasets correspond to possible uncertainties in the calibration of the different bands.We must remember that neither the filters nor their response are exactly the same across the different facilities, so small changes could be expected.The error bars of the color profile are obtained by summing in quadrature the errors on the surface brightness profiles in the g and r bands.
With the goal of interpreting the meaning of the radial color profile of NGC 1052-DF2, we explored what implications it might have for understanding the age and metallicity variations along the structure of the galaxy.We explored two scenarios.In the first (left panel), the age of the stellar population is fixed at 8.5 Gyr and we plot the expected color as a function of metallicity (dashed lines) using MILES models (Vazdekis et al. 2016).In the middle panel we repeat the procedure, but this time fixing the metallicity at [Fe/H] = −1.18 and plotting the expectations for different ages.The age and metallicity chosen in this exercise are motivated by the values found in the center of the galaxy using spectroscopy (see, e.g., Ruiz-Lara et al. 2019).The Gemini data are consistent with a moderate decrease in the A99, page 9 of 21 age and/or metallicity of the galaxy toward the outskirts.We performed a linear fit of the form (g − r) 0 (R) = m × R + c to the Gemini radial color profile.We find m = −(8.5 ± 1.2) × 10 −4 and c = 0.61 ± 0.03.The null hypothesis (i.e., m = 0) is rejected at 7σ.This analysis suggests that the observed color gradient in the galaxy could be: an age gradient, a variation in stellar metallicity, or a combination of both.However, as previously pointed out by Montes et al. (2021), it is difficult to distinguish between age and metallicity variations based on a single color.Nevertheless, the lack of HI detection in this galaxy (Chowdhury 2019) suggests that the galaxy may have a uniform age throughout its structure.In addition, previous observations of MUSE data of NGC 1052-DF2 by Fensch et al. (2019) indicate a possible metallicity gradient.Assuming a uniform age for the entire galaxy, the observed color gradient implies a possible increase in metal-poor stellar populations toward the outer regions.
Finally, in the right panel of Fig. 6 we show the stellar surface mass density radial profile of NGC 1052-DF2.To derive this profile we used the prescriptions given in Bakos et al. (2008).We used the g-band profile as the reference band, and we estimated the (M/L) g at each radial distance using the g − r profile (extinction corrected) using Roediger & Courteau (2015).We assumed a Chabrier (2003) initial mass function (IMF).Changing the IMF will scale the stellar mass density profile up and down, affecting the total stellar mass, but will not change the global shape of the profile.In addition to the Gemini data, we also show the stellar surface mass density using the DECaLS data (fixed ellipses, same as Gemini) and the INT data (fixed ellipses with PA and ellipticity found in Montes et al. 2021).The slight vertical offsets between the different profiles correspond to the small offsets in the g − r color that produce different (M/L) g .In general, there is a nice agreement up to 70 .Beyond this radius, the Gemini profile shows a strong indication of a truncation (edge) in the stellar mass profile.This is not so obvious in the case of the INT data.The reason for this is the shape of the surface brightness profile in the g band in the case of the INT data.As shown in Montes et al. (2021), the outer part of NGC 1052-DF2 is brighter in this band than what we find in the Gemini data.This is not observed in the r band where both works agree very well on the surface brightness profiles.The excess light in the INT g-band data is also responsible for the strong bluing in the g − r profile for that data.Interestingly, if the clear truncation on the Gemini profile is correct, the stellar surface mass density where it happens (i.e., ∼0.1 M pc −2 ) is within the expectation for the values found for the location of the truncation by Chamba et al.

NGC 1052-DF4
Figure 7 illustrates the gradual emergence of tidal features around NGC 1052-DF4 as the depth of the data increases.The four stamps shown, each covering a square region of 4 × 4 around the galaxy NGC 1052-DF4, correspond to observations made at different surface brightness limiting depths using a variety of facilities.While the presence of tidal streams remains invisible in the case of the SDSS and DECaLS, both the IAC80 and Gemini data show a clear excess of light in the outskirts of the object that we identify as tidal tails.

Radial surface brightness profiles
To obtain the surface brightness profile of NGC 1052-DF4, we followed the same strategy as for NGC 1052-DF2.First, we carefully masked the sources of contamination that we find in the field.We used exactly the same tools and codes explained before.The mask applied to the data is shown in the right panel of In the case of NGC 1052-DF4, it is particularly important to be careful when determining the background level, because the tidal tails of the galaxy cover a significant fraction of the pixels of the image.For this reason, only circular regions (with radius equal to 10 ) in the lower part of the image are used for this purpose.Next, we ran ellipse to derive the ellipticity and PA radial profiles.Our starting radius for fitting the shape of the isophotes is 30 .In Fig. 8 we compare these ellipticity and PA profiles with those obtained from other images (DECaLS in this work) and from the literature: IAC80 (Montes et al. 2020) and Dragonfly (Keim et al. 2022).The latter are the published ones and therefore correspond to those obtained from the sum of the g + r + i and g + r images, respectively.In the case of Gemini and DECaLS, we used only the r band.There is a reasonable agreement between the different datasets on the radial profile of For the PA profile, the agreement between the different telescopes is less clear.All data agree on a slight increase in the PA up to R ∼ 40 .Beyond this radius, while all the works converge on an increase in the PA, the significance of this change depends on the dataset.We think that the differences between the different images could be related to the removal of the contamination caused by the nearest star and NGC 1035.At these very low surface brightness levels, the size of the residuals after subtracting the light from these objects is crucial.Fortunately, as we can see in Fig. B.4, the effect on the shape of the surface brightness profile is not large and does not affect the interpretation of the tidal tails of the galaxy.
Following the same strategy as for NGC 1052-DF2, we measured a representative value of axis ratio and PA for the outer isophotes of NGC 1052-DF4 (40 < R < 80 ; see Fig. 8).The values we find are: b/a = 0.83 ± 0.05 and PA = 36 ± 4 deg.With these values fixed, we extracted the surface brightness radial profile of the galaxy.This is shown in Fig. 9. Beyond R = 40 , and in parallel with the increase in the ellipticity of the isophotes that we found earlier associated with the tidal tails of the galaxy, the surface brightness profile of the galaxy changes its slope, showing an excess of flux (an upward bend) with respect to the exponential decrease in the innermost region of the galaxy.This transition occurs at a surface brightness of about 27.5 mag arcsec −2 (r band).To compare the profile with previously published data where the ellipticity and PA of the isophotes are left free, we show in Fig. B.4 the surface brightness profile of NGC 1052-DF4 using Gemini data and isophotes free.The results are qualitatively the same.
Very deep imaging of this system also allows us to study the effect of the appearance of the tidal tails on the measurement of the effective radius of the galaxy.Using the very deep Gemini data, the effective radius derived by integrating the light profile of the galaxy in the r band is: R eff,r = 24.5 ± 0.1 (using fixed ellipses) and R eff,r = 25.8 ± 0.1 (free ellipses).This slight increase in the effective radius is as expected, since the free ellipses better follow the light distribution in the outer parts of the object.In shallower data such as SDSS, the effective radius is significantly smaller (R eff,r = 15.3 ± 1.3 ) because the data are not sensitive to the presence of tidal tails.For DECaLS we get R eff,r = 21.5 ± 0.5 .With respect to the values found in the literature, using Dragonfly, Keim et al. (2022) derived R eff,DF4 = 19.8 by applying a Sérsic fit with n = 0.85, while using IAC80 deep data, Montes et al. (2020) obtained 18 by applying a Sérsic fit with n = 0.86.In both cases, while the fit does a good job of describing the inner parts of the galaxy (R < 40 ), it falls short to characterizing the outer part of the object affected by the tidal tails.The different symbols correspond to data obtained with different telescopes.These profiles were all obtained using a fixed ellipticity and PA (based on the Gemini data) representative of the outer part of the galaxy (see the main text for details).

Discussion
As we mentioned in the Introduction, ultra-deep imaging of both NGC 1052-DF2 and NGC 1052-DF4 is expected to greatly constrain their formation scenario and shed light on their low dark matter content.In the following, we discuss what the data obtained in this work, in addition to the information obtained in previous works, tell us about the nature of these objects.

NGC 1052-DF2: No signatures of tidal disruption
All the data we analyzed and collected from the literature agree very well with the ellipticity and PA radial profiles of NGC 1052-DF2.The ellipticity of the object varies only slightly throughout its structure, with values fluctuating between 0.1 and 0.2.The PA does vary along the galaxy's structure, but beyond R = 40 it remains fairly constant around 50 • .These results on ellipticity and PA variation were already known and discussed in the literature.Despite the remarkable agreement between the different datasets, the interpretation is very different depending on the research group.For Montes et al. (2021), the radial profile of the PA can be interpreted as a transition from a pressure-supported region (R < 40 ) to a disk-like outer part.For Keim et al. (2022), however, the change in the PA plus the change in the surface brightness profile (with transitions from a low Sérsic index n < 0.6 to an exponential shape) is associated with a tidal disturbance.
The exponential decrease in the radial surface brightness profile beyond R > 40 is not the only argument that Montes et al. (2021) use to support the idea of a rotationally dominated structure in the outer part of NGC 1052-DF2.In fact, these authors argue that the dynamics of the GC system of this galaxy is also consistent with a rotating disk (Lewis et al. 2020).
The new Gemini ultra-deep imaging is useful to further reinforce or reject the two contending scenarios (i.e., rotating disk vs. tidal disturbance).In this work we find no signatures that support NGC 1052-DF2 being tidally disrupted.It should be emphasized that the change in PA occurs at a relatively bright surface brightness (µ g ∼ 27 mag arcsec −2 ).If this feature was related to the development of tidal tails in the galaxy, the Gemini data should show the S-shaped structure that we can see in NGC 1052-DF4, for example.This is not the case.While it is arguable that the surface brightness profile of Keim et al. (2022) shows an excess of brightness in the outer part compared to the rest of the data (see, e.g., Fig Another way to look at this is to make a direct comparison between the surface brightness profile of NGC 1052-DF2 (with no signatures of tidal distortions) and NGC 1052-DF4 (with obvious distortions).This is done in Fig. 10.The two profiles are very different.NGC 1052-DF2 shows an exponentially decreasing surface brightness profile up to R = 60 , with a hint of a break in the profile at R = 60 .In contrast, NGC 1052-DF4 shows an excess of light, an upward bend, in the outer regions (R > 40 ) of the surface brightness profile.
The absence of tidal features around NGC 1052-DF2 up to the surface brightness limits of the image (µ g ∼ 31 mag arcsec −2 ; 3σ; 10 × 10 boxes) is useful to rule out a number of formation scenarios for this galaxy.For example, either an ongoing (Katayama & Nagamine 2023) or a close flyby (Moreno et al. 2022) of NGC 1052-DF2 around NGC 1052 should have produced visible signatures of the interaction in the outskirts of this galaxy in images as deep as the ones we have here.Therefore, the likelihood that either of these scenarios has occurred should be low.The absence of tidal distortions is also useful to place strong constraints on the possible association of NGC 1052-DF2 with its closer (in projection) massive galaxies: NGC 1052 and NGC 1042.This has been done in detail by Montes et al. (2021).According to these authors, under the hypothesis that NGC 1052-DF2 is physically connected to NGC 1052, it should suffer from a tidal field that should have produced distortions well inside the galaxy body (55 ± 15 ).This is rejected by the observations.However, if the galaxy is associated with NGC 1042, no signatures of tidal distortions are expected down to ∼150 .This last scenario cannot be ruled out even with this deeper dataset.
Another alternative way of testing the tidal disruption scenario for NGC 1052-DF2 is by looking at the spatial distribution of its GC system.The object has eleven spectroscopically confirmed GCs (Shen et al. 2021) with an additional four proposed in Montes et al. (2021).In Fig. 11 (left panel) we show the distribution of the GCs orbiting NGC 1052-DF2.It is remarkable that except for two, the vast majority of the GCs are directly on top the body of the object.This is suggestive that the object has not undergone a tidal disruption process.One way of quantifying this statement is by comparing the effective radius of the galaxy with the radius enclosing half of the GCs of the system (R GC ).
We determined R GC by computing the median of the distances of the GCs in logarithmic units.If we consider only the spectroscopically confirmed GCs (Shen et al. 2021), we obtain R GC = 25.9 .Adding the GC candidates suggested in Montes et al. (2021) changes this parameter to R GC = 20.6 .The effective radius of NGC 1052-DF2 varies from 22.1 (SDSS; shallower data) to 25.6 (Gemini; deeper data).Consequently, considering all these possible variations, we find that 0.8 < R GC /R e < 1.2.These low values for the ratio of these two quantities are also found in other ultra diffuse galaxies (UDGs) that show no evidence of tidal distortions (Saifollahi et al. 2021(Saifollahi et al. , 2022)).In short, both very deep imaging from Gemini plus the GCs distribution do not favor a scenario where the low dark matter content of NGC 1052-DF2 can be associated with a removal of the dark matter by tidal forces.

Which galaxy is disrupting NGC 1052-DF4?
There is a consensus in the community that the excess of light in the outer parts of NGC 1052-DF4 found by Montes et al. (2020) can be interpreted as this galaxy being tidally perturbed by a nearby massive neighbor.In support of this scenario, the analysis of the GC distribution of this object also shows that it is significantly extended compared to other galaxies that show no sign of being tidally perturbed (see the right panel of Fig. 11).In fact, the ratio of the radius containing half of the confirmed (van Dokkum et al. 2019) and candidate (Montes et al. 2020) GCs to the effective radius of this galaxy (R eff,r from Gemini data using fixed ellipses in this work) is R GC /R e ∼ 1.7.If we assume only spectroscopically confirmed GCs we obtain the same R GC = 41.7 , which would lead to the same ratio of 1.7.
The above ratio is a lower limit.We can use the effective radius measured in shallower data such as the SDSS (see Sect. 3.2) as a more representative measure of the bulk of the galaxy that is not affected by the light from the object in the tidal tails.When we do this, the ratio between the two radii increases to R GC /R eff ∼ 2.7.This is around a factor of 2.5 larger than what we observe in UDGs that do not show tidal distortions (Saifollahi et al. 2022).
While there is no debate about the excess of light in the outer parts of NGC 1052-DF4, there is still no consensus about which galaxy is producing the tidal field responsible for perturbing the A99, page 13 of 21  2022) suggest that the disturbance is generated by the massive elliptical NGC 1052.There is no way to answer this question definitively without a precise measurement of the distance to all galaxies in the NGC 1052 group (not just NGC 1052-DF2 and NGC 1052-DF4).A comprehensive analysis of the distances to the galaxies in the line of sight of NGC 1052 by Monelli & Trujillo (2019) suggests that there are two groups in the projection.The study suggests that both NGC 1052-DF2 and NGC 1052-DF4 belong to the foreground group, at 13.5 Mpc, with the more massive galaxies NGC 1042 and NGC 1035.This would imply that the disruptor is NGC 1035.Further evidence for this can be found in Román et al. (2021).This result is not helped by the latest published distance to NGC 1052-DF4, which is D = 20.0 ± 1.6 Mpc (Danieli et al. 2020), used by Keim et al. (2022) to assign it to NGC 1052.However, to make the scenario even more complex, this distance could be in tension with the recent accurate distance estimate to the most massive galaxy in the group NGC 1052 (D = 17.9 +0.3 −0.6 Mpc; Jacoby et al. 2023), which implies a relative distance of ∼2.1 Mpc, larger than the virial radius of NGC 1052 (390 kpc; Forbes et al. 2019).Consequently, the relationship between NGC 1052 and NGC 1052-DF4, as well as between NGC 1035 and NGC 1052-DF4, is still not clear.The distance debate is therefore far from settled and will require further work in the future.To help resolve this debate, deep data such as those presented in this paper, but with a larger area coverage that includes NGC 1035 in the same pointing, would also be useful to test whether NGC 1035 also shows signs of perturbations that could be related to a direct gravitational interaction with NGC 1052-DF4.

Conclusion
The low velocity dispersion of both GCs and stars in NGC 1052-DF2 and NGC 1052-DF4 (Shen et al. 2023) may indicate that the dark matter content of these galaxies is low or zero.One possible explanation for this is that the dark matter in these galaxies has been removed via interactions with other objects (Ogiya 2018;Jackson et al. 2020;Macciò et al. 2021;Ogiya et al. 2022).These mechanisms should leave an excess of brightness in the form of tidal tails in the outer part of the galaxies (e.g., Johnston et al. 2002;Moreno et al. 2022).With sufficiently deep imaging (µ V ∼ 30.5 mag arcsec −2 ), these tidal tails could be detected (see, e.g., Moreno et al. 2022;Katayama & Nagamine 2023).In the current work, we explore this scenario using new ultra-deep images of these galaxies taken with the Gemini telescopes, 1 mag deeper than those previously published.
The new ultra-deep images of NGC 1052-DF4 confirm the presence of tidal tails already observed in previous works (Montes et al. 2020;Keim et al. 2022).However, the debate continues as to which galaxy is responsible for the tidal destruction of this object.To solve this puzzle it would be necessary to know the exact distance to all the galaxies in the FoV of NGC 1052.In the case of NGC 1052-DF2, the new Gemini images show no evidence that the galaxy is being destroyed, to a surface brightness limit of µ g = 30.9mag arcsec −2 (3σ; 10 × 10 boxes).This observation, together with the fact that the distribution of GCs in this galaxy is compact (R GC /R eff ∼ 1), suggests that the possible absence of dark matter in this galaxy is not due to its removal.Therefore, either there are no known mechanisms that can sustain the dynamics of this galaxy without dark matter, or the absence of dark matter could be due to an incorrect estimation of the distance to this galaxy (van Dokkum et al. 2018a,b;Trujillo et al. 2019;Monelli & Trujillo 2019;Zonoozi et al. 2021) or an incorrect hypothesis about the dynamical equilibrium of this galaxy (Lewis et al. 2020;Montes et al. 2021).Further analysis in this regard may help resolve this conflict.and Employment, through the Regional Budget of the Autonomous Community.G.G. acknowledges support from IAC project P/302304 and through the PID2022-140869NB-I00 grant from the Spanish Ministry of Science and Innovation which is partially supported through the state budget and the regional budget of the Consejería de Economía, Industria, Comercio y Conocimiento of the Canary Islands Autonomous Community.M.M. and I.T. acknowledge support from IAC project P/302302.Facilities: Gemini: GMOS-S (Program ID: GS-2021B-FT-104), GMOS-N (Program ID: GN-2021B-FT-108).Software: SExtractor (Bertin & Arnouts 1996), SCAMP (Bertin 2006), SWarp (Bertin 2010), Gnuastro (Akhlaghi & Ichikawa 2015), Photutils (Bradley et al. 2019), Matplotlib (Hunter 2007)   We developed a special observing strategy to deal with vignetting caused by the arm of the guide star and to minimize scattered light effects on GMOS cameras.This strategy involves implementing a dithering pattern with an offset of 10 , which is applied both over the galaxy and over an adjacent region of the sky free of bright objects.For every six images of the galaxy, we repositioned the telescope to take five images of a nearby region of the sky.We then returned the telescope to the original location of the galaxy.This approach allowed us to build the master flat field, a critical component of deep imaging, by combining only these sky images.By adopting this observing strategy, we aim to mitigate vignetting problems, reduce scattered light contamination, and optimize the quality of our imaging data obtained with the Gemini telescopes.
In the upper panel of Variations in airmass correspond to differences in image quality, which we can characterize by the standard deviation of the background pixels (in units after photometric calibration, i.e., nanomaggies in our case).The left plot in  (Keim et al. 2022), and Gemini (this work).In all cases the ellipticity and PA are allowed to vary freely.We also include the Gemini profile (green dots) resulting from fixing the ellipticity and PA (see the main text for details).The IAC80 data correspond to images where the g-, r-, and i-band filters have been combined, and for Dragonfly the g and r bands have been combined, while in the case of DECaLS and Gemini we show only the r-band results.The IAC80 and Dragonfly profiles have been shifted vertically by 0.1 mag and 0.25 mag, respectively, to match the r-band profiles only, to facilitate comparison between different telescopes.

Fig. 1 .
Fig. 1.Three stages of subtraction of the scattered light emitted by the star HD 16873 near NGC 1052-DF2.The panels depict the complete mosaic of the field (5 × 5 ) using g-filter Gemini data.The left panel illustrates the original stacked image, while the middle panel shows the model fitted to the starlight.Subtracting the star model from the original image yields the residual image shown in the right panel.

Fig. 2 .
Fig. 2. Three stages of subtraction of the scattered light emitted by NGC 1035 and the star (Gmag = 14 mag) near NGC 1052-DF4 (similar to Fig. 1).The panels depict the mosaic of the field (5 × 5 ) using a color image generated with r-filter Gemini data.The left panel shows the original stacked image, and the middle panel illustrates the models of the star and the galaxy NGC 1035.The right panel shows the residual image generated by subtracting the models from the original image.
The mask is created to avoid contamination from foreground and background sources.A first step in generating the mask is to run Noisechisel and Segment.As a second step, to improve the detection of the contamination on top of our galaxies, we applied the technique of unsharp masking to remove the light of the galaxy and thus easily detect the individual contamination sources on top of it.Unsharp masking consists of smoothing the original image with a Gaussian filter of the width of the point-like sources and subtracting the smoothed image from the original.Then we applied Noisechisel and Segment again to mask the remaining undetected objects.The mask applied to NGC 1052-DF2 field is shown in the left panel of Fig. C.1.
Fig.3.Appearance of NGC 1052-DF2 at different surface brightness depths.The FoV of each stamp is 3 × 3 .The limiting surface brightness of the r-band images (SDSS, DECaLS, INT, and Gemini) is indicated on the panels using a common metric of 3σ in 10 × 10 boxes.The color of the images is a combination of the g and r filters.The black and white background is from the g-band filter.As the image gets deeper, the apparent size of the galaxy also increases, with no obvious signatures of distortions in its periphery.

Fig. 4 .
Fig. 4. Ellipticity (top) and PA (bottom) radial profiles of NGC 1052-DF2 using different datasets.The profiles for the DECaLS and Gemini images were derived by our group using ellipse, starting the fit at 40 .The radial profiles for INT and Dragonfly are those published in Montes et al. (2021) and Keim et al. (2022), respectively.The green horizontal lines correspond to the average ellipticity (top) and PA (bottom) values in the radial range between 40 and 80 .

Fig. 5 .
Fig. 5. Surface brightness profiles of NGC 1052-DF2 in the g band (upper panel), r band (middle panel), and the combination of the two filters (lower panel).The different symbols correspond to data obtained with different telescopes.These profiles were obtained using a fixed ellipticity and PA representative of the outer parts of the galaxy (see the main text for details).The dashed line in the middle panel corresponds to the Sérsic fit to the F606W band of the HST ACS camera (n = 0.6 and R e = 22.6 ; van Dokkum et al. 2018a).
Fig. B.2. Interestingly, the feature is less obvious when the ellipticity and PA are allowed to vary when deriving the surface brightness profiles (see Fig. B.2).The significance of this is discussed in Appendix B. We also include in the middle panel of Fig. 5 the Sérsic fit to the Hubble Space Telescope (HST) F606W band of the HST ACS camera made to the galaxy in van Dokkum et al. (2018a).This fit is a good representation of the light distribution of the galaxy in the innermost region, but differs from what has been observed with deeper data beyond 40 .

Fig. 6 .
Fig. 6.Color (g − r) 0 and stellar surface mass density profiles of NGC 1052-DF2 using data from Gemini, INT (Montes et al. 2021), and DECaLS.In the left and middle panels, the horizontal dashed lines represent the expected color values based on the Vazdekis et al. (2016) models for different metallicities at a fixed age of 8.5 Gyr (left panel) and for different ages at a fixed [Fe/H] of −1.18 (middle panel).The right panel shows the stellar mass density profile of NGC 1052-DF2.

Fig. 7 .
Fig. 7. NGC 1052-DF4 appearance as the image depth increases.The four stamps presented, each covering a 4 × 4 region around the galaxy NGC 1052-DF4, correspond to observations obtained at different surface brightness limiting depths in the r band (3σ, 10 × 10 ) using various facilities.The tidal streams become visible when reaching surface brightness levels fainter than µ r ∼ 29 mag arcsec −2 .In all stamps, to represent different pixel intensities, the white corresponds to faint pixels and black represents brighter features.The black and white background is generated used the r band only.The color part of the stamps is created using rescaled HiPERCAM data from Montes et al. (2020) using the g, r, and z filters.

Fig. 8 .
Fig. 8. Ellipticity (top)  and PA (bottom) radial profiles of NGC 1052-DF4 using different datasets.The profiles for the DECaLS and Gemini images were derived by our group using ellipse, starting the fit at 30 .The radial profiles for IAC80 and Dragonfly are those published inMontes et al. (2020) andKeim et al. (2022), respectively.The green horizontal lines correspond to a fit to the average ellipticity (top) and PA (bottom) values in the radial range between 40 and 80 .

Fig. 9 .
Fig.9.Surface brightness profiles of NGC 1052-DF4 in the Sloan r band.The different symbols correspond to data obtained with different telescopes.These profiles were all obtained using a fixed ellipticity and PA (based on the Gemini data) representative of the outer part of the galaxy (see the main text for details).

Fig. 10 .
Fig. 10.Surface brightness profiles in the r band of NGC 1052-DF4 (light blue) and NGC 1052-DF2 (orange) obtained from fixed ellipses with Gemini ultra-deep imaging.The shaded regions correspond to the error bars associated with the measurements.The two profiles are strikingly different.NGC 1052-DF4, which is clearly tidally disturbed, shows an excess of light (R > 40 ) that is not seen in the case of NGC 1052-DF2.
. B.2), this excess can be understood as the result of a different sky background subtraction.The background subtraction used in the present work, in Montes et al. (2021), as well as in Keim et al. (2022, see their Appendix B) is done by removing a constant value.Therefore, it cannot destroy a possible S-like shape around NGC 1052-DF2.

Fig. 11 .
Fig. 11.GC distributions around NGC 1052-DF2 (left panel) and NGC 1052-DF4 (right panel).Red circles indicate the positions of the spectroscopically confirmed GCs (van Dokkum et al. 2019; Shen et al. 2021), while the positions of the GC candidates found in Montes et al. (2020, 2021) are denoted with orange circles.The radii of the green circles (R GC ) represent the spatial location where half of spectroscopically confirmed GCs are found.The radii of the light blue circles instead, represent the spatial location where half of all GCs candidates are found.

Fig
Fig. A.2. Similar to Fig. A.1, but for NGC1052-DF4.Upper panel: NGC1052-DF4 field over an area of 18 × 18 using data from the SDSS r band.The dimension of the GMOS camera and the centers of the dithering pattern around the galaxy and for the sky images are highlighted.Lower-left panel: Final stack of the galaxy in black and white background, using the GMOS r filter.The color part of the Gemini stamp is created using rescaled HiPERCAM data from Montes et al. (2020) using the g, r, and z filters.Lower-right panel: Sky region that was observed to obtain the flat field of the data.
Fig. A.1 we present the SDSS mosaic of the sky region surrounding NGC1052-DF2 (18 × 18 ).The background image is created using the SDSS r filter downloaded from the DR12 Science Archive Server (SAS).The CCDs of the Gemini camera (total FoV = 5 × 5 ) are outlined in black on top of the galaxy.Five representative central positions of the dithering pattern over the galaxy are marked with blue crosses.The central coordinates of each sky frame are highlighted in red.The lower-left panel shows the final Gemini stack of galaxy frames, while the lower-right panel shows the final distribution of sky regions used to construct the flat field of the images.In the same way, Fig. A.2 illustrates the observing strategy of NGC1052-DF4.
Fig. A.3.Airmass and sky variation during observations.Left panel: Relationship between image quality (characterized by the standard deviation of the background pixels in nanomaggies units) and the airmass of the observations.Poorer quality images (larger values of the background standard deviation) are observed in those images with larger airmass values.Right panel: Variation in the airmass during the observation run.

Figure C. 1 Fig. C. 1 .
Figure C.1 shows the masks used in both Gemini images to compute the surface brightness profiles.

Table 1 .
Characteristics of the Gemini observing run conducted to achieve ultra-deep imaging of NGC 1052-DF2 and NGC 1052-DF4.ObjectCamera Filter t exp,obj n frames,obj t exp,sky n frames,sky t total,obj µ limit (3σ; 10 × 10 boxes)

Table 2 .
Surface brightness limiting depths for the different datasets analyzed in this paper.
of Universities for Research in Astronomy (AURA) under a cooperative agreement with the National Science Foundation on behalf of the Gemini Observatory partnership: the National Science Foundation (USA), National Research Council (Canada), Agencia Nacional de Investigación y Desarrollo (Chile), Ministerio de Ciencia, Tecnología e Innovación (Argentina), Ministério da Ciência, Tecnologia, Inovações e Comunicações (Brazil), and Korea Astronomy and Space Science Institute (Republic of Korea).The Legacy Surveys consist of three individual and complementary projects: the Dark Energy Camera Legacy Survey (DECaLS; Proposal ID 2014B-0404; PIs: David Schlegel and Arjun Dey), the Beijing-Arizona Sky Survey (BASS; NOAO Prop.ID 2015A-0801; PIs: Zhou Xu and Xiaohui Fan), and the Mayall z-band Legacy Survey (MzLS; Prop.ID 2016A-0453; PI: Arjun Dey).DECaLS, BASS and MzLS together include data obtained, respectively, at the Blanco telescope, Cerro Tololo Inter-American Observatory, NSF's NOIRLab; the Bok telescope, Steward Observatory, University of Arizona; and the Mayall telescope, Kitt Peak National Observatory, NOIRLab.Pipeline processing and analyses of the data were supported by NOIRLab and the Lawrence Berkeley National Laboratory (LBNL).The Legacy Surveys project is honored to be permitted to conduct astronomical research on Iolkam Du'ag (Kitt Peak), a mountain with particular significance to the Tohono O'odham Nation.NOIRLab is operated by the Association of Universities for Research in Astronomy (AURA) under a cooperative agreement with the National Science Foundation.LBNL is managed by the Regents of the University of California under contract to the US Department of Energy.M.M. acknowledges support from the Project PCI2021-122072-2B, financed by MICIN/AEI/10.13039/501100011033,and the European Union "NextGenerationEU"/RTRP. I.T. 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 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, financed by the Ministry of Science and Innovation, through the State Budget and by the Canary Islands Department of Economy, Knowledge , NumPy (van der Walt et al. 2011), Astropy (Astropy Collaboration 2018).