© The Authors 2024

1. Introduction

In the currently accepted Λ cold dark matter (CDM) scenario, the Universe is almost flat and expanding, and the expansion is accelerating (Riess et al. 1998; Perlmutter et al. 1999) ≈30% (Ω_m ≈ 0.3) of the Universe is composed of baryons and dark matter (DM), and the cosmic acceleration is due to the remaining ≈70% (Ω_Λ ≈ 0.7), represented by so-called dark energy. Our understanding of dark energy is very poor. We believe that it exerts a negative pressure, and has an equation of state with w ≈ −1, where w is defined as the ratio between pressure and energy density, p/ρc². The ΛCDM model successfully describes the Universe at large scales (≳1 Mpc), but has difficulty explaining some properties related to the formation of structures at smaller scales, like the sub-halo structures in galaxy clusters (Grillo et al. 2015; Meneghetti et al. 2020) and the value of the inner slope of DM halos (e.g., Gnedin et al. 2004; Newman et al. 2013a,b; Martizzi et al. 2012). Moreover, in extended cosmological models, currently available data cannot accurately measure the values of Ω_m and Ω_Λ, and different models can be reconciled with the observations (e.g., Motta et al. 2021). This motivates the continuous investigation and testing of the ΛCDM model, through new projects and observations. In this context, the use of different and independent cosmological probes is crucial, as they are subject to different systematics and degeneracies. Thus, they can offer valuable help in investigating the current tensions in cosmology (Verde et al. 2019; Moresco et al. 2022).

Strong gravitational lensing is an extremely powerful tool for extragalactic and cosmological studies (e.g., Bartelmann 2010; Treu 2010). Among the many applications, such as characterizing the total and DM mass distributions of clusters of galaxies (Acebron et al. 2022; Bergamini et al. 2023a; Granata et al. 2023) and galaxies (Vegetti et al. 2012; Schuldt et al. 2019; Ballard et al. 2024) acting as lenses, gravitational lensing can be used to probe the geometry of the Universe. Beside being particularly effective at measuring the value of the current expansion rate of the Universe (the Hubble constant, H₀) through observations of multiply lensed variable sources, such as quasars or supernovae, both on galaxy and cluster scales (e.g., Refsdal 1964; Suyu et al. 2017; Grillo et al. 2018, 2020, 2024; Birrer et al. 2019; Wong et al. 2020; Rusu et al. 2020; Shajib et al. 2023), it also allows one to measure the values of Ω_m, w, and Ω_k (the latter parametrizing the curvature of the Universe, Ω_k = 0 in a flat geometry), when kinematic data for lens galaxies are available (e.g., Grillo et al. 2008; Cao et al. 2012) or in systems where two or more sources are multiply imaged by the same deflector (Tu et al. 2009; Collett & Auger 2014; Tanaka et al. 2016; Smith & Collett 2021, and with clusters of galaxies).

Massive clusters of galaxies can produce tens to hundreds of multiple images from several background sources and in the last decade they have been employed to measure the values of Ω_m, Ω_Λ, w, and Ω_k (e.g., Jullo et al. 2010; Caminha et al. 2016, 2022; Grillo et al. 2024), thanks also to the advent of very deep integral field spectroscopic (IFS) observations, which represent the most effective way to spectroscopically confirm and discover lensed sources. For instance, the number of spectroscopically confirmed multiple images lensed by the Hubble Frontier Field galaxy cluster MACS J0416.1–2403 has increased, in less than 10 years, from 10 sources observed in 30 multiple images (Grillo et al. 2015) to 88 sources lensed into 237 multiple images (Bergamini et al. 2023b).

Unfortunately, the total mass distribution of clusters of galaxies is usually very complex, and to be properly modeled it requires several mass components, whose parameters may be degenerate with the cosmological ones. On the other hand, galaxy-scale systems are in general easier to model and the lens can often be described with an effective single total mass profile, but the limited number of background sources makes them prone to be affected by the mass-sheet degeneracy (Schneider 2014). Galaxy or group-scale systems with a larger number (≥3) of lensed background sources might represent the best compromise between these two regimes to learn about the geometry of the Universe. They allow one to measure the values of Ω_m and w, independently of that of H₀, but only a few of these systems are known to date and are suitable for cosmological studies (e.g., Smith & Collett 2021). Their analysis, even individually, can offer competitive estimates of the cosmological parameter values, and will pave the way for the exploitation of a larger number of systems of this kind, which are foreseen to be discovered with Euclid and the Vera C. Rubin Observatory – Legacy Survey of Space and Time (LSST) (Collett & Auger 2014; Li et al. 2024).

In this paper, we extend the study published in Bolamperti et al. (2023) (hereafter, B23) on SDSS J010049.18+181827.7 (hereafter, SDSS J0100+1818), a strong lensing system (Fig. 1) included in the Cambridge And Sloan Survey Of Wide ARcs in the skY (CASSOWARY) survey (Belokurov et al. 2009; Stark et al. 2013) as a candidate fossil system at z = 0.581 (Johnson et al. 2018a). In our previous work, we developed a strong lensing model of the system from the observed positions of the four multiple images of sources A and B (spectroscopically confirmed) and of the two multiple images of source C (with its redshift as a free parameter), and from the extended surface brightness distributions of the multiple images from the Hubble Space Telescope (HST) data. We employed best-fit models to measure the cumulative total mass profile of the deflector, disentangle the DM and baryonic mass distributions, and reconstruct the background sources. Now, we leverage new data taken with the Multi Unit Spectroscopic Explorer (MUSE) on the Very Large Telescope (VLT). These IFS observations allow us to measure the redshift of source C, discover two additional strongly lensed objects (E and F), spectroscopically confirm the group members in the MUSE field of view, and measure the stellar velocity dispersion profile of the brightest galaxy. We develop enhanced strong lensing models by including this information, which also allow us to measure the posterior probability distributions of the cosmological parameters Ω_m and w.

Fig. 1. HST F160W image of the SDSS J0100+1818 strong lensing system studied in this work. Letters label multiple images of the same background source. E and F do not show a continuum in the HST image and the orange and brown crosses mark the positions of their Lyα emission line peaks, detected by MUSE. The dashed gray circle represents an aperture of 50 kpc, the approximate average Einstein radius of the system.

This paper is organized as follows. In Sect. 2, we summarize the currently available data for SDSS J0100+1818, focusing on the MUSE data reduction and spectra extraction. In Sect. 3, we describe the SDSS J0100+1818 system, characterizing the multiply imaged sources A, B, C, E, and F (excluding D, which is not a secure system), the deflector with its velocity dispersion profile, and the group members. In Sect. 4, we describe the enhanced strong lensing models developed and the relative results. In Sect. 5, we show the strong lensing models with variable values of the cosmological parameters Ω_m and w. We discuss the results in Sect. 6 and draw conclusions in Sect. 7. Unless otherwise specified, as in Sect. 5, throughout this work we assume H₀ = 70 km s⁻¹ Mpc⁻¹, Ω_m = 0.3, and Ω_Λ = 0.7. In this model, 1 arcsec corresponds to a linear size of 6.585 kpc at the deflector redshift of z = 0.581. All magnitudes are given in the AB system (Oke 1974) and are measured in the HST F160W filter, unless otherwise specified.

2. Observations and data reduction

2.1. Ancillary data

SDSS J0100+1818 benefits from a large sample of both photometric and spectroscopic observations, detailed, with the adopted data reduction procedures, in B23. In particular, we observed SDSS J0100+1818 with the HST Wide Field Camera 3 (WFC3; program GO-15253; PI: R. Cañameras), spending one orbit in each of the two F438W and F160W filters, resulting in PSF FWHMs of 0.086″ and 0.187″ in F438W and F160W, respectively. We complemented photometric observations with VLT/X-Shooter (Vernet et al. 2011) spectroscopy (program 091.A-0852, PI: L. Christensen), obtaining the spectrum of the main deflector and of three out of four multiple images of families A and B, spectroscopically confirming them, and measuring the lens stellar velocity dispersion.

2.2. Integral field spectroscopy observations with VLT/MUSE

We added IFS data of SDSS J0100+1818 by observing it for 4.8 hours with VLT/MUSE, from October to December 2022 (program 0110.245R, PI: A. Bolamperti). The observations were taken with seeing ≲1″, clear sky conditions, an airmass ≲1.4, and with the wide-field mode, resulting in a 1′ × 1′ field of view and a spatial sampling of 0.2″ pix⁻¹. We divided the observations into six different observation blocks composed of two exposures of 1440 s each. Each exposure was taken with four different rotation angles (90° apart) and applying small offsets to the center, to optimize the final data quality.

We reduced the data with the standard MUSE pipeline v2.8.9 (Weilbacher et al. 2020), following the procedure detailed in Caminha et al. (2019), through the ESO Recipe Execution Tool (EsoRex; ESO CPL Development Team 2015) pipeline. In summary, we corrected all the raw exposures using the associated BIAS, FLAT, and ILLUMINATION calibrations. We then applied wavelength and flux calibrations and created a PIXTABLE relative to each exposure. At the end, we combined all of them into a final stacked datacube, and defined the astrometry with respect to the HST F160W image. To enhance the background sky subtraction, we made use of the Zurich Atmosphere Purge (ZAP; Soto et al. 2016) tool.

