© The Authors 2026

1. Introduction

The circumgalactic medium (CGM), that is, the diffuse, multiphase gas surrounding galaxies, is a critical component of galaxy evolution that serves as the fuel for ongoing star formation and as the receptacle for material expelled by galactic winds (Tumlinson et al. 2017; Faucher-Giguère & Oh 2023). In particular, cool (T ∼ 10⁴ K) gas traced by low-ionization transitions such as Mg II offers unique insights into the baryon cycle at intermediate redshifts, since it is sensitive to inflows of recycled material and outflows of metal-enriched winds. Therefore, Mg II spectroscopy adds original insights that help us to understand better how galaxies evolve over cosmic time.

Over the past decade, two complementary observational strategies have emerged for absorption studies. First, large quasar–galaxy (or galaxy–galaxy) pair surveys based on observations obtained with the Multi Unit Spectroscopic Explorer (MUSE), such as the MUSE Analysis of Gas around Galaxies (MAGG) survey (Dutta et al. 2020), and the MUSE Gas Flow and Wind (MEGAFLOW) survey (Schroetter et al. 2016), as well as catalog-based studies such as the Mg II Absorber–Galaxy Catalog (MAGIICAT) (Nielsen et al. 2013b), and other efforts (e.g., Chen et al. 2010; Bordoloi et al. 2011; Kacprzak et al. 2012) have exploited bright background quasars to probe Mg II absorption out to hundreds of kiloparsecs around z ∼ 0.5 − 1.5 galaxies. These studies have mapped the radial and azimuthal dependence of cool gas around ∼L^* galaxies, revealing anisotropies linked to galactic outflows and accretion flows (e.g., Zabl et al. 2019, Wendt et al. 2021). Second, “down-the-barrel” studies of galaxies themselves using (stacked) rest-frame UV spectra of star-forming galaxies to directly measure gas kinematics against the galaxy continuum detected blueshifted Mg II absorption from outflows or redshifted signatures of inflows, thereby reducing the line-of-sight ambiguities inherent in transverse quasar probes. These analyses have revealed strong trends of an increasing absorption strength with stellar mass and star formation rate and have constrained outflow velocities and covering fractions (e.g., Weiner et al. 2009; Rubin et al. 2012; Martin et al. 2012; Rubin et al. 2018). However, these large and wide-area surveys necessarily reach only moderate depths, with continuum detection limits of m_AB ∼ 24–25 that correspond to characteristic stellar masses M_⋆ ≳ 10^9.5 M_⊙ at redshift z ∼ 1.

To directly explore the regime of lower-mass systems, where feedback might expel gas into the halo more efficiently (Dekel & Silk 1986), deeper observations are required. In this context, integral-field units (IFUs) have expanded the accessible parameter space by reaching down to low-mass systems and providing a complementary, spatially resolved view of emission line structures and gas dynamics beyond the reach of traditional slit spectroscopy (Chang et al. 2025). An example of the advantage of integral-field spectroscopy, which is sensitive to low-mass galaxies, is offered by the MUSE eXtremely Deep Field (MXDF). This ≳100 h ultradeep field has demonstrated the power of very long exposures in characterizing the kinematics and multiphase nature of the CGM around M_⋆ ∼ 10⁹ M_⊙ galaxies (Guo et al. 2023). Likewise, the MUSE Ultra Deep Field (MUDF), with its ∼140 h integration in a field of 1.5 × 1.5 arcmin² (approaching 100 h in the deepest region) and complementary HST and ground-based imaging, offers exquisite sensitivity to cool CGM gas down to stellar masses as low as M_⋆ ∼ 10⁸ M_⊙. In this study, we specifically focus on characterizing the cool phase of the CGM through detailed Mg II absorption down-the-barrel and transverse analysis in the MUDF to investigate how deep integral-field spectroscopy constrains the properties of galaxy halos at low masses.

The paper is organized as follows. In Sect. 2 we describe the MUDF sample and our data reduction. Section 3 presents the results of the down-the-barrel analysis, compares them with the literature, and interprets them through radiative transfer modeling. In Sect. 4 we present and discuss the results of the transverse analysis via galaxy pairs, providing a comparison with previous results and theoretical predictions. We summarize our main findings in Sect. 5. Throughout this paper, we adopt the standard ΛCDM cosmology with H₀ = 70 km s⁻¹ Mpc⁻¹, Ω_m = 0.3, and Ω_λ = 0.7.

2. Observations and sample selection

2.1. The MUSE Ultra Deep Field dataset

Our study relies on extremely deep observations conducted with the MUSE instrument at the Very Large Telescope (VLT) (Bacon et al. 2010) in the 1.5 × 1.5 arcmin² MUDF (ESO PID 1100.A-0528, PI: M. Fumagalli) with an integration time of ∼140 hours in wide-field mode with extended wavelength coverage (4600–9350 Å) and spectral resolution R ∼ 2000 − 4000. Thanks to the GALACSI adaptive optics system, the final image quality reaches a width at half-maximum of ∼0.73 arcsec for point sources. The full description of the MUSE observations and the data reduction strategy can be found in Lusso et al. (2019) and Fossati et al. (2019). Hubble Space Telescope (HST) follow-up observations with 90 orbits of Wide Field Camera 3 (WFC3) F140W imaging and WFC3 G141 infrared grism spectroscopy (PID 15637, PIs: M. Rafelski & M. Fumagalli) are also available, providing a deep galaxy catalog and a characterization of the spectral energy distributions (SEDs) from ultraviolet (UV) to near-infrared (NIR) wavelengths. The analysis of these HST observations is detailed in Revalski et al. (2023).