The resulting datacube spans in wavelength from 4750 Å to 9350 Å, with a constant 1.25 Å pix⁻¹ sampling. The final reduced MUSE data cube has a median PSF FWHM of 0.8″, covers a square of approximately 400 kpc on a side at the redshift of the main deflector, and is centered on it. This pointing includes all the multiple images previously identified, as well as the most likely group members from our photometric analysis.

We identified and measured the source redshifts following the multistep procedure described in Caminha et al. (2019). Firstly, we ran SExtractor (v2.28.0, Bertin & Arnouts 1996) on a cutout of the HST F160W image that covers the final MUSE datacube field of view (shown in Fig. 2), to detect the position of all the sources. We then extracted a spectrum from the datacube within a circular aperture with a radius of 0.8″ centered on each detected position. This aperture was chosen to represent a good compromise between collecting as much signal as possible and not including nearby contaminants in the most crowded regions. For some particular systems of distorted arcs hosting multiple images, or with faint and extended objects, we assumed specifically designed regions to extract the spectra, which allowed us to obtain a better signal-to-noise ratio (S/N). We complemented this catalog by adding sources that present clear emission lines in the MUSE datacube, but are not detected in the HST F160W continuum, through visual inspection and the use of the Cube Analysis and Rendering Tool for Astronomy (CARTA; Comrie et al. 2021) tool. We extracted their spectra within circular apertures with a radius of 0.8″ centered on the luminosity peak of the detected emission lines.

Fig. 2. Group members (squares) identified as those galaxies whose redshift is within 0.581 ± 0.011, equivalent to the mean redshift of the group with a spread of 2000 km s−1 (rest-frame). The colored regions represent the apertures within which we extracted the spectra of the multiple images of A and B (blue), C (green), D (purple), E (orange), and F (brown). The spectra of E and F, which do not present an HST-detected continuum, are extracted within circular apertures with a radius of 0.8″. The dashed gray curve shows a circle with a radius of 50 kpc, which represents approximately the average Einstein radius of the system.

We measured the redshift values of the objects in the catalog with a dual automatic and visual procedure, identifying clear spectral features like emission and/or absorption lines, and continuum breaks. We made use of the software Marz (Hinton et al. 2016), which performs an automatic fitting of each spectrum with different templates. We found the automatic procedure reliable in the high S/N regime (e.g., to identify the bright group members), while we visually inspected the faint objects or those with only emission lines detected. We assigned for each redshift measurement a quality flag (QF), defined to be QF = 1 for unreliable measurements, QF = 2 for possible measurements based on faint spectral features, QF = 3 for secure measurements from more than one absorption and/or emission lines, and QF = 9 for measurements based on a single narrow line emission. Most of the objects with QF = 9 present features that allowed us to identify their nature, like the Lyα line or a doublet, and their reliability is similar to that of the QF = 3 objects (see Grillo et al. 2015; Balestra et al. 2016). The catalog used in this work contains 133 redshift measurements with QF ≥ 2 (65 with QF ≥ 3), divided into 3 (3) stars, 22 (19) group members, 27 (10) foreground galaxies, 80 (32) background galaxies, and 1 (1) background quasar.

3. The SDSS J0100+1818 system

SDSS J0100+1818, (RA, Dec) = (01:00:49.18, +18:18:27.79), was introduced in a later release of the CASSOWARY catalog. In our previous work B23, we determined a robust lens redshift, z = 0.581, and a stellar mass value of (1.5 ± 0.3) × 10¹² M_⊙ for the most luminous galaxy. We assumed different mass distributions for the main lens and, through a strong lensing analysis, we measured the total mass profile of the deflector. We have recently spotted a flaw in our code that caused the deflector total mass profile to be slightly underestimated. In the following, and in particular in Sect. 6.1, we refer to the corrected values. They consist of a total mass value of (1.16 ± 0.01) × 10¹³ M_⊙ within 42 kpc and a stellar over total mass fraction of (38 ± 9)% at the half-light radius (R_e = 9.3 kpc) of the main lens galaxy. These values are consistent with the previous estimates, and do not affect the results and discussion presented in B23, which remain valid.

3.1. Multiple images of A and B

We confirm the joint redshift measurement of z = 1.880 for the two source components forming image families A and B (Fig. 1) that we found in B23. Previously, based on X-Shooter data, we performed a joint analysis of the targeted multiple images, A1/B1, A3/B3, and A4/B4. We made use of lines detected at about 10740 Å in the binned 2D spectra, whose width is consistent with the [OII]λ3727, 3729 doublet, of a faint detection of [OIII]λ5007 in A3/B3, and of the lack of additional line detections over the spectral range covered by X-Shooter. In the MUSE cube, we extracted the spectra over four regions, shown in Fig. 2, which include the compact and extended emission from all four images. We detected the Mg IIλ2800 doublet at about 8065 Å, five absorption features of FeII (Fe IIλ2344, 2374, 2382, 2586, 2600), three of Al (Al IIλ1671, AlIIIλ1854, Al IIIλ1862), and a tentative detection of He Iλ3188 and [CII]λ2326. These features are shown, for each of the A1/B1, A2/B2, A3/B3, A4/B4 images and stacked, in Fig. 3.

Fig. 3. MUSE 1D spectra of the multiple images A1+B1, A2+B2, A3+B3, A4+B4, and stacked (top row), extracted within the blue apertures highlighted in Fig. 2. The spectra were smoothed with a ∼ 2.5 Å boxcar filter. We zoom in around five wavelength ranges where we observe the main absorption and emission lines used to confirm the redshift measurement of 1.880. The blue lines represent the observed spectra in units of 10−18 erg s−1 cm−2 Å−1, and the gray regions indicate the 1σ uncertainties.

3.2. Multiple images of C

By combining observations in the HST F160W and F438W filters with the first strong lensing models of the system, we previously identified another candidate lensed background source, with two multiple images labeled as C1 and C2 in Fig. 1. Since C1 and C2 have similar colors and their positions were correctly predicted by our initial strong lensing models, it was highly likely that a single source was observed multiple times. Although we lacked spectroscopic confirmation, as they were not included in the X-Shooter pointing, we considered the redshift value of source C as a free parameter in the strong lensing analysis. The two models we assumed to be the most reliable ones in B23 predicted, in the point-like approximation and in the extended source reconstruction, best-fit redshift values of, respectively, 1.72 and 1.69. They also showed that the value of z_C was degenerate with those of other parameters, mainly the Einstein radius (θ_E) and the 3D logarithmic total mass density slope (γ′) of the deflector. We extracted the spectra of C1 and C2 from the apertures shown in Fig. 2, and measured a redshift value of 1.698, from the detection of the MgIIλ2800 at about 7550 Å, with some evidence of a P-Cygni line profile of the doublet, five absorption features of FeII (Fe IIλ2344, 2374, 2382, 2586, 2600), and the emission lines of CIII]λλ1907, 1909 and [OII]λ2470, shown in Fig. 4.

Fig. 4. MUSE 1D spectrum of the multiple images C1 and C2 (stacked), extracted within the green apertures highlighted in Fig. 2. The spectrum was smoothed with a ∼2.5 Å boxcar filter. We zoom in around three wavelength ranges where we observe the main absorption and emission lines used to confirm the redshift measurement of 1.698. The green line represents the observed spectra in units of 10−18 erg s−1 cm−2 Å−1, and the gray regions indicate the 1σ uncertainties.

3.3. Multiple images of E and F

With MUSE, we identify two additional multiply imaged background sources, labeled as E and F in Fig. 1, which do not show any stellar continuum counterpart in HST F160W. They were detected through a strong emission line at ≈7231 Å, as it often happens in deep MUSE observations of similar (e.g., Collett & Smith 2020) and blank (e.g., Bacon et al. 2023) fields. Based on the asymmetric line profile with a clear blue cutoff, this line is interpreted as the Lyα line at z = 4.95. E and F have four multiple images each, whose spectrum around the observed emission line is shown in Fig. 5. Their position coordinates are reported in Table 1.

Fig. 5. MUSE Lyα 2D contours and 1D spectra of the multiple images of families E and F. Left: Contours from a continuum-subtracted narrow band image, centered on the emission line detected at ≈7231 Å, identified as the Lyα emission at z = 4.95. The contours are superimposed on the HST F160W image, to show that continuum counterparts are not detected. Center: Zoom-in of the MUSE 1D spectra of the multiple images E1, E2, E3, E4, and stacked (thicker line), around the detected emission line. Right: Same as in the center for F1, F2, F3, and F4. The spectra were extracted within the orange (E) and brown (F) apertures in Fig. 2, consisting of circular apertures with a radius of 0.8″, and were smoothed with a ∼2.5 Å boxcar filter.

3.4. Discussion of system D