We relied on the catalog of continuum-detected sources provided by Revalski et al. (2023), which contains galaxies with redshifts constrained with HST and MUSE spectroscopy. This catalog is 50% complete at 27.6 mag (F140W) and covers the entire area observed with MUSE, as described by Revalski et al. (2023). A quality flag represents the confidence in the redshift measurement based on the presence of high signal-to-noise (S/N) lines. We only relied on high-quality redshifts (redshift quality flag of at least 3 or above), which is assigned to 25% of the sources in the catalog. The physical properties of the galaxies, such as stellar mass (M_⋆) and star formation rate (SFR), were obtained by fitting the multiwavelength photometry and the MUSE spectra simultaneously with the Monte Carlo Spectro-Photometric Fitter code (see Fossati et al. 2019; Revalski et al. 2023, 2024). For each galaxy, we also extracted a spectrum from the MUSE cube using the HST segmentation map convolved with the MUSE point spread function and registered the wavelength to the vacuum frame.

2.2. Sample selection

Our primary galaxy sample was drawn from the MUDF catalog described above and consisted of 213 sources, covering the Mg II doublet within the MUSE wavelength range, with redshifts in the range 0.6 <  z <  2.4. From this sample, we derived two distinct subsamples: one sample on which we performed a down-the-barrel analysis (i.e., using the galaxy continuum to probe foreground gas in absorption), and the other sample of projected galaxy pairs (i.e., where a background galaxy was used to investigate the halo gas of the foreground galaxy in the transverse direction).

For the first set, we selected spectra with a continuum S/N ≥ 5 per pixel, measured in a wavelength window of ∼25 Å on either side of the expected Mg II doublet (see Fig. 1). We also removed sources in which the Mg II spectral window was strongly affected by skyline residuals. This subsample comprised 84 galaxies. The physical properties (mass and SFR as a function of redshift) are summarized in Fig. 2. To ensure that the subsample was representative of the parent sample, we compared the stellar mass and redshift distributions of the two samples. Fig. 3 shows that the down-the-barrel subsample (cyan points) follows the mass-redshift distribution of the parent sample (blue points), although it is biased toward higher stellar masses due to the selection criterion on the S/N. For the second subsample of projected galaxy pairs we used to analyze the cool gas in the transverse direction, we selected galaxies following these main criteria:

1. The redshift offset between the foreground (f/g) and background (b/g) galaxies satisfies $c(z_{b/g}-z_{f/g})/(1+z_{\mathrm{pair}})>1000{\mathrm{km}\mathrm{s}}^{-1}$, with z_pair = (z_b/g + z_f/g)/2, to ensure that the background galaxy was not associated or interacting with the foreground one. This choice provided a conservative separation well above typical peculiar velocities in bound systems, excludes galaxies in the same structure, and minimizes contamination from correlated environments in absorption studies.

2. We excluded pairs in which strong lines of the b/g galaxies, such as Fe II, Al II, and Al III, might contaminate the Mg II absorption lines of the f/g galaxies.

3. We only selected b/g spectra with a continuum S/N ≥ 5 at wavelengths where we expected Mg II absorption in the rest frame of the f/g galaxy. The S/N distribution is shown in the right panel of Fig. 1.

Fig. 1. Distribution of the continuum S/N for the MUSE spectra of the 84 galaxies included in the down-the-barrel subsample (left panel) and of the 360 galaxies used for the transverse analysis, divided into bins of impact parameter: b ≤ 40 kpc (P1), 40 kpc <  b ≤ 65 kpc (P2), 65 kpc <  b ≤ 90 kpc (P3), and 90 kpc <  b ≤ 110 kpc (P4; right panel).

Fig. 2. Distributions of redshift (left), stellar mass (center), and SFR (right) of galaxies selected for the down-the-barrel analysis. As described in the text, these sources all have S/N > 5 in the spectral region surrounding Mg II. The dashed lines represent the median values of each distribution.

Fig. 3. M⋆ − z diagram of the parent sample (blue dots) and the selected subsample for the down-the-barrel analysis (cyan dots). The horizontal dashed black line at log (M⋆/M⊙)=8.5 represents the limit we adopted for building a mass-limited sample (see Sect. 3.1).

With this selection, we identified 360 pairs, which we further split into four bins of impact parameter b, defined to optimize the size and S/N of measurements for each subsample: P1, with b ≤ 40 kpc, P2, with 40 kpc <  b ≤ 65 kpc, P3, with 65 kpc <  b ≤ 90 kpc, and P4 with 90 kpc <  b ≤ 110 kpc. The final pair samples included 54, 90, 108, and 108 pairs for P1, P2, P3, and P4, respectively. In Fig. 4 we show the distributions of the f/g galaxy properties for each subsample. The mass distribution (middle panel) of the closest-pair sample (P1, yellow) appears to be slightly deviant with respect to the other bins. As galaxy pairs are selected only based on chance alignment in the plane of the sky and the foreground systems share a similar redshift distribution, we do not anticipate physical reasons for a different mass distribution. We instead hypothesize that this is an effect of the small number of galaxies in this bin. To explicitly assess whether the apparent deviation is driven by limited statistics and not an intrinsic difference in the galaxy properties, we performed two-sample Kolmogorov-Smirnov (KS) tests between P1 and the other bins and between P1 and the combined sample (P2 + P3 + P4). In each case, the resulting p value indicated that the distributions were statistically consistent, with at most marginal differences that were not significant at conventional thresholds. In addition, we carried out a Monte Carlo exercise to further test the role of small-number statistics. We generated 1000 realizations by drawing samples of equal size to P1 from a Gaussian distribution centered at log(M_⋆/M_⊙)=8.5, with the dispersion matching that of the observed sample. For each realization, we performed a KS test against P1. The resulting p-value distribution shows that 95% of realizations are statistically consistent with P1 (p value > 0.05), supporting the interpretation that the observed deviations are compatible with sampling fluctuations and do not reflect intrinsic differences.

Fig. 4. Distributions of redshift (left), stellar mass (center), and SFR (right) of the f/g galaxy sample.

2.3. Continuum normalization

To measure absorption lines, we normalized each galaxy spectrum to its continuum level, which was determined by a linear fit to the data in the spectral region on either side of the MgII doublet. In particular, we fit spectra in the windows (2747−2767) × (1+z)Å and (2830−2850) × (1+z) Å, where z corresponds to the redshift of the galaxy itself for the down-the-barrel analysis, while it is the redshift of the f/g galaxies for the measurements on the pairs. For the pair sample, we confirmed whether these regions in each b/g galaxy spectrum coincided with the gap caused by the Na laser for AO (5760−6010 Å) and/or included strong emission or absorption features of the b/g sources. When this was the case, we made minor adjustments to the continuum window to avoid incorporating extraneous features. As the spectra of galaxy pairs were analyzed through stacking, we also interpolated the spectra in the region between −2000 and 2000 km s⁻¹ on a common velocity scale with a bin size of 40 km s⁻¹, given an original rest velocity resolution between 40 km s⁻¹ and 80 km s⁻¹. The zero velocity was assigned to the λ2796 line of the Mg II doublet.

3. The cool gas along the line of sight to MUDF galaxies

3.1. Observed Mg II stacking spectra

First, we studied the relation between the Mg II lines and the physical properties of galaxies through a down-the-barrel analysis. To do this, we split the 84 galaxies into two bins of stellar mass and two bins of SFR. To preserve a comparable number of galaxies in each bin, the thresholds were defined as the median properties of the sample at log(M_⋆/M_⊙)∼9.15 and log(SFR/M_⊙ yr⁻¹)∼0.65.

We built stacked spectra, obtained by computing the median at each wavelength, focusing on the regions around the Mg II doublet for the entire sample of 84 galaxies and for each mass and SFR subsample. The resulting spectra are shown in Figs. 5 and 6 as a function of the rest-frame velocity scale, where the zero reference is on the bluer line of the doublet. The errors were estimated as the standard deviation of the flux distributions derived by performing a bootstrap analysis with 1000 resampled spectra. In the entire sample (Fig. 5), asymmetric blueshifted absorption in the K (λ2796) and H (λ2803) lines is readily visible. The K line peaks at a velocity of ∼ − 350 km s⁻¹ and shows a ∼10% absorption dip from the continuum in the center of the line. The H line reaches a velocity of ∼−100 km s⁻¹ and is about twice as strong. An emission signature at ∼ + 100 km s⁻¹ is also evident between the absorption lines.

Fig. 5. Stack of all Mg II down-the-barrel spectra (gray data points) and best-fit simulated spectra from RT-scat, assuming the gas as a smooth medium (blue line) and as a clumpy medium (red line).

Fig. 6. Mg II down-the-barrel stacked spectra (gray data points) and best-fit simulated spectra from RT-scat for the different mass (upper panels) and SFR (lower panels) bins, assuming the gas as a smooth medium (blue lines) and as a clumpy medium (red lines).

The high M_⋆ and high SFR spectra in the right panels of Fig. 6 show a deeper absorption of 30%, and 40% at the core of the K and H lines, respectively. The line core is also more aligned with zero velocity, implying that the gas giving rise to the absorption is closer to the systemic velocity. Moreover, the emission signature disappears almost entirely in the high-mass case, but is significantly reduced (but not absent) in the high SFR subsample. In the low M_⋆ and SFR cases (left panels), the emission feature instead becomes more intense, and a second weaker emission line appears at ∼ + 100 km s⁻¹ relative to the H line, giving the profiles the characteristic P-Cygni shape. To ensure that the trends observed in the stacked spectra were driven by stellar mass and not redshift, we applied a mass cut at log(M_⋆/M_⊙) = 8. We then divided the subsample into low- and high-mass galaxies using the median stellar mass of this restricted sample as a threshold. The resulting stacked spectra show consistent trends with those obtained from the original division, confirming that our results are not significantly affected by a redshift dependence on the stellar mass.

Similar features have also been observed in previous studies. For example, Weiner et al. (2009) analyzed a stacked sample of star-forming galaxies at z ∼ 1.4 from the Deep Extragalactic Evolutionary Probe 2 (DEEP2) survey, Martin et al. (2012) used Keck Low Resolution Imaging Spectrometer (LRIS) average spectra of galaxies at 0.4 <  z <  1.4, and Feltre et al. (2018) worked on stacked spectra of galaxies between 0.70 <  z <  2.34 drawn from the MUSE Hubble Ultra Deep Survey. These authors reported an increasing trend of the absorption strength with stellar mass and SFR and the appearance of a P-Cygni emission line in the Mg II λ2796 transition within the composite spectra of the lowest-mass galaxies in their sample. On the more theoretical side, Chang et al. (2025) modeled the observed Mg II spectra of a comprehensive dataset of star-forming galaxies at z ∼ 1 from the MAGG and MUDF programs, dividing the sample into different stellar mass bins. Consistent with our results, they found Mg II emission in low-mass galaxies (M_⋆/M_⊙ <  10⁹), while strong absorption becomes more prominent in higher-mass galaxies (M_⋆/M_⊙ >  10¹⁰).