We identify an additional system with four multiple images, labeled as D1-4 in Fig. 1. D1 and D4 are distorted into an extended arc southeast of the main elliptical deflector, D2 is faint in the northeast direction, and D3 lies between B3 and F3 in the western direction, angularly close to a group member with a F160W magnitude of 21.32 mag. We extracted the D spectrum from an elongated region that includes the D1/D4 arc, shown in Fig. 2. We also extracted the spectra of D2 and D3 within circular apertures, but they revealed to be too faint and too contaminated, respectively. We do not detect any clear spectral feature that allows us to obtain a secure redshift measurement for D. The lack of strong emission lines in MUSE and the positions of the multiple images suggest that D lies at a redshift between approximately 1.5 and 2, but further observations will be needed to confirm the nature of this system. System D was also hypothesized in our previous study, but did not enter into the analysis. We anticipate here that we shall continue to exclude the multiple images of source D in the strong lensing modeling described in the following for different reasons. The lack of a secure redshift measurement, the fact that it is not possible to clearly identify a position and a brightest pixel for image D2, and the position of D3 that is strongly perturbed by the mass distribution of a group member angularly very close, make the inclusion of family D strongly uncertain. Furthermore, our simple best-fit strong lensing model, detailed in Sect. 4, which can reproduce very well the positions of the spectroscopically confirmed multiple images, currently disfavors the inclusion of this system. To improve on that, we would need to include additional mass components that are not fully justified by the insecure nature of D. For these reasons, we leave the inclusion of system D to future strong lensing modeling efforts, should its redshift and its spatial positions be confirmed with deeper observations.

3.5. Group members

In B23, given the large estimated value of the deflector total mass, we considered the possibility of SDSS J0100+1818 being a group of galaxies, and we found 53 galaxies in a 2.5′×2.5′ field of view with photometric redshifts consistent, at the 2σ level, with the spectroscopic redshift of the main elliptical galaxy. Thanks to the MUSE data, we securely select and consider as part of the same group all the galaxies with spectroscopic redshift in the 0.581 ± 0.011 range, corresponding to the mean redshift of the overdensity at z ≈ 0.6, with a spread of ±2000 km s⁻¹ (rest-frame). With this method, we select 19 group members (with QF ≥ 3), including the main elliptical galaxy, hereafter referred to as the brightest group galaxy (BGG), whose positions are shown in Fig. 2. They are distributed through the entire MUSE field of view. The closest and the furthest group members lie approximately 4″ (27 kpc) and 33″ (220 kpc) away in projection from the BGG, respectively. Three group members are located within 7″ (approximately the physical Einstein radius of system AB), eight within 15″, and thirteen within 20″. We measure their F160W Kron magnitudes with SExtractor, ranging from 17.05 mag to 22.96 mag. We note that the second brightest member has a magnitude of 19.94, almost three magnitudes fainter than the BGG.

3.6. Velocity dispersion profile

With the new MUSE data, we are able not only to confirm the extremely large stellar velocity dispersion value, σ_⋆, which distinguishes the main deflector as one of the most massive galaxies known at intermediate redshift, but also to measure a stellar velocity dispersion profile, considering the first aperture of 0.4″ of radius and then different annuli with the same center (see Fig. 6). The spectral S/N is larger than 15 in all the selected bins and the corresponding velocity dispersion values are partially correlated, given the observational PSF FWHM. We measured the value of the stellar line-of-sight velocity dispersion of each bin following the procedure presented in Granata et al. (2023). We used pPXF (penalized pixel-fitting, Cappellari & Emsellem 2004; Cappellari 2017, 2023) to perform a full-spectrum fit, comparing the observed spectra with a combination of stellar templates. These were chosen from a set of 463 UVB stellar spectra with S/N > 100 Å⁻¹ from the X-Shooter Spectral Library (XSL) DR2 (Gonneau et al. 2020), degraded to the instrumental resolution of MUSE, and then combined with additive 12th-degree Fourier polynomials and convolved with a Gaussian line-of-sight velocity distribution. The obtained global stellar velocity dispersion presents a value that is consistent with, but systematically lower than, the previous estimate of (451 ± 37) km s⁻¹ with X-Shooter B23. The difference can be explained mainly by the different S/N of the two datasets (significantly higher for MUSE). In the following, we shall refer to the MUSE estimate, given that it is more robust and will moreover be consistent with the BGG mass obtained from the strong lensing analysis (see Sects. 4 and 5).

Fig. 6. Stellar velocity dispersion profile of the main deflector. The apertures (the first 0.4″-radius aperture and then different annuli with the same center) used to extract the spectra relative to each bin are shown in the inset. In the bottom panel, we show the spectral S/N of each bin, which is larger than 15 in all of them, peaking at around 60. Given the value of the observational PSF FWHM and of the chosen bin widths, the σ⋆ values of the profile are partially correlated.

4. Lensing modeling

We performed our strong lensing modeling with the Gravitational Lens Efficient Explorer (GLEE; Suyu & Halkola 2010) software, which allows one both to optimize the parameter values and to infer their probability distributions, through Bayesian analyses like simulated annealing and Markov chain Monte Carlo (MCMC), respectively. GLEE supports several types of mass and light profiles, whose parameters can be optimized by reproducing the observed positions, fluxes, and/or time delays of multiple images. In this study, we employed the positions of the 18 multiple images from the five background sources A, B, C, E, and F, listed in Table 1. For A, B, and C, we considered each multiple image position as that of its brightest pixel in the HST F160W image, with an uncertainty of one HST pixel (0.066″). For E and F, which do not show continuum in the HST images, we considered the brightest pixel in a narrow-band image centered around the Lyα emission line. For these images, we adopted a positional uncertainty of 0.15″. In this work, we only performed a point-like source modeling, given that the images of the sources E and F do not present any HST continuum. We leave the extended source modeling approach to a future work, in which we shall properly treat the weighing of a combination of extended and point-like sources.

The best-fit values of the parameters of each model were estimated with a simulated annealing technique, structured in several steps, first minimizing the χ² value on the source plane, and then on the deflector plane. We assumed different total mass profiles for the deflector, and evaluated their goodness by considering different statistical estimators, often employed in similar strong lensing studies (see, e.g., Acebron et al. 2017; Mahler et al. 2018). We took into account: (1) the number of degrees of freedom (ndof), defined as the difference between the number of observables and the number of free parameters, (2) the value of the minimum χ², defined as

$$\begin{array}{c}\hfill {\chi }^2={\displaystyle \sum _{i=1}^N}\frac{{|{\mathit{\theta }}_i^{\mathrm{obs}}-{\mathit{\theta }}_i^{\mathrm{pred}}|}^2}{{\sigma }_i^2},\end{array}$$(1)

where θ_i^obs are the positions of the observed images, θ_i^pred are their positions predicted by the model, N is the total number of multiple images, and σ_i is the positional uncertainty relative to the i-th image, (3) the value of the root-mean-square difference (rms) between the observed and the model predicted positions of the multiple images, defined as

$$\begin{array}{c}\hfill \mathrm{rms}=\sqrt{\frac{1}{N}{\displaystyle \sum _{i=1}^N}{|{\mathit{\theta }}_i^{\mathrm{obs}}-{\mathit{\theta }}_i^{\mathrm{pred}}|}^2},\end{array}$$(2)

(4) the Bayesian information criterion (BIC, Schwarz 1978), given by

$$\begin{array}{c}\hfill \mathrm{BIC}=kln{N}_{\mathrm{obs}}+{\chi }^2,\end{array}$$(3)

where k is the number of free parameters and N_obs is the number of data points (=2N, for the x and y coordinates of the multiple images), and (5) the corrected Akaike information criterion (AICc, Akaike 1974; Cavanaugh 1997), defined as

$$\begin{array}{c}\hfill \mathrm{AICc}=2k+{\chi }^2+\frac{2k(k+1)}{N_{\mathrm{obs}}-k-1}·\end{array}$$(4)

The BIC and AICc estimators penalize models with an increasing number of free parameters, to contrast with overfitting. Thus, models with lower BIC and AICc (as well as rms and minimum χ²) values are preferred.

Additionally, we estimated median values and uncertainties for the model parameters from MCMCs of 10⁷ steps, with acceptance rates between 20% and 30%, and rejecting the first 10% burn-in steps. These chains are the final step of a sequence in which each intermediate chain is used to estimate the covariance matrix of the model parameters and to extract the starting point for the following one. To obtain sensible confidence intervals for the values of the parameters from the chains, we rescaled the adopted positional uncertainties until the χ² value was similar to the number of degrees of freedom.

Thanks to the relatively simple structure of group-scale systems, we described the total mass distribution of SDSS J0100+1818 and disentangled the different mass contributions by assuming a straightforward composite model, with different mass components describing those of the BGG, of the group members, and of the extended DM halo. In the following subsections, we detail the different parametrizations considered, which give rise to the four models presented in Sect. 4.3.

4.1. Inclusion of group members