3.2. Radiative transfer modeling

Taking advantage of available RT modeling (see also Prochaska et al. 2011), we fit our stacked spectra using the simulated spectra generated by the 3D Monte Carlo simulation RT-scat (further information on the Monte Carlo implementation can be found in Sect. 2.5 of Chang & Gronke 2024). Explicitly, this model assumes two different geometries for the Mg II halo, a smooth sphere and a clumpy sphere at T = 10⁴ K, with varying expanding velocities and a point source in the center. Chang & Gronke (2024) assumed two types of intrinsic spectral distributions near the Mg II doublet: a flat spectrum, representing the stellar continuum near 2800 Å, and Gaussian-like emission with a fixed flux ratio of Mg II K and H lines equal to 2 characterized by a width of v_emit and an intrinsic equivalent width EW_int. The radial outflow velocity is proportional to the distance from the central source, that is, v(r)=v_expr/R_H, where R_H denotes the outer radius. For the smooth model, the fit was performed considering five parameters: the Mg II column density N_MgII, the outflow velocity v_exp, the random motion of cold gas v_ran, the intrinsic Mg II emission width v_emit, and the intrinsic equivalent width of Mg II emission EW_int. For the clumpy model, it also included the clumps’ random velocity v_cl, and the covering factor f_cl, defined as the mean number of clumps per line of sight (for a detailed discussion of the model geometry and parameter space explored, we refer to Sect. 2 of Chang & Gronke 2024).

The results of this analysis are reported in Tables 1 and 2, where we summarize the best-fit parameters obtained with a minimum χ² for the stack of all the galaxy spectra and the stacks in each bin of mass and SFR. In Figs. 5 and 6, we also show the best fit superimposed on the stack spectra. Examining these results, we see that the RT-scat model is capable of capturing the main features of the observed spectra with both the medium geometries, including a varying degree of blueshift in the absorption lines as a function of stellar mass (and SFR), and the emergence of marked emission lines at lower masses. The model consistently reveals a higher column density for more massive (higher SFR) galaxies than lower masses, reflected in a more prominent absorption in the former subsample. However, we caution that the reported value is likely to be biased high by the intrinsic absorption in the interstellar medium (ISM) of the galaxies, which is included in the observed spectra, but not in the model (see, e.g., Rubin et al. 2012). We also find that low–mass (low–SFR) galaxies exhibit significantly higher v_exp. This suggests that their outflows are more capable of expelling gas beyond the virial radius. In contrast, in more massive systems, the outflowing material tends to stall within the halo, producing stronger absorption signatures. This was also argued by Chang et al. (2025).

3.3. Spectral features and galaxy properties

Furthermore, inspecting the individual spectra, we visually classified sources with evident Mg II absorption, with emission, or displaying a P-Cygni feature, obtaining 14, 18, and 8 galaxies in each class. In Fig. 7 we color-code the galaxies of each group in the stellar mass versus SFR plane. Consistent with the stacking analysis, we note that the median values of these properties neatly separate the class of Mg II absorbers (with higher M_⋆ and SFR) from the Mg II emitters (characterized by lower M_⋆ or SFR). P-Cygni profiles are more scattered throughout the stellar mass/SFR plane, but are concentrated in the middle of the distribution. Various authors have conducted similar exercises in previous works on different datasets, but reported comparable results: Kornei et al. (2013) in star-forming galaxies at z ∼ 1 from the DEEP2 Survey, and Erb et al. (2012) on a dataset obtained through a survey conducted with the LRIS-B spectrograph on the Keck I Telescope.

Fig. 7. Distribution of galaxies in the stellar mass and SFR plane classified as detected Mg II absorption (blue), emission (red), or P-Cygni profile (yellow). The dashed black lines represent the median of the physical sample properties.

In Fig. 8 we further built stacked spectra for the three absorption, emission, and P-Cygni profile classes. When stacking, we resorted to two ways to determine the zero velocity of each spectrum. The first way relied on the systemic galaxy redshift as inferred from other emission lines, for example, [OII] (black lines). The second method used the Mg II K line, taking the minimum for the absorption features and the maximum for the emission and P-Cygni features (filled gray lines). Within each sample, we performed a bootstrap resampling to estimate the errors of the median fluxes and to compute the distribution of the equivalent widths of the Mg II doublet, W_MgII, integrated over the velocity range between ∼ − 750 and ∼ + 1250 km s⁻¹ relative to the K line. We note that partial covering of the background continuum, commonly observed in down-the-barrel Mg II absorption studies (e.g., Rupke et al. 2005; Martin & Bouche 2009), primarily affects the inference of the column densities because it alters the depth and shape of the absorption troughs. However, the integrated equivalent width is purely an observable quantity that is also well defined in the presence of nonunity covering fractions. Since our analysis exclusively relied on W_MgII as a direct measurement, the quantities we derived do not depend on any assumptions regarding the covering fraction of the absorbing gas. The implications and possible systematic effects associated with partial covering when comparing observations with models are discussed in more detail in Sect. 4.3.

Fig. 8. Stack of spectra with detected Mg II absorption (top panel), emission (middle panel), and P-Cygni (bottom panel) profiles. The black lines represent the stack of all the galaxies of each class relative to the systemic redshift inferred from other ISM lines. The filled gray spectra instead refer to stacks on the minimum of the K line for the absorption sample and the maximum for the emission for the emitters and the P-Cygni samples.