We modeled the dimensionless surface mass density (convergence, κ) of each group member as a dual pseudoisothermal elliptical mass distribution (dPIE; Elíasdóttir et al. 2007). In GLEE, this profile is described by six parameters: the x and y coordinates of the center, the semiminor (b) to semimajor (a) axis ratio q = b/a, the position angle, θ (measured counterclockwise from +x), the Einstein radius, θ_E, the core radius, r_core, and the truncation radius, r_T. Throughout this work, we have assumed the total mass distribution of the members to be spherical (q = 1), and without a core radius (r_core = 0). Within these approximations, the convergence becomes

$$\begin{array}{c}\hfill {\kappa }_{\mathrm{dPIE}}(x,y)=\frac{{\theta }_{\mathrm{E}}}{2}(\frac{1}{R}-\frac{1}{\sqrt{R^2+r_{\mathrm{T}}^2}}),\end{array}$$(5)

where $R=\sqrt{x^2+{\mathit{y}}^2}$, and

$$\begin{array}{c}\hfill {\theta }_{\mathrm{E}}=\frac{4\pi {\sigma }^2}{c^2},\end{array}$$(6)

with σ representing the effective velocity dispersion of the deflector. In our analysis, we imposed the total mass of each group galaxy based on its luminosity, L, as is commonly done in cluster-scale strong lensing modeling (e.g., Grillo et al. 2015; Caminha et al. 2019; Acebron et al. 2022; Bergamini et al. 2023b). We estimated L from the Kron magnitude measured with SExtractor in the HST F160W band, and we linked the i-th galaxy total mass, described by the values of the Einstein radius, θ_E,  i, and truncation radius, r_T,  i, and its luminosity via the following scaling relations:

$$\begin{array}{c}\hfill {\theta }_{\mathrm{E},i}={\theta }_{\mathrm{E},\mathrm{ref}}{\left(\frac{L_i}{L_{\mathrm{ref}}}\right)}^{\alpha }\mathrm{and}{r}_{\mathrm{T},i}=r_{\mathrm{T},\mathrm{ref}}{\left(\frac{L_i}{L_{\mathrm{ref}}}\right)}^{\beta },\end{array}$$(7)

where θ_E,  ref and r_T,  ref are reference values relative to a galaxy at the redshift of the group. In particular, we selected the BGG as the reference galaxy, when it is included in the scaling relations, and the second brightest group member at (RA, Dec) = (1:00:48.68, +18:18:12.67) (m = 19.94 mag), when the BGG is parametrized separately; that is, without its parameter values following the scaling relations. In this latter case, the models will be labeled with “+BGG” at the end of the name.

We fixed α = 0.7 and β = 0.5 (e.g., Grillo et al. 2015; Wang et al. 2022). Given that for a dPIE profile the total mass, M_T, is proportional to θ_E and r_T, the total mass-to-light ratio, for the assumed values of α and β, increases with the luminosity as

$$\begin{array}{c}\hfill \frac{M_{\mathrm{T}}}{L}\sim \frac{{\theta }_{\mathrm{E}}{r}_{\mathrm{T}}}{L}\sim \frac{L^{0.7}{L}^{0.5}}{L}\sim L^{0.2},\end{array}$$(8)

corresponding to the “tilt of the fundamental plane” observed in elliptical galaxies (Faber et al. 1987; Bender et al. 1992; Ciotti et al. 1996; Bernardi et al. 2003; Grillo & Gobat 2010).

4.2. Dark matter halo

We completed the total mass modeling of the group with the inclusion of an additional mass component, representing, in particular, the group-scale DM halo and all the contributions not associated with the member galaxies. We explored two different mass profiles describing the extended and smooth mass distribution of the group. The first one is a pseudo-isothermal elliptical mass distribution (PIEMD; Kneib et al. 1996). In GLEE, it is defined by six parameters: the x and y coordinates of the center, the semi-minor (b) to semi-major (a) axis ratio, q = b/a, the position angle, θ (measured counterclockwise from +x), the Einstein radius, θ_E, and the core radius, r_core. We stress that the value of θ_E is defined for a source at z = ∞ and does not correspond to that of the Einstein radius of the system, which should be nearly independent of the mass modeling details. The value of θ_E is a parameter that describes the lens strength and enters the dimensionless surface mass density, κ_PIEMD, as

$$\begin{array}{c}\hfill {\kappa }_{\mathrm{PIEMD}}(x,y)=\frac{{\theta }_{\mathrm{E}}}{2\sqrt{r_{\mathrm{core}}^2+\frac{x^2}{{(1+e)}^2}+\frac{y^2}{{(1-e)}^2}}},\end{array}$$(9)

where the ellipticity is $e=\frac{1+q}{1-q}$. The second profile is a singular power law elliptical mass distribution (SPEMD; Barkana 1998). In GLEE, it is characterized by seven parameters: the first six are in common with the PIEMD profile, and the seventh is the slope, g, which is related to the three-dimensional logarithmic mass density slope γ′= − d log[ρ(r)]/d log(r) (i.e., ρ ∝ r^−γ′) through γ′ = 2g + 1 (i.e., an isothermal profile corresponds to γ′ = 2 and g = 0.5). In the following, we refer to the physical parameter γ′. Similarly to the PIEMD case, θ_E is a parameter of the mass distribution introduced in the dimensionless surface mass density, κ_SPEMD, as

$$\begin{array}{c}\hfill {\kappa }_{\mathrm{SPEMD}}(x,y)={\theta }_{\mathrm{E}}{(x^2+\frac{y^2}{q^2}+\frac{4r_{\mathrm{core}}^2}{{(1+q)}^2})}^{-\frac{{\gamma }^{\prime }-1}{2}}.\end{array}$$(10)

Following the definitions implemented in GLEE, the values of the θ_E parameters adopted in the PIEMD and in the SPEMD cases differ by a factor of (1 + q).

4.3. Mass models and results

We tested four different mass models, exploring two alternative profiles for the smooth, extended group mass distribution (PIEMD and SPEMD) and two options for the BGG in terms of scaling relations. Given the nature of the group, with the BGG being almost three magnitudes brighter than the second brightest galaxy and likely including a large fraction of the total mass of the group members, we tried two different options for the BGG: one in which the BGG follows the scaling relations of Eq. (7), and another one in which the total mass of the BGG is still parametrized by a dPIE profile, but independently of the other members. This second option allows the model to associate with more freedom the total mass contribution related to the BGG and to the other group members. To employ the lowest number of free parameters, we fixed the BGG position to its luminosity center, q = 1, and r_core = 0, leaving its total mass described only by θ_E and r_T. The four considered models are thus:

– PIEMD: This model includes a PIEMD mass distribution for the extended halo, and dPIEs for all the group members. The BGG is included in the scaling relations, so that the total mass of each group member is related to that of the BGG. This mass model is described by 8 free parameters: 6 relative to the extended halo and 2 to the scaling relations. Given that for all the models the number of observables is equal to 36 (x and y coordinates of each of the 18 total multiple images), this model has ndof = 18.

– SPEMD: Similar to the PIEMD model, this model includes a SPEMD mass distribution for the extended halo, and dPIEs for the group members, including the BGG. It is described by 9 free parameters (adding γ′), and has ndof = 17.

– PIEMD+BGG: This model includes a PIEMD mass distribution for the extended halo, an independent dPIE for the BGG, and dPIEs for the group members, linked through the scaling relations. This mass model is described by 10 free parameters: 6 relative to the extended halo, 2 relative to the dPIE of the BGG, and 2 to the scaling relations, and has ndof = 16.

– SPEMD+BGG: This model considers a SPEMD mass distribution for the extended halo, an independent dPIE for the BGG, and dPIEs for the group members linked through the scaling relations. It is described by 11 free parameters and has ndof = 15.

The best-fit values of the parameters corresponding to each model are reported in Table 2. By comparing the best-fit parameter values of these different models, we noticed that the center of the total mass approximately coincides with the luminosity center, when the BGG is modeled inside the scaling relations, while it is shifted by ≈0.2″ in the −x-direction, when it is not. The value of the axis ratio, q, of the extended halo does not vary significantly among all the different models; nor does that of the position angle, θ. The PIEMD and SPEMD models predict similar values of the Einstein radii, which describe the strength of the lens, both for the extended halo and the BGG (θ_E,  ref). This latter value of approximately 4.7″ is also found when the BGG is modeled outside of the scaling relations (in the PIEMD+BGG and SPEMD+BGG models), strengthening its robustness. The core radius of the diffuse halo increases by ≈15%, when it is parametrized with a SPEMD mass distribution, whose best-fit slope (γ′∼2.2) is slightly steeper than isothermal. The BGG has a truncation radius (r_T,  ref in PIEMD and SPEMD, and r_T in PIEMD+BGG and SPEMD+BGG) larger than the distance over which the multiple images appear (see Table 1).

In summary, we have found that the best-fit values of the parameters are stable when the BGG is modeled outside of the scaling relations. We note that, by reducing the number of degrees of freedom from 18 of the PIEMD model to 15 of the SPEMD+BGG one, the χ_min² value decreases from 86.63 to 76.82, with a corresponding reduction in the rms value from 0.23″ to 0.21″.

4.4. Mass profiles of the deflectors