We report the measured W_MgII and error, defined as the mean and standard deviation of these distributions. For the absorption case (top panel), we note substantial agreement between the two stacking methods, reinforcing the finding discussed above: this class is mainly characterized by static or slowly outflowing gas, with no strong infalling/outflowing components. In contrast, for the emitter case, a shift of ∼100 km s⁻¹ is visible, which is indicative of a noticeable effect induced by outflows. As expected, velocity shifts contribute to a dilution effect on the line strength. However, this has little effect on the equivalent width, which remains nearly unchanged for the emission and P-Cygni features when compared to the stacked spectra relative to the ISM systemic redshift.

4. The cool gas in the halos of MUDF galaxies via galaxy pairs

4.1. Absorption profile in galaxy halos

After exploring the cool-gas properties of MUDF galaxies along their line of sight, we studied the cool halo gas using b/g galaxies paired with f/g ones. The depth of the MUDF and the high density of the detected galaxies that is possible with MUSE enabled this study in a small area as well. Starting from the four bins of the impact parameter defined above, we built composite spectra of b/g galaxies, shown in Fig. 9. Within each sample, we performed a bootstrap resampling to estimate the errors of the median fluxes and to compute the distribution of the equivalent widths of the K line W₂₇₉₆ integrated in the velocity range of ±200 km s⁻¹. We report the mean impact parameter ⟨b⟩, and the measured W₂₇₉₆ and σ_W, defined as the mean and standard deviation of these distributions for P1 and P2. Admittedly, the detection of P2 is only marginal, but replacing the detection with a corresponding upper limit does not alter our conclusions. For P3 and P4 stacks, we instead report the 2σ upper limit, since the W₂₇₉₆ values are consistent with zero within their uncertainties. Similarly to the down-the-barrel analysis (Sect. 3.1), an imposed mass cut at log(M_⋆/M_⊙)=8 in the transverse sample yields stacked measurements consistent within 1σ with the original values. In line with previous studies (e.g., Bordoloi et al. 2011; Nielsen et al. 2013a; Rubin et al. 2018; Dutta et al. 2020), we observe a decreasing EW with radius at least up to ∼80 kpc, at which point nondetections modulate the signal (See Fig. 10).

Fig. 9. Comparison of the absorption lines in composite spectra of background galaxies in pairs from P1, P2, P3, and P4 samples, from top to bottom. The standard deviation and the confidence interval for each wavelength step were estimated by making a bootstrap within each sample.

Fig. 10. Relation of EW (W2796) vs. impact parameter (b). Measurements and upper limits (2σ) from individual pairs of f/g and b/g galaxies are plotted as grey dots and triangles, respectively. The solid purple line shows the best-fitting log-linear relation. The shaded regions indicate the 1σ and 2σ bands. The blue dots and triangles represent the detections and upper limits resulting from the stack in impact parameter bins. Various best-fit models from the literature are plotted (Nielsen et al. 2013a; Rubin et al. 2018; Dutta et al. 2020; Cherrey et al. 2025). The dashed lines represent variations in the radiative transfer model prediction for W2796, normalized with the down-the-barrel observations and extrapolated to the halo.

4.2. Individual detections and EW distribution

We further complemented the study of the mean EW from the stacking analysis with the analysis of individual measurements in the spectra of b/g galaxies, as shown in Fig. 10. The 99 individual detections scatter in the interval W₂₇₉₆ ∼ 0.3 − 3 Å and the relative fraction of the upper limits (2σ) increases with increasing impact parameters. However, significant detections ≳1 Å persist at these large distances, an effect that has been attributed to the increasing cross section around galaxies in groups (Dutta et al. 2020, 2021). About 76% of the f/g galaxies in our sample are in groups with ≥3 members, and the absorption appears to be slightly stronger in these denser environments, with a median logW₂₇₉₆ ∼ 0.09 Å at b >  90 kpc, compared to ∼ − 0.01 Å for isolated galaxies.

To characterize the resulting dependence of the EW on the impact parameter, we modeled our dataset assuming a log-linear dependence of W₂₇₉₆ on b as in Nielsen et al. (2013a), for instance,

$$\begin{array}{c}\hfill log\frac{W_{2796}}{\mathrm{\AA }}=q+m\times \frac{b}{\mathrm{kpc}}·\end{array}$$(1)

Following the method of Chen et al. (2010), the likelihood function for this model can be written as

$$\begin{array}{cc}\hfill \mathcal{L}\left(\overline{W}\right)& =({\displaystyle \prod _{i=1}^n}\frac{1}{\sqrt{2\pi s_i^2}}exp\{-\frac{1}{2}{\left[\frac{(W_i-\overline{W})}{s_i}\right]}^2\})\hfill \\ \hfill & \times ({\displaystyle \prod _{i=1}^m{\int }_{-\infty }^{W_i}}\frac{dW\prime }{\sqrt{2\pi s_i^2}}exp\{-\frac{1}{2}{\left[\frac{(W\prime -\overline{W})}{s_i}\right]}^2\}),\hfill \end{array}$$

where W_i represents the logW₂₇₉₆ value for each measurement i, and $\overline{W}$ equals the value of logW₂₇₉₆ given by the model at each b. The first product includes all detected systems, and the second product includes the systems for which our constraint on logW₂₇₉₆ is an upper limit. We further assumed that the variance of the data around the mean relation is a Gaussian with two components