We measured the cumulative total mass (M_T) profile for each model by randomly extracting 1000 parameter value sets from the last MCMC chains of 10⁷ steps described above. GLEE can create convergence maps, which we converted into total mass maps. Then, we summed up the contribution of all the pixels within circular apertures centered on the brightest pixel of the main lens galaxy, with a step of 0.5 pixels, to obtain the cumulative total mass profiles that are presented in the following. For each aperture, we estimated the 1σ uncertainties as the 16^th and 84^th percentile values of the distribution of the total mass values of all the 1000 random models. We separated the total mass contribution of the extended halo (PIEMD or SPEMD) from that of the group members (dPIE), by considering the relevant free parameters in the MCMC chains. The resulting cumulative mass profiles are shown in Fig. 7. The four models predict fully consistent total mass profiles, with mean relative uncertainties of only 2%. The total mass value projected within 50 kpc, approximately equivalent to the Einstein radius of the system, is (1.55 ± 0.01)×10¹³ M_⊙, while within the projected distance of the furthest multiple image (R ≈ 77 kpc) is (2.78 ± 0.04)×10¹³ M_⊙.

Fig. 7. Cumulative projected total mass profiles for the PIEMD (solid black), SPEMD (solid blue), PIEMD+BGG (dashed black), and SPEMD+BGG (dashed blue) models with ±1σ uncertainties (differently hatched shaded areas). The curves on top, labeled with MT, are relative to the total mass, and they are divided into the halo component (Mhalo, left panel) and the galaxy member component (Mgal, right panel), identified by the same linestyle. For comparison, we leave Mgal and Mhalo as faint lines in the left and right panels, respectively. The vertical lines close to the x axis locate the distances from the lens galaxy center of the different multiple images, color-coded following Fig. 1. The black arrow at R = 9.32 kpc shows the effective radius of the BGG.

The mass component associated with the diffuse halo presents a slightly larger mean uncertainty of around 10% in all the four models. For this component, the models predict a projected total mass value of $9.2_{-3.1}^{+1.2}\times {10}^{12}{\mathrm{M}}_{\odot }$, at R = 50 kpc, and of $1.9_{-0.7}^{+0.2}\times {10}^{13}{\mathrm{M}}_{\odot }$, at R = 77 kpc. As is visible in Fig. 7, the profiles show that, in the central regions, the total mass is dominated by the group member contribution, while from R ≈ 35 kpc outward the extended halo prevails. The different models predict consistent relative weights for the different mass components, with the models with the BGG outside of the scaling relations assigning to the extended halo a slightly lower contribution.

We measured the stellar mass (M_*) profile of the BGG in B23 from its reconstructed luminosity profile, by assuming a constant stellar mass-to-light ratio estimated through a spectral energy distribution (SED) fitting. We assume here the same value for the stellar mass-to-light ratio of the other galaxy members and convert the HST F160W image into a M_* map, by multiplying it by the BGG mass-to-light ratio and by a mask designed to include only the pixels associated with the 19 group members (including the BGG). This is justified by the fact that all the group members show similar spectra and colors in the MUSE and photometric data, and thus they are likely early-type galaxies hosting similar stellar populations. From the resulting image, we measured the cumulative projected stellar mass profile of the group by summing the contribution of the pixels within circular apertures centered on the brightest pixel of the main lens galaxy, with a step of 0.5 pixels. The profile is shown in the top panel of Fig. 8. The BGG clearly dominates the light (and thus the stellar mass) distribution, and the inclusion in the external apertures of the other group members is visible as steps in the profile. Given that we can measure the stellar mass distribution of the group all across the MUSE field of view, it is possible to explore the outer regions, up to R ≈ 30″ ≈ 180 kpc. We note that the total stellar mass value of the group at R = 180 kpc is approximately equal to 1.5 times that of the BGG and, at the largest distance where the multiple images appear (R ≈ 77 kpc), the total stellar mass of the group is (1.7 ± 0.3)×10¹² M_⊙, and the BGG contributes with almost 90% of it, namely with (1.5 ± 0.3)×10¹² M_⊙. We remark that the cumulative stellar mass profile is between 1.5 and 5 times smaller than the total one associated with the galaxy member components, shown in the right panel of Fig. 7. This difference is mainly due to the presence of DM, in the form of galaxy-scale halos, in the group members.

Fig. 8. Cumulative projected total mass and stellar-over-total mass fraction profiles. Top: Cumulative projected total mass profiles for the PIEMD (solid black), SPEMD (solid blue), PIEMD+BGG (dashed black), and SPEMD+BGG (dashed blue) models with ±1σ uncertainties (differently hatched shaded areas), as in Fig. 7, compared with the cumulative projected stellar mass profile with ±1σ uncertainties, represented by the dashed light blue curve and shaded area, respectively. Bottom: Cumulative projected stellar-over-total mass fraction profiles, color-coded for each model according to the top panel. In both panels, the vertical lines close to the x axis locate the distances from the lens galaxy center of the different multiple images, color-coded following Fig. 1. The black arrow at R = 9.32 kpc shows the effective radius of the BGG.

In the bottom panel of Fig. 8, we plot the cumulative projected stellar-over-total mass fraction profiles relative to each model. They are all in very good agreement, differing only by less than 2% in the inner region, and by a fraction of percent from R > 20 kpc outward. We measured a stellar-over-total mass fraction value of $45.6_{-8.3}^{+8.7}\%$ at the lens galaxy effective radius (R_e = 9.32 kpc), decreasing to (11.1 ± 1.9)% at R = 50 kpc and to (6.6 ± 1.1)% at R ≈ 77 kpc.

5. Strong lensing models with variable cosmology

5.1. Single-plane models

The observed positions of the lensed multiple images of a system depend both on the total mass distribution of the lens and on the cosmological angular-diameter distances between the observer, the lens, and the sources. Systems where multiple sources at different redshifts are strongly lensed by the same foreground deflector allow one to break the mass-cosmology degeneracy, and thus to learn about the geometry of the Universe. By measuring the ratio of different Einstein radii in a multiple-source system, and assuming a cosmological model, it is possible to estimate the values of the cosmological parameters defining the global geometry of the Universe, without any dependence on the value of H₀ (Collett et al. 2012). The same parameter values can be measured thanks to other cosmological probes (Eisenstein et al. 2005; Komatsu et al. 2014; Treu et al. 2022). In light of the results of the statistical estimators adopted and discussed in the previous sections, we decided to consider in the following only the PIEMD and SPEMD models. We assumed a flat (Ω_Λ + Ω_m = 1) ΛCDM cosmology and reoptimized the PIEMD_Ω_m and SPEMD_Ω_m models, where Ω_m was introduced as a free parameter. We further reconsidered the PIEMD_Ω_mw and SPEMD_Ω_mw models in a flat wCDM cosmology, where both Ω_m and w were free to vary. We reran the same procedure described above, starting from a χ² minimization on the lens plane and then obtaining a sequence of MCMC chains, resulting in a final one composed of 10⁷ steps. We adopted uniform priors on the values of Ω_m and w, respectively, between 0 and 1, and −2 and 0.

To fully exploit our observational dataset and alleviate possible degeneracies bewteen the values of the cosmological and lens total mass distribution parameters, we used the available kinematic information on the BGG to impose a Gaussian prior on the value of its Einstein radius θ_E,  ref. As we wished to probe the inner kinematics of the BGG, we weighed the MUSE cube with its observed surface brightness in the HST F160W band, degraded and re-binned to the PSF and pixel-scale of MUSE. We then extracted a spectrum within an aperture with a radius of 3″, centered on the BGG. This light-weighed extraction procedure yields measured velocity dispersion values equivalent to those that would be obtained from a non-weighed spectrum extracted within the galaxy effective radius (Granata et al., in prep.), while increasing the S/N of the extracted spectrum compared to a small straightforward aperture. By doing so, we measured a stellar velocity dispersion value of (380.5 ± 4.4) km s⁻¹. We converted this value, with a conservative uncertainty of 10 km s⁻¹, to the corresponding Einstein radius for an isothermal profile, through Eq. (6). We note that the value obtained in this way is very similar to the best-fit values of θ_E,  ref and θ_E of the BGG reported in Table 2, as well as consistent within about 1σ with the results from the respective posterior probability distributions (e.g., (409 ± 29) km s⁻¹ for the PIEMD+BGG model). Furthermore, we checked that the introduction of this prior does not have any significant impact on the lens mass profiles reconstructed through the different models. In detail, none of the profiles shown in Figs. 7 and 8 vary appreciably.

We report in Table 3 the median values, with the 68% confidence level (CL) uncertainties, of the cosmological parameters Ω_m and w, and show in Fig. 9 their posterior probability distributions for the PIEMD_Ω_mw and SPEMD_Ω_mw models. The average 1σ uncertainty is on the order of 0.1 for Ω_m and 0.4 for w. The Ω_m distributions are skewed toward the lower limit of our uniform prior (i.e., zero), and thus in the table we also include the 95.4% quantile upper limits. The models with w as a free parameter predict median values of w smaller than −1, with relative uncertainties of approximately 30%. In particular, we found $\mathit{w}=-1.{27}_{-0.48}^{+0.43}$ and $\mathit{w}=-1.{38}_{-0.41}^{+0.38}$ for the PIEMD_Ω_mw and SPEMD_Ω_mw models, respectively. Given that the w distributions decrease quite slowly toward small w values, the resulting median values and lower limits are sensitive to the extent of the (flat) prior. We remark that our results were obtained with a flat and quite large prior (i.e., w ∈ [ − 2, 0]). We additionally report here the more robust estimates of the 95.4% quantile upper limits for w, of −0.53 and −0.79 for the PIEMD_Ω_mw and SPEMD_Ω_mw models, respectively. The considered statistical estimators did not reveal a strongly preferred model, as all of them reproduced the observed positions of the multiple images with the same rms of 0.23″ (and thus had very similar χ_min² values). Consequently, the PIEMD_Ω_m model, given its smallest number of free parameters, is moderately preferred by both the BIC and AICc metrics, followed by the PIEMD_Ω_mw model.

Fig. 9. Probability density distributions of the cosmological parameters Ωm and w in the PIEMD_Ωmw (black) and SPEMD_Ωmw (blue) models. We show the marginalized 1D histograms of each parameter and their joint 2D probability distribution. The vertical dashed lines in the 1D histograms represent the 16th, 50th, and 84th percentiles, while the solid lines in the 2D distributions represent the 0.68, 0.95 and 0.99 contour levels.

5.2. Multiplane models

Given the presence of galaxies (i.e., mass concentrations) at different redshifts, we additionally adopted a multiplane strong lensing approach (Blandford & Narayan 1986; Schneider et al. 1992), to take into account the fact that the light rays coming from the furthest sources are deflected not only once, by the main group at z = 0.581, but instead multiple times by other secondary deflectors at different redshifts. In particular, we explored the impact of such an approach on the inferred values of the cosmological parameters.

In the multiplane formalism, the total deflection angle is a suitable sum of all the deflection angles relative to the different deflector planes (see, e.g., Gavazzi et al. 2008). This approach has been used to model both different foreground and background deflectors at distances different from that of the primary lens (e.g., Chirivì et al. 2018; Schuldt et al. 2024; Acebron et al., in prep.) and to assign mass to multiply imaged background sources, where several sources at different redshifts are lensed by the same deflector (e.g., Gavazzi et al. 2008; Schuldt et al. 2019; Collett & Smith 2020; Wang et al. 2022).

The lens models described above predict that the sources A, B, E, and F are located angularly close to the optical axis (i.e., the center of the BGG), while C lies ∼2″ away. Considering the mutual angular distances, and the pixelated source reconstructions developed in B23, it is likely that the second largest contribution to the total deflection is that happening at the plane of A and B on the light rays coming from sources E and F. Thus, we developed the MP-PIEMD_Ω_m and MP-PIEMD_Ω_mw models, which extended the PIEMD_Ω_m and PIEMD_Ω_mw models, respectively, by adding an additional PIEMD mass distribution at the redshift of A and B, z = 1.880. We imposed it to be spherical (q = 1) and with a vanishing core radius. Thus, this mass structure is described by three parameters: the coordinates (x₂, y₂) of the center of mass and θ_E2, where the subscripts indicate that these parameters are referred to the secondary lens. We repeated the procedure described above, minimizing the χ² on the lens plane, and then running a set of MCMC chains, resulting in a final one of 10⁷ steps. The value of θ_E2 was extracted with a uniform large prior, while the values of x₂ and y₂ were linked to the position of source A (that is brighter and more extended than B). The MP-PIEMD_Ω_m and MP-PIEMD_Ω_mw models are able to reproduce the observed positions of the multiple images with rms values of 0.23″ (as for the PIEMD_Ω_m and PIEMD_Ω_mw models), but they are penalized by their BIC and AICc scores, due to the inclusion of the additional free parameter (the position of source A was optimized also in the previous models). From the posterior probability distributions, we observed a clear degeneracy between the parameters describing the centroids of the main halo at z = 0.581, (x, y), and of the secondary lens, (x₂, y₂). No clear degeneracy is observed with the cosmological parameters, which present consistent median values and uncertainties similar for w and boosted by a factor of ∼1.5 for Ω_m, due to the more extended tails toward large values (see Fig. 10 and the bottom rows of Table 3). We also notice that these extended models predict cumulative total mass profiles for the main deflector that are fully consistent with those shown in Fig. 7.

Fig. 10. Marginalized 1D probability density distributions of the cosmological parameters Ωm and w in the PIEMD_Ωm (black), SPEMD_Ωm (blue), and MP-PIEMD_Ωm (red) models (top row), and PIEMD_Ωmw (black), SPEMD_Ωmw (blue), and MP-PIEMD_Ωmw (red) models (bottom row). The vertical dashed lines represent the 16th, 50th, and 84th percentiles for each distribution.

6. Discussion

6.1. An enhanced model for SDSS J0100+1818

We compare here our results with those obtained in B23, highlighting the improvements and the overall consistency. The new MUSE data allowed us to:

– spectroscopically confirm the redshift of source C (z_C = 1.698). This value is fully consistent with the best-fit measurements of 1.72 and 1.69 by B23. This spectroscopic confirmation enabled us to reduce some degeneracies found in the previous work between z_C and the parameters (in particular θ_E and γ′) describing the total mass distribution of the deflector;

– characterize the deflector plane. We identified 18 group members (excluding the BGG) that were included in the strong lensing model as individual mass components, instead of as an external shear term, and we measured the BGG stellar velocity dispersion profile with an accuracy of ≲8% (approximately 3% in the four inner bins);

– reduce the relative uncertainties on the cumulative total mass profiles from 17%, 1%, and 5%, to 3%, < 1%, and 1%, at R ≈ 15, 42, and 63 kpc, respectively, thanks to the increased number of multiple images from additional sources at different redshifts;

– enlarge the radial interval probed by strong gravitational lensing with the discovery of the systems E and F, from the closest multiple image (C2, at R ≈ 15 kpc) to the furthest (E3 at R ≈ 77 kpc, previously C1 at R ≈ 63 kpc). This allowed us to robustly reconstruct the total mass profile at previously unexplored physical distances from the group center with a precision of just a few percent;

– confirm the measured median and ±1σ cumulative total mass profiles (Fig. 7), as well as the stellar-over-total mass profiles (Fig. 8), found (after correction) by B23 over the explored radial distances. In particular, we measure (B23 measured) a projected total mass value of (1.20 ± 0.01)×10¹³ M_⊙ ((1.16 ± 0.01)×10¹³ M_⊙) within 42 kpc, and a stellar-over-total mass fraction of $45.6_{-8.3}^{+8.7}\%$ at the lens galaxy effective radius, decreasing to (13.6 ± 2.3)% at R = 42 kpc ((38 ± 9)% decreasing to (12 ± 2)%).

6.2. SDSS J0100+1818 as a fossil system

The hierarchical evolution of structures predicted by the ΛCDM scenario implies that, in the case of some compact groups of galaxies, all the main galaxies of a group could merge into a single one (Ponman & Bertram 1993). As this central galaxy cannibalizes more galaxies, the luminosity gap between it and the remaining galaxies increases with time. Jones et al. (2003) firstly gave an operational definition of the so-called “fossil” groups, namely systems where the magnitude gap in the r-band between the first two brightest galaxies, Δm₁₂, lying at a projected distance smaller than half the value of the virial radius, is larger than 2. Later, Dariush et al. (2010) proposed to consider the magnitude gap between the first and the fourth brightest galaxies, requiring Δm₁₄ > 2.5. This criterion was shown to be more robust in classifying the groups with a small number of members, where the serendipitous infall of a bright galaxy would change the fossil status. Simulations show that the fossil system status is a phase in the evolution of a group rather than a class itself, as it can be lost if a bright galaxy is included in the system, but also acquired if the BGG merges with other members. Both definitions also include the presence of a diffuse X-ray emission, with L_X ≥ 10⁴² erg s⁻¹. The growing mechanism of fossil groups, given the typical angular momentum loss rates via dynamical friction, occurs over long timescales, and thus the discovery of fossil systems as old, undisturbed groups at intermediate redshift (approximately half of the age of the Universe) can act as a probe for galaxy evolution models. Moreover, fossil systems have a mass concentration in their center that is higher than that of normal galaxy groups (Johnson et al. 2018b), and this makes them particularly efficient as strong gravitational lenses.