$$\begin{array}{c}\hfill s_i^2={\sigma }_i^2+{\sigma }_C^2,\end{array}$$

where σ_i is the measurement uncertainty in W_2796, i, and σ_C is an additional factor that accounts for intrinsic scatter in the relation.

We sampled the posterior probability density function (PPDF) for this model using the Markov chain Monte Carlo (MCMC) method implemented in emcee (Foreman-Mackey et al. 2013). We adopted uniform probability densities for the intervals −5.0 <  m <  5.0, −10.0 <  q <  10.0, and 0 <  σ_C <  10.0 as priors. Markov chains generated by 100 walkers taking 6000 steps each (discarding the first 1000 steps) provide an exhaustive sampling of the PPDF in each parameter dimension (see Rubin et al. 2018). We defined the best value of each parameter and its uncertainty interval as the maximum and the 16th and 84th % iles of the inferred PPDFs.

The resulting relation is flat, with a slope of m = −0.003 ± 0.002, an intercept of q = −0.3 ± 0.1, and an intrinsic scatter of σ_C = 0.5 ± 0.1. In Fig. 10 we show the fit result (in purple) with the 1σ and 2σ bands compared with the fit previously performed on other datasets from quasar (dotted lines, Nielsen et al. 2013a; Dutta et al. 2020; Cherrey et al. 2025) or galaxy b/g (solid line, Rubin et al. 2018) sightlines. The inferred dependence of the EW on the impact parameter appears to be consistent with the result of the stacking analysis discussed above. Our trend is consistent with the trend reported by Rubin et al. (2018), whose results were also based on stacking background galaxy spectra, but are flatter than the others. As noted by Rubin et al. (2018), differences in slope among the datasets might arise from selection effects, such as the inclusion of isolated or group-member galaxies. This effect was tested by Dutta et al. (2020), who found a flattening relation when group galaxies were included. However, this difference seems marginal in their analysis. As noted by Rubin et al. (2018), other factors such as the interval of b considered for the fit and the different nature of quasar versus galaxy spectroscopy also play a role in modulating the trend. In particular, background quasar sightlines probe narrow pencil beams (of a few parsecs, smaller than typical Mg II cool clouds), while background galaxies sample regions that are several kiloparsecs across, averaging over multiple clouds. As a result, the covering fraction measured and the corresponding transverse profile from galaxy–galaxy and galaxy-quasar pairs should be compared with caution.

4.3. Comparison with model predictions

Next, we compared the prediction of the halo model described above, constrained down the barrel, with the observations of the EW distribution measured in the transverse direction. To do this, we reconstructed the Mg II column density and outflow velocity with the “all” parameters from Table 1 (smooth model) and from Table 2 (clumpy model), setting the remaining parameters to the values fixed by Chang et al. (2025). R_H was set to 120 kpc to include the maximum impact parameter we observed. The predicted spectra were constructed as follows.

For the smooth model, we computed the intersection of the line of sight with the spherical volume at each selected impact parameter. This path was then divided into segments of equal length. For each segment, we calculated the line-of-sight velocity and assigned an optical depth along with a Gaussian velocity profile centered on that velocity. The width of each Gaussian was set by the thermal broadening of Mg II at 10⁴ K. The complete absorption spectrum was obtained by summing the contributions from all segments along the line of sight.

For the clumpy model, the approach was similar, with a few additions. We estimated the column density of each clump as $N_{\mathrm{MgII},\mathrm{cl}}=\frac{3N_{\mathrm{MgII}}}{4f_c}$ and determined the number of clumps intersected at each impact parameter as $2f_c\sqrt{1-{\left(\frac{b}{R_H}\right)}^2}$. For each clump, we calculated the line-of-sight velocity. Different from the smooth case, we added a random motion term to the radial outflow velocity, modeled as a Gaussian deviate 𝒩(0,1) scaled by the clump random motion ν_cl obtained from the fit. Then, we assigned an optical depth and a Gaussian velocity profile to each clump, following the same method as in the smooth case.

The dashed lines in Fig. 10 represent the W₂₇₉₆ versus b relation derived from these models. In general, the halo model predicts an almost flat profile, consistent with the fit result, followed by a drop at the assumed outer edge of the halo. However, this cutoff at the external radius of the halo is an inherent feature of the model. In contrast, the CGM is not expected to have sharp edges, and the contribution of additional halos becomes progressively more relevant with increasing radius.

The smooth model lies in the upper envelope of the detected systems, but exceeds the best fit of the observed EW. The clumpy model, which describes the expected Mg II distribution in the halo better (e.g., Schroetter et al. 2019), predicts a lower value of EW, but lies consistently above the fit. This excess likely arises from simplifying assumptions in the model.

First, the halo models omit an ISM component, which is present in the data. Second, they assume either a uniform gas density or a constant number of clumps, whereas in reality, the clumpiness and the overall density of cold gas both decline with galactocentric radius. Beyond the missing ISM contribution and radial density gradient, the biconical nature of the outflows differs from the spherical one assumed by the model and also depends on the inclination of the line of sight with the galactic plane. Previous observations have indeed revealed an anisotropic distribution of circumgalactic cold gas (e.g., Schroetter et al. 2019; Bouché et al. 2013; Kacprzak et al. 2012; Bouché et al. 2012; Ho & Martin 2019), which likely amplifies the discrepancy between spherical-model predictions and sightline measurements. As a result of these combined effects, we expect a broad distribution in column density with a lower mean than measured in the down-the-barrel line of sight, which more readily crosses the outflow component. Similar geometric considerations are also expected to affect the effective velocity along the line of sight.