Given the dominance of the BGG in its center, the SDSS J0100+1818 group was identified as a candidate fossil system in the CASSOWARY survey, but the lack of observations did not allow one to draw firm conclusions. Even now, we would need additional data to securely assign it the status: future deep X-ray data would be necessary to accurately measure the virial radius and confirm whether SDSS J0100+1818 is a dynamically evolved system immersed in a fairly undisturbed hot gas halo, as is imposed by the criterion. This X-ray condition was originally imposed to guarantee the presence of an extended halo in which the group is immersed, and this is clearly the case from the outcomes of the strong lensing modeling and from the total mass parametrization assumed (see also the studies Grillo et al. 2013; Wang et al. 2022, about CSWA31, a similar system included in the CASSOWARY survey). Additionally, SDSS J0100+1818 fulfills both the Jones et al. (2003) and Dariush et al. (2010) criteria, presenting Δm₁₂ = 2.89 and Δm₁₄ = 3.35. The confirmation of SDSS J0100+1818 as a fossil system would be decisive in putting new constraints on the mechanisms that drive the formation and evolution of fossil systems and of extremely massive early-type galaxies, given that it would be the most distant system known so far (Khosroshahi et al. 2006; Aguerri & Zarattini 2021, and references therein). SDSS J0100+1818 redshift is similar to that of the system studied by Ulmer et al. (2005), which does not satisfy the condition of Jones et al. (2003), with Δm₁₂ = 1.93.

6.3. The values of Ω_m and w from multiplane strong lensing systems

We compare here the cosmological measurements obtained in this work with those from the galaxy-scale multi-source-plane strong lensing system SDSS J0946+1006 (also known as the “Jackpot”; Gavazzi et al. 2008; Smith & Collett 2021). The latter is composed of a galaxy at z_l = 0.222 that acts as a strong lens on three background sources, at redshifts 0.609, 2.035, and 5.975. Even if the Jackpot deflector redshift is lower than that of SDSS J0100+1818, the lensed sources span very similar redshift intervals, with the D_ds/D_s ratios ranging from 0.54 to 0.72 for SDSS J0100+1818 and from 0.60 to 0.89 for the Jackpot (in the standard cosmological model). As is shown in Collett et al. (2012), a double-source system with the first source at a redshift relatively close to that of the deflector and the second one at a significantly higher redshift is particularly effective at estimating the value of w. More generally, different combinations of source redshifts lead to a higher sensitivity in different parts of the Ω_m − w plane, making it particularly important in the future to combine the cosmological inferences from a large sample of multisource strong lensing systems. We find that our measurements of Ω_m and w with SDSS J0100+1818 are comparable, in terms of values and accuracy, to those obtained from the Jackpot in Collett & Smith (2020) and Smith & Collett (2021), where they were combined though with priors from Planck (Planck Collaboration XVI 2014).

Given the group nature of SDSS J0100+1818 and the relatively large number of background sources and their multiple images, we also compare our results with those obtained with lens clusters of galaxies. These systems can rely on a large number of multiple images, from several sources located at different redshifts, which help to break the mass-sheet and mass-cosmology degeneracies, making them excellent cosmological probes. On the other hand, cosmographic measurements are hindered by the complex total mass distribution of clusters: in order to accurately reproduce the observed positions of many multiple images, distributed over wide spatial regions, several mass components, tracing both the baryonic (stars and hot-gas) and DM distributions, are usually needed. This complexity can introduce some systematic uncertainties associated with the specific choice of the total mass parametrization. On the contrary, galaxy- and group-scale systems have the advantage that the lens is usually successfully modeled with a few mass components, with an immediate physical interpretation, like the simple DM halo, BGG, and group member contributions in SDSS J0100+1818.

Caminha et al. (2022) exploited five strong lensing clusters (with redshift from 0.234 to 0.587), with 7−20 multiply imaged sources per system, with redshifts from 0.68 to 6.85. They combined the results from the lens clusters in the sample and found ${\mathrm{\Omega }}_{\mathrm{m}}=0.{30}_{-0.11}^{+0.09}$ and $\mathit{w}=-1.{12}_{-0.32}^{+0.17}$ but, when considering the individual clusters, they obtained 1σ uncertainties on Ω_m ranging from 0.2 to 0.4. This relative uncertainty is fully consistent with the one obtained from SDSS J0100+1818, showing that the lower number of multiple images is balanced by the simpler total mass distribution. Moreover, Caminha et al. (2022) also noticed that the typical “L-shape” they observed in the Ω_m − w plane (also visible for SDSS J0100+1818 in Fig. 9) is less pronounced, as the redshift of the deflector increases. Only one of the clusters they studied is at a redshift larger than that of SDSS J0100+1818. This means that our system allows us to exploit still relatively uncommon deflector and source redshift combinations, and thus to explore different regions in the Ω_m − w plane. This evidence suggests that, together with the comparable accuracy level reached, massive galaxy- and group-scale systems can complement cosmographic studies through strong lensing performed with clusters.

7. Summary and conclusions

In this paper, we have exploited new VLT/MUSE data to expand the study of Bolamperti et al. (2023) (B23) about the strong lensing system SDSS J0100+1818. This system consists of a group of galaxies, dominated by a central massive early-type galaxy (BGG) at redshift z = 0.581, which acts as a strong gravitational deflector for five background sources, all of which have been spectroscopically confirmed. The sources are lensed into a total of 18 multiple images, located between 2.25″ (≈15 kpc) and 11.70″ (≈77 kpc) in projection from the BGG.

We have built, with the software GLEE, different strong lensing models of the deflector, consisting of an extended component (whose total mass distribution has been parametrized with either a PIEMD or a SPEMD profile), representing the DM halo in which the group is immersed, and several dPIE components for the member galaxies. We have measured the best-fit parameter values, and their posterior probability distribution, by reproducing the observed position of all 18 multiple images.

Taking advantage of the multiplane nature of the system, with several sources at different redshifts lensed by the same deflector, we have developed some additional models that also vary the geometry of the Universe. In fact, by considering the distance ratios between the observer, the lens, and the different sources, and assuming a model of the Universe (flat ΛCDM or wCDM models), the values of Ω_m and w can be inferred.

The main results of this work can be summarized as follows:

• We have confirmed the redshift of the background sources A and B, at z = 1.880. We have also measured the redshift of source C, at z = 1.698. Furthermore, we have discovered two additional background sources, E and F, at z = 4.95, which do not show any HST continuum but which do have clear Lyα emission lines.

• We have securely measured the redshifts of 65 objects in the field of view, identifying 19 group members (including the BGG) on the deflector plane, 3 stars, 10 foreground galaxies, 32 background galaxies, and 1 background quasar.

• We have measured the stellar velocity dispersion profile of the BGG over six bins with S/N > 15 up to an aperture of 3″ radius (≈20 kpc, approximately twice the value of the effective radius). The profile has a peak of (386 ± 6) km s⁻¹, and slightly decreases to (350 ± 25) km s⁻¹ in the outer regions.

• We have developed different strong lensing models for the system. We have modeled the projected dimensionless surface mass density of each group member as a dPIE mass distribution, and weighed the total mass of each group galaxy based on its luminosity. We have included an additional mass component, representing all the contributions not associated with the group galaxies. In particular, this main group-scale DM halo has been parametrized with either a PIEMD or a SPEMD total mass distribution. Given the dominance of the BGG, we have also considered models in which the total mass of the BGG was still parametrized by a dPIE profile, but not weighed through the scaling relations adopted for all the other members. The results of these models, summarized in Table 2, show that we could reproduce the observed positions of the 18 multiple images with an rms value between 0.21″ and 0.23″.

• We have accurately measured a total mass profile for the deflector with an average relative uncertainty of only 2% between 15 and 77 kpc. In particular, we have measured a total mass of (1.55 ± 0.01)×10¹³ M_⊙ ((2.78 ± 0.04)×10¹³ M_⊙) within an aperture of R = 50 kpc (R ≈ 77 kpc). The total masses of the extended halo within the same apertures are, respectively, $9.2_{-3.1}^{+1.2}\times {10}^{12}{\mathrm{M}}_{\odot }$ and $1.9_{-0.7}^{+0.2}\times {10}^{13}{\mathrm{M}}_{\odot }$.

• By adopting the stellar mass-to-light ratio estimated in B23, we have measured the group cumulative stellar mass profile up to R ≈ 30″ ≈ 180 kpc. We have employed this profile to obtain the stellar over total mass fraction profiles for each model. They have resulted in a value of $45.6_{-8.3}^{+8.7}\%$ at the BGG effective radius, decreasing to (11.1 ± 1.9)% at R = 50 kpc and to (6.6 ± 1.1)% at R ≈ 77 kpc.

• We have developed strong lensing models with the values of Ω_m and w free to vary, in flat ΛCDM and wCDM models of the Universe. The average 1σ uncertainty is on the order of 0.1 for Ω_m, and of 0.4 for w. We have measured ${\mathrm{\Omega }}_{\mathrm{m}}=0.{14}_{-0.09}^{+0.16}$ (with a 95.4% upper limit of 0.47) in the flat ΛCDM case, and ${\mathrm{\Omega }}_{\mathrm{m}}=0.{19}_{-0.10}^{+0.17}$ (with a 95.4% upper limit of 0.59) and $\mathit{w}=-1.{27}_{-0.48}^{+0.43}$ (with a 95.4% upper limit of −0.53) in the flat wCDM case.