Following this idea, using the same parameters as obtained with the clumpy spherical model, we considered a bipolar wind geometry with an opening angle of 30 deg (Chang et al. 2025; Guo et al. 2023). By generating random lines of sight for each impact parameter, we found a covering fraction of ∼14%, which we used as a scaling factor for the column density. The result, shown as the dashed blue line in Fig. 10, aligns better with the observations.

A further caveat concerns the role of the covering fraction in shaping the predicted absorption strength, as introduced in Sect. 3.3. Since our radiative-transfer framework assumes a homogeneous statistical distribution of clumps for the clumpy model by construction, the effective area-covering fraction is encoded in the mean number of intercepts per sightline and is therefore fixed at all radii. The question arises whether the derived equivalent width distributions depend on the form of the background beams (either extended, as for galaxies, or compact, like in quasars). To assess this effect, we carried out a simple numerical experiment in which, instead of using single skewers, we simulated beam-like sightlines sampling a finite transverse area. We then compared the resulting equivalent widths with those obtained from infinitesimal skewers.

For the smooth model, the two approaches yield identical results. This is expected, as averaging over multiple rays in a medium with uniform density is equivalent to evaluating a single line of sight through that medium. For the clumpy model, where the number of clumps per sightline is drawn from the mean parameter f_c, beams produce EW values that are only marginally lower than those from skewers. This small difference arises from the averaging over the finite beam area, but remains negligible because the clump population is defined through a mean statistical covering that already captures the effective probability of intercepting gas. In this modeling framework, the beam-covering fraction is therefore inherently contained in the definition of f_c, and its impact on the mean predicted absorption is minimal.

This test confirms that within the assumptions of the present RT setup, adopting beams rather than skewers does not significantly modify the predicted equivalent-width profiles. However, it also highlights that the model has simplistic characteristics, that is, it does not capture any radial variation in the covering fraction, which is expected observationally. A more sophisticated treatment of the spatial distribution of cool gas in the models will be the subject of future work.

5. Summary and conclusions

Our study leveraged deep MUSE observations from the MUDF survey to investigate the cool-gas properties in and around galaxies at 0.5 ≲ z ≲ 2. Using a catalog of continuum-detected sources with high-quality redshifts from Revalski et al. (2023), we constructed two samples of galaxies with masses between 10⁷ and 10^11.5 M_⊙ and an SFR between 10⁻⁷ and 10³ M_⊙ yr⁻¹: one sample for a down-the-barrel analysis to probe foreground gas in absorption through the galaxy continuum, and the other sample for projected galaxy pairs to investigate halo gas in the transverse direction.

For the down-the-barrel analysis, we selected 84 galaxies with S/N ≥ 5 and minimum skyline contamination to examine their Mg II absorption features. In the stacked spectra of two mass and SFR bins, we detected asymmetric blueshifted absorption features Mg II indicative of outflows. The absorption strength increases with stellar mass and SFR, whereas P-Cygni profiles are more prevalent in lower-mass, lower-SFR galaxies. We modeled the spectra with the 3D Monte Carlo radiative transfer simulation RT-scat (Chang & Gronke 2024), finding that high-mass galaxies exhibit higher column densities and lower outflow velocities. In contrast, low-mass galaxies show weaker absorption, but faster outflows, suggesting more efficient gas expulsion in low-mass systems (Chang et al. 2025).

For the projected pair sample, we analyzed 360 galaxy pairs, with b/g spectra S/N ≥ 5 grouped by impact parameter and avoiding the contamination of b/g strong lines, to study the transverse distribution of Mg II absorption in the CGM. The Mg II absorption strength exhibits a declining radial profile around galaxies. The dependence of W₂₇₉₆ (including 2σ upper limits for nondetections) on the impact parameter is flatter than reported in previous studies using quasar spectroscopy (e.g., Nielsen et al. 2013a; Dutta et al. 2021), but agrees with the analysis of galaxy spectroscopy by Rubin et al. (2018). Various effects related to the environment of galaxies, the interval of the impact parameter considered for the fit, and the different nature of quasar versus galaxy spectroscopy contribute to this discrepancy.

Next, we compared the Chang & Gronke (2024) halo model predictions for a smooth and a clumpy medium, constrained through a down-the-barrel measurement, with observations of the EW distribution in the transverse direction. Using the best estimates of the column density and outflow velocity for the entire sample, we reconstructed the halo profile. The model predicts an almost flat EW in line with the data, but the predicted EW exceeds the observed values, especially for the smooth model. While the clumpy model agrees better, it still overshoots the observations considerably. We interpret this discrepancy as caused by the lack of an ISM component in the model and differences in the outflow geometry and galactic plane inclination. Indeed, we showed that the resulting equivalent width shifts in the direction of the observations when we included the geometric dilution induced by a biconical outflow with an opening angle of 30 deg.

Our analysis offers new insight into the distribution of cool gas in the CGM of galaxies, considering the absorption properties along the line of sight to their star-forming regions and in the transverse direction inside their halos, using a deep dataset that reaches masses as low as M_⋆ ≲ 10⁸ M_⊙. Thus, this study complements the work of Chang et al. (2025), who employed a similar dataset to combine absorption diagnostics along the line of sight with emission measurements in the halos. We can thus directly compare the findings obtained by relying on two different probes inside the halo.