• We have considered the MP-PIEMD_Ω_m and MP-PIEMD_Ω_mw models, by adopting a multiplane strong lensing modeling approach. These models have included an additional PIEMD mass distribution at z = 1.880, representing the total mass term associated with sources A and B, and contributing to the total deflection angle. We have found consistent median values for the cosmological parameters Ω_m and w, with a similar uncertainty on w and an increased uncertainty on Ω_m by a factor of approximately 1.5.

• We have checked that the accuracy of the cosmological measurements obtained with SDSS J0100+1818 is similar to that obtained with the galaxy-scale multi-source-plane strong lensing system known as the “Jackpot” (Collett & Auger 2014; Collett & Smith 2020; Smith & Collett 2021). Given the group nature of SDSS J0100+1818 and the quite large number of multiple images, we have also compared the accuracy of our estimates with those provided by lens clusters of galaxies. By comparing with the sample presented in Caminha et al. (2022), we have found that SDSS J0100+1818 performs as a typical cluster of galaxies, showing that the lower number of multiple images is balanced by the simpler total mass distribution. As an advantage, SDSS J0100+1818 has a lens redshift that is larger than the average one for the lens clusters, and this reduces the degeneracy between the values of Ω_m and w.

The main goals of this study were to make public a new spectroscopic catalog in the SDSS J0100+1818 field of view, obtained from recent MUSE data, with a particular focus on the deflector group members and on the multiply lensed background sources, and to present the capabilities of this system through a detailed strong lensing analysis to characterize the deflector total mass and infer the values of some cosmological parameters. We plan to enhance this study with already available and future observations. For instance, the kinematics of the BGG and of the brightest group members could be included to refine the strong lensing analysis, new models could be developed considering not just the point-like position but also the extended structure of the multiple images, and the background sources could be further characterized, by using their morphology and kinematics. Moreover, the extended radial profile shown in Fig. 6 will allow us to develop a joint strong lensing and dynamics model, to better characterize the total mass distribution in the inner regions of the group. Additional new deep rest-frame infrared spectroscopic observations could allow us to measure the redshift of source D, making it the first group-scale system, dominated by a central BGG, with 6 lensed sources on 4 different redshift planes, and to possibly discover other distant lensed line emitters, like E and F. Furthermore, deep X-ray observations could confirm the fossil nature of SDSS J0100+1818, and enable us to put constraints on the extended halo in which the group is immersed. A detailed study of SDSS J0100+1818 will pave the way for the characterization of tens of such rare galaxy- and group-scale systems that will be employed for cosmography and that will be discovered with new and upcoming facilities like Euclid and Rubin-LSST.

Acknowledgments

This research is based on observations collected at the European Organisation for Astronomical Research in the Southern Hemisphere under ESO programme 0110.A-0248(A). It is also based on observations collected under ESO programme 091.A-0852(A) and on observations made with the NASA/ESA Hubble Space Telescope obtained from the Space Telescope Science Institute, which is operated by the Association of Universities for Research in Astronomy, Inc., under NASA contract NAS 5–26555. These observations are associated with program GO-15253. We acknowledge the use of numpy (Harris et al. 2020), matplotlib (Hunter 2007), astropy (Astropy Collaboration 2022), pandas (The pandas development team 2023), corner.py (Foreman-Mackey 2016), and CARTA (Comrie et al. 2021). A.B. and A.Z. acknowledge support from the INAF minigrant 1.05.23.04.01 “Clumps at cosmological distance: revealing their formation, nature, and evolution”. C.G. and G.G. acknowledge support through grant PRIN-MIUR 2020SKSTHZ. S.H.S. thanks the Max Planck Society for support through the Max Planck Research Group. L.C. is supported by DFF/Independent Research Fund Denmark, grant-ID 2032–00071.

References

All Tables

All Figures

Fig. 1. HST F160W image of the SDSS J0100+1818 strong lensing system studied in this work. Letters label multiple images of the same background source. E and F do not show a continuum in the HST image and the orange and brown crosses mark the positions of their Lyα emission line peaks, detected by MUSE. The dashed gray circle represents an aperture of 50 kpc, the approximate average Einstein radius of the system.

In the text

Fig. 2. Group members (squares) identified as those galaxies whose redshift is within 0.581 ± 0.011, equivalent to the mean redshift of the group with a spread of 2000 km s−1 (rest-frame). The colored regions represent the apertures within which we extracted the spectra of the multiple images of A and B (blue), C (green), D (purple), E (orange), and F (brown). The spectra of E and F, which do not present an HST-detected continuum, are extracted within circular apertures with a radius of 0.8″. The dashed gray curve shows a circle with a radius of 50 kpc, which represents approximately the average Einstein radius of the system.

In the text

Fig. 3. MUSE 1D spectra of the multiple images A1+B1, A2+B2, A3+B3, A4+B4, and stacked (top row), extracted within the blue apertures highlighted in Fig. 2. The spectra were smoothed with a ∼ 2.5 Å boxcar filter. We zoom in around five wavelength ranges where we observe the main absorption and emission lines used to confirm the redshift measurement of 1.880. The blue lines represent the observed spectra in units of 10−18 erg s−1 cm−2 Å−1, and the gray regions indicate the 1σ uncertainties.

In the text

Fig. 4. MUSE 1D spectrum of the multiple images C1 and C2 (stacked), extracted within the green apertures highlighted in Fig. 2. The spectrum was smoothed with a ∼2.5 Å boxcar filter. We zoom in around three wavelength ranges where we observe the main absorption and emission lines used to confirm the redshift measurement of 1.698. The green line represents the observed spectra in units of 10−18 erg s−1 cm−2 Å−1, and the gray regions indicate the 1σ uncertainties.

In the text

Fig. 5. MUSE Lyα 2D contours and 1D spectra of the multiple images of families E and F. Left: Contours from a continuum-subtracted narrow band image, centered on the emission line detected at ≈7231 Å, identified as the Lyα emission at z = 4.95. The contours are superimposed on the HST F160W image, to show that continuum counterparts are not detected. Center: Zoom-in of the MUSE 1D spectra of the multiple images E1, E2, E3, E4, and stacked (thicker line), around the detected emission line. Right: Same as in the center for F1, F2, F3, and F4. The spectra were extracted within the orange (E) and brown (F) apertures in Fig. 2, consisting of circular apertures with a radius of 0.8″, and were smoothed with a ∼2.5 Å boxcar filter.

In the text

Fig. 6. Stellar velocity dispersion profile of the main deflector. The apertures (the first 0.4″-radius aperture and then different annuli with the same center) used to extract the spectra relative to each bin are shown in the inset. In the bottom panel, we show the spectral S/N of each bin, which is larger than 15 in all of them, peaking at around 60. Given the value of the observational PSF FWHM and of the chosen bin widths, the σ⋆ values of the profile are partially correlated.

In the text

Fig. 7. Cumulative projected total mass profiles for the PIEMD (solid black), SPEMD (solid blue), PIEMD+BGG (dashed black), and SPEMD+BGG (dashed blue) models with ±1σ uncertainties (differently hatched shaded areas). The curves on top, labeled with MT, are relative to the total mass, and they are divided into the halo component (Mhalo, left panel) and the galaxy member component (Mgal, right panel), identified by the same linestyle. For comparison, we leave Mgal and Mhalo as faint lines in the left and right panels, respectively. The vertical lines close to the x axis locate the distances from the lens galaxy center of the different multiple images, color-coded following Fig. 1. The black arrow at R = 9.32 kpc shows the effective radius of the BGG.

In the text

Fig. 8. Cumulative projected total mass and stellar-over-total mass fraction profiles. Top: Cumulative projected total mass profiles for the PIEMD (solid black), SPEMD (solid blue), PIEMD+BGG (dashed black), and SPEMD+BGG (dashed blue) models with ±1σ uncertainties (differently hatched shaded areas), as in Fig. 7, compared with the cumulative projected stellar mass profile with ±1σ uncertainties, represented by the dashed light blue curve and shaded area, respectively. Bottom: Cumulative projected stellar-over-total mass fraction profiles, color-coded for each model according to the top panel. In both panels, the vertical lines close to the x axis locate the distances from the lens galaxy center of the different multiple images, color-coded following Fig. 1. The black arrow at R = 9.32 kpc shows the effective radius of the BGG.

In the text

Fig. 9. Probability density distributions of the cosmological parameters Ωm and w in the PIEMD_Ωmw (black) and SPEMD_Ωmw (blue) models. We show the marginalized 1D histograms of each parameter and their joint 2D probability distribution. The vertical dashed lines in the 1D histograms represent the 16th, 50th, and 84th percentiles, while the solid lines in the 2D distributions represent the 0.68, 0.95 and 0.99 contour levels.

In the text

	Fig. 10. Marginalized 1D probability density distributions of the cosmological parameters Ωm and w in the PIEMD_Ωm (black), SPEMD_Ωm (blue), and MP-PIEMD_Ωm (red) models (top row), and PIEMD_Ωmw (black), SPEMD_Ωmw (blue), and MP-PIEMD_Ωmw (red) models (bottom row). The vertical dashed lines represent the 16th, 50th, and 84th percentiles for each distribution.
In the text