Applying radiative transfer modeling to the down-the-barrel absorption/emission spectrum and the Mg II surface brightness in emission inside halos, Chang et al. (2025) constrained the key halo parameters in Tables 1 and 2. They observed that low-mass galaxies display prominent Mg II emission and high outflow velocities in the core (b <  10 kpc) and nearest the halo regions (10 <  b <  30 kpc), and they concluded that these features can be explained by fast-moving outflowing gas components. Conversely, high-mass galaxies tend to show strong core absorption, with lower outflow velocities indicative of abundant, relatively static cold gas-accompanied by more extended halo emission.

Similar conclusions are confirmed by our down-the-barrel analysis alone: low-mass galaxies are better modeled by a lower column density and a higher velocity than more massive systems. This agrees well with the picture that more massive galaxies retain more cool gas in their halos, whereas lower-mass objects can more effectively push gas outward under the effect of winds. It might be concluded that down-the-barrel information alone can constrain the system fully. However, we showed that extrapolating the information recovered along the line of sight to galaxies to their halos yields a higher EW than observed. This discrepancy can be mitigated by the effects of an asymmetric wind geometry, as often invoked in Mg II absorption studies (Schroetter et al. 2019). The transverse component in absorption is thus still critical to add information about the spatial distribution and geometry of the cool phase of the outflow (see also Guo et al. 2023, for a complementary approach in emission).

Using complementary techniques in deep datasets, we conclude that the mass is a fundamental parameter in shaping the properties of cool gas inside halos, but geometry still plays a critical role. As various diagnostics in absorption and emission become more available for the same samples, the constraining power of refined radiative transfer models will yield a comprehensive picture of the cold CGM.

Acknowledgments

We thank Kate H.R. Rubin, Rongmon Bordoloi, and Davide Gerosa for discussions that improved this work. This work is based on observations collected under ESO programme ID 1100.A-0528. Based on observations with the NASA/ESA Hubble Space Telescope obtained from the MAST Data Archive at the Space Telescope Science Institute, which is operated by the Association of Universities for Research in Astronomy, Incorporated, under NASA contract NAS5-26555. Support for program numbers 15637 and 15968 was provided through a grant from the STScI under NASA contract NAS5-26555. These observations are associated with program numbers https://archive.stsci.edu/proposalsearch.php?mission=hstid=6631, https://archive.stsci.edu/proposalsearch.php?mission=hstid=15637, and https://archive.stsci.edu/proposalsearch.php?mission=hstid=15968.

References

All Tables

All Figures

	Fig. 1. Distribution of the continuum S/N for the MUSE spectra of the 84 galaxies included in the down-the-barrel subsample (left panel) and of the 360 galaxies used for the transverse analysis, divided into bins of impact parameter: b ≤ 40 kpc (P1), 40 kpc <  b ≤ 65 kpc (P2), 65 kpc <  b ≤ 90 kpc (P3), and 90 kpc <  b ≤ 110 kpc (P4; right panel).
In the text

	Fig. 2. Distributions of redshift (left), stellar mass (center), and SFR (right) of galaxies selected for the down-the-barrel analysis. As described in the text, these sources all have S/N > 5 in the spectral region surrounding Mg II. The dashed lines represent the median values of each distribution.
In the text

	Fig. 3. M⋆ − z diagram of the parent sample (blue dots) and the selected subsample for the down-the-barrel analysis (cyan dots). The horizontal dashed black line at log (M⋆/M⊙)=8.5 represents the limit we adopted for building a mass-limited sample (see Sect. 3.1).
In the text

	Fig. 4. Distributions of redshift (left), stellar mass (center), and SFR (right) of the f/g galaxy sample.
In the text

	Fig. 5. Stack of all Mg II down-the-barrel spectra (gray data points) and best-fit simulated spectra from RT-scat, assuming the gas as a smooth medium (blue line) and as a clumpy medium (red line).
In the text

	Fig. 6. Mg II down-the-barrel stacked spectra (gray data points) and best-fit simulated spectra from RT-scat for the different mass (upper panels) and SFR (lower panels) bins, assuming the gas as a smooth medium (blue lines) and as a clumpy medium (red lines).
In the text

	Fig. 7. Distribution of galaxies in the stellar mass and SFR plane classified as detected Mg II absorption (blue), emission (red), or P-Cygni profile (yellow). The dashed black lines represent the median of the physical sample properties.
In the text

Fig. 8. Stack of spectra with detected Mg II absorption (top panel), emission (middle panel), and P-Cygni (bottom panel) profiles. The black lines represent the stack of all the galaxies of each class relative to the systemic redshift inferred from other ISM lines. The filled gray spectra instead refer to stacks on the minimum of the K line for the absorption sample and the maximum for the emission for the emitters and the P-Cygni samples.

In the text

	Fig. 9. Comparison of the absorption lines in composite spectra of background galaxies in pairs from P1, P2, P3, and P4 samples, from top to bottom. The standard deviation and the confidence interval for each wavelength step were estimated by making a bootstrap within each sample.
In the text

Fig. 10. Relation of EW (W2796) vs. impact parameter (b). Measurements and upper limits (2σ) from individual pairs of f/g and b/g galaxies are plotted as grey dots and triangles, respectively. The solid purple line shows the best-fitting log-linear relation. The shaded regions indicate the 1σ and 2σ bands. The blue dots and triangles represent the detections and upper limits resulting from the stack in impact parameter bins. Various best-fit models from the literature are plotted (Nielsen et al. 2013a; Rubin et al. 2018; Dutta et al. 2020; Cherrey et al. 2025). The dashed lines represent variations in the radiative transfer model prediction for W2796, normalized with the down-the-barrel observations and extrapolated to the halo.

In the text