© The Authors 2023

1. Introduction

It is now clear that the evolution of galaxies must be linked to their interaction with their surrounding gaseous halos, the circumgalactic medium (CGM). Of particular importance is the CGM cool phase (T ∼ 10⁴ K, see Tumlinson et al. 2017), which traces the flows of material in and out of galaxies, known as the baryon cycle, and represents a huge reservoir of fuel for star formation. However, granted that this medium holds the key to understanding how galaxies grow and evolve, many of its properties remain unclear.

Because of its low emissivity, observing the cool CGM directly in emission is challenging and currently most of these observations are focused on the high-redshift Universe (z ≳ 2), where there seems to be an ubiquity of extended Lyα nebulae illuminated by the central quasars (QSOs) or galaxies (e.g., Cantalupo et al. 2014; Farina et al. 2019). At z ≲ 1, there are instead only a handful of cases where the cool CGM has been observed in emission, either around starbursts and galaxy groups (e.g., Chen et al. 2019; Burchett et al. 2021; Leclercq et al. 2022) or thanks to stacking techniques (e.g., Zhang et al. 2018). The vast majority of the cool CGM detections therefore come from absorption studies, which have been extensively performed for decades both at high and low redshift using absorption lines as tracers, namely those of hydrogen and low-ionization ion species such as Mg II, C II, and Si II in the spectra of background quasars (e.g., Bergeron 1986; Churchill et al. 2000; Zhu & Ménard 2013; Keeney et al. 2017; Zahedy et al. 2019; Qu et al. 2023).

A major issue in understanding the elusive properties of this material is the presence of typically only one single line of sight per galaxy, meaning that the properties of the cool CGM are usually studied in a statistical way, combining detections from surveys of tens or hundreds of galaxy–QSO pairs (e.g., Lanzetta et al. 1995; Chen et al. 2010a; Werk et al. 2012; Wilde et al. 2021). This has proven very useful in unveiling some of the properties of this medium, such as its absorption strength or its covering fraction as a function of the distance from the central galaxy (e.g., Nielsen et al. 2013; Huang et al. 2021). However, due to the single-pencil-beam limitation, there are still very strong uncertainties on the general dynamics and structure of this gas. High-spectral-resolution observations have shown that this medium is likely composed of multiple clumps of gas, or clouds, which produce multiple kinematic components along the same line of sight and are bound to the galaxy halo (e.g., Borthakur et al. 2015). However, whether these clouds are primarily inflowing toward (e.g., Bouché et al. 2013) or outflowing (e.g., Schroetter et al. 2019) away from the central galaxies is still a matter of debate.

From a theoretical perspective, the cool CGM properties across the galaxy halo can be studied using cosmological and zoom-in simulations. However, given the low resolution in the CGM (kiloparsec scale), such simulations are not yet converged in terms of the cool gas dynamics and structure (e.g., van de Voort et al. 2019) and the picture remains uncertain. On the other hand, high-resolution (parsec-scale) hydrodynamical simulations focused on the motion of a single cool cloud interacting with the surrounding hotter CGM phase (T ∼ 10⁶ K) have shown that the dynamics and the survival (e.g., Marinacci et al. 2010; Gronke & Oh 2020a) of such clouds highly depend on their initial size and mass, as found also by analytical models (e.g., Nipoti & Binney 2007; Afruni et al. 2019; Fielding & Bryan 2022). Constraining these sizes is therefore crucial to advancing our understanding of the cool CGM dynamics and its connection with the central galaxy. Nevertheless, estimating these values directly from observations has proven extremely challenging.

Several studies have resorted to photoionization models, such as CLOUDY (e.g., Ferland et al. 1998), to estimate the thickness of the cool CGM clouds (e.g., Werk et al. 2014), but these estimates are subject to large uncertainties (up to several orders of magnitude) due to the inherent assumptions of the models (Haehnelt et al. 1996). A more direct way to estimate the spatial extension of these clouds is to use multiply lensed quasars (e.g., Smette et al. 1995; Lopez et al. 1999, 2007; Ellison et al. 2004; Chen et al. 2014; Zahedy et al. 2016; Rubin et al. 2018a), which allow the absorbers to be probed through more than just one line of sight. However, the results on absorber sizes are inconclusive because of the small number of bright lenses and the small number (2 or 3) of lines of sight per system, preventing a broader view of the gas structure.

The present work is largely motivated by the promising approach of using extended background sources (Steidel et al. 2010; Bordoloi et al. 2014; Rubin et al. 2018b). Using a sample of 27 galaxy-galaxy pairs taken from the PRIMUS survey (Coil et al. 2011), Rubin et al. (2018b) estimated the scale over which the absorption strength of the cool CGM – traced by Mg II absorption – is expected to vary. Through their analysis, these latter authors were able to put a lower limit on this quantity of about 2 kpc. This means that they did not observe variations in the absorption strength on scales smaller than this threshold and that even if the single clouds have smaller sizes, their distribution and kinematics result in a constant absorption strength on scales of at least a few kiloparsec. This coherence scale therefore provides an idea of the clustering scale of the cool medium. In the present work, we aim to estimate this quantity with greater accuracy through the use of gravitational arcs.

In recent years, the ARCTOMO Collaboration¹ has developed a new technique to spatially resolve the cool CGM in absorption across individual galaxy halos. The technique, referred to as “arc-tomography”, uses resolved spectroscopy of giant gravitational arcs as background sources, which are observed with integral-field spectrographs such as VLT/MUSE (Bacon et al. 2010). The observations enable studies of the cool CGM distribution and kinematics (as traced by Mg II absorption) of single galaxies (e.g., Lopez et al. 2018, 2020; Tejos et al. 2021; Fernandez-Figueroa et al. 2022) in an unprecedented fashion by directly mapping its properties across the arc. In the present work, we use this unique dataset (hereafter “ARCTOMO data”) combined with the predictions from empirical models to constrain the cool CGM coherence scale, which is defined in the same way as in Rubin et al. (2018b).

In Sect. 2, we summarize our data, describing the sample selection, data reduction, and absorption line search, and define the main observational constraints that we use in this work. In Sect. 3, we outline the creation of our fiducial CGM model and the technique used to compare the model outputs with the data. In Sect. 4, we report our findings regarding the cool CGM coherence scale, in Sect. 5 we discuss our interpretation of the coherence scale, the limitations of our analysis, a comparison with previous studies, and the implications of our findings for theoretical models, while in Sect. 6 we summarize the main conclusions of this study. Throughout the paper, we assume a flat ΛCDM cosmology with H₀ = 70 km s⁻¹ Mpc⁻¹, Ω_Λ = 0.7 and Ω_M = 0.3.

2. Observational data

Table 1 lists the summarised properties of the three systems studied in this work, that is, the sample of spatially independent Mg II detections and their identified absorbing galaxies (hereafter referred to as “G1”). These systems were found in giant gravitational-arc fields observed by us with VLT/MUSE. The selection for this work was based on absorbing galaxies that (1) have a large number of spatially independent Mg II detections in their CGM; (2) share similar redshifts, stellar masses, luminosities, and star formation rates; and (3) are considered “isolated”. The “isolation” character only strictly holds for G1 toward PSZ1 G311.65−18.48, where no galaxies at the same redshift (Δv ≲ 1000 km s⁻¹) are found in the MUSE field (Lopez et al. 2020), while for G1 toward SGAS J1226+2152 one galaxy at the same redshift is found over 225 kpc away (delensed projected distance; Tejos et al. 2021) and for G1 toward RCS2 032727−132623 two galaxies are found at the same redshift over 200 kpc away (Lopez et al. 2018). In addition RCS2 032727−132623 G1 itself is resolved into three galaxies in the HST images (Lopez et al. 2018) and therefore integrated properties are provided in Table 1. Although these systems have been presented individually in the indicated references, the present work includes improved data reduction and a repeat analysis of the absorption line profiles, the details of which are provided in the following.

2.1. MUSE observations and data reduction

VLT/MUSE integral-field observations of these three ARCTOMO fields, PSZ1 G311.65−18.48, RCS2 032727−132623, and SGAS J1226+2152, were carried out in service mode within ESO programs 297.A-5012(A) (PI: Aghanim), 098.A-0459(A) (PI: Lopez), and 0101.A-0364(A) (PI: Lopez), respectively. All observations were carried out in wide-field mode, which provides a field of view of ≈1′×1′, a spaxel size of 0.2″, and spectra of resolving power R = 2000 − 4000 in the range of 4650−9300 Å at a dispersion of 1.25 Å. In addition, SGAS J1226+2152 was observed using the adaptive-optics extended wavelength mode.

Observations proceeded during dark time and under good weather conditions. To minimize instrumental and sky artifacts, short (640−700 s) individual exposures were taken with small spatial ditherings and cumulative 90° rotations. The number of exposures and the total integration times varied according to arc brightness, resulting in relatively homogeneous continuum signal-to-noise ratios (S/Ns) near the relevant Mg II absorption. For more details on the individual observing conditions, we refer the reader to the references shown in Table 1.

Data of all three fields were reduced in a homogeneous fashion using the MUSE pipeline (Weilbacher et al. 2012) using the EsoRex (ESO Recipe Execution Tool) tool. The wavelength solution was calibrated to vacuum directly from the pipeline. The effective PSF FWHM in the combined cubes ranges between 0.7″ and 0.8″. Residual sky contamination was removed using the Zurich Atmosphere Purge code (ZAP v2.1; Soto et al. 2016) on the final combined cubes, using the default parameters.

2.2. Automated absorption line analysis

To obtain spatially resolved rest-frame EW maps, the absorption line analysis proceeds automatically in two basic stages: (i) extraction of arc spectra and (ii) absorption line fitting. For each data cube, the specific tasks are as follows.

We first rebinned the cube spatially. This task serves two purposes: account for the cross-talk between neighboring MUSE native spaxels due to the seeing, and increase the S/N. In a sense, this is equivalent to extracting spectra of many resolved sources in a coherent fashion. In this work, we use “3 × 3 binning”, meaning that each binned cube spaxel is 0.6″ on a side and contains a weighted average spectrum out of nine native spaxels. This choice maximizes spatial sampling while maintaining a low level of seeing-induced cross-correlation between neighboring spaxels (Tejos et al. 2021). We then proceeded to select binned spaxels in a masked region that encompasses the arc based on S/N. The S/N is calculated in a featureless spectral region close to the expected Mg II absorption using the actual flux RMS. We use S/N = 2 as a threshold. This limit is arbitrary; it serves the purpose of preselecting arc spectra whilst ruling out spaxels that are contaminated by the sky, especially at the arc borders. The final Mg II selection used in this work is done at the absorption-line fitting level.

To obtain EWs, we normalized the spectra and fitted two Gaussians that account for the Mg IIλλ2796, 2803 doublet. Describing the profiles with Gaussians assumes the absorption signal is dominated by the instrumental profile. Therefore, the four free parameters are two amplitudes², one common line width, and one common redshift. Fits are considered successful if the EW measurement is larger than two times its 1σ error in both doublet lines. The scheme yields 91 detections among the three fields. The fitted parameters are reported in Appendix A. The present work is based on the (rest-frame) EW of the 2796 doublet line, hereafter referred to as EW; these are assigned to the binned spaxels to create the maps shown in Fig. 1 (both in the image and the absorber planes). Only successful fits are considered. Also, as we are using an updated selection criterion, these maps are different from those in Lopez et al. (2018, 2020) and (Fernandez-Figueroa et al. 2022).

Fig. 1. Maps of EWs of the 2976 Mg II doublet line for all the detections found in the three fields analyzed in this work: PSZ1 G311.65−18.48 (top), RCS2 032727−132623 (center), and SGAS J1226+2152 (bottom). Both the image (observed) plane and the reconstructed absorber plane (dashed-line rectangle) are shown, with the crosses showing the positions of the three identified absorbing galaxies, referred to collectively as G1. In SGAS J1226+2152, the absorption on top of G1 is shown for completeness but not used in the analysis.

2.3. Transverse distances in the absorber plane

At the absorber redshifts, all right ascension and declination (RA, Dec) coordinates are affected by lensing. We transform those to “delensed” coordinates using the lens models described in detail in the references listed in Table 1. This transformation produces a map of delensed spaxels in the absorber plane (also referred to here as the “reconstructed” plane; see Fig. 1) at the redshift of the absorber. The delensed coordinates are used to measure projected physical separations between spaxel centers and also between spaxels and the central galaxy. We call this latter quantity the “impact parameter”, or R.

Lens models are precise to around 5% (Lopez et al. 2020) and this uncertainty propagates to R. We verified that this uncertainty is negligible for the purposes of the present study and does not affect our findings. However, the models are also subject to systematic uncertainties depending on how well they are constrained observationally, and so the solutions might not be unique. For instance, our model for the SGAS J1226+2152 field considers two perturbers, while the model presented in Mortensen et al. (2021) considers just one, resulting in two different reconstructed absorber planes for this field. However, the two different reconstructed spaxel maps can be reasonably well matched with each other by a translation of coordinates only. Hence, R and transverse distances between spaxels are unlikely to be strongly affected by systematic uncertainties in the lens model (Tejos et al. 2021).

2.4. Main observational constraints

In this work, we combine EWs and spaxel positions with respect to the central galaxy of the three fields described above. The underlying assumption is that the cool CGM has self-similar properties across our sample of galaxies. This is common practice in absorption studies of the CGM, but while it is normally assumed for samples of tens or even hundreds of galaxies, here we focus instead on the halos of only three galaxies, with similar properties (see Table 1). To take into account the slight differences in stellar mass between the three galaxies, we define a rescaled version of R, as

$$\begin{array}{c}\hfill R^{\prime }=R/{(M_{\ast }/M_{\ast ,0})}^a,\end{array}$$(1)

where M_*, 0 = 10^10.3 M_⊙ and a = 0.26. This definition was introduced in Huang et al. (2021), who conducted a survey with MagE of Mg II absorbers around almost 400 galaxies at z < 0.4 using more than 150 quasars as background sources. In particular, the values of M_*, 0 and a reported above are derived from fitting the Mg II EW as a function of the impact parameter and the stellar mass (see below) for their sample of isolated star-forming galaxies, which are similar to the objects analyzed in the present study. We then investigate how the Mg II EW varies in our data as a function of R′ and as a function of the spaxel separation, as described below. We emphasize that, for simplicity, we use only confirmed Mg II detections. We consider that most nondetections are due to spectra that (despite our selection based on S/N; see Sect. 2.2) are at the arc borders and are therefore contaminated by the sky. Given that the sky contamination cannot be trivially implemented in our models (see Sect. 3), we plan to address this issue (and to include also the nondetections) in future works.

In the left panel of Fig. 2, we show the Mg II EWs of the three ARCTOMO fields in blue with their respective 1σ uncertainties, all as a function of R′: we can see that the EW decreases with increasing distance from the central galaxy. This distribution represents the first observational constraint used in this study. As a comparison, we also report the data coming from the survey of Huang et al. (2021) in orange (considering only the sample of isolated star-forming galaxies), with the corresponding best-fit curve in black, given by

$$\begin{array}{c}\hfill log\mathrm{EW}=a_0+a_{1}log{R}^{\prime },\end{array}$$(2)

Fig. 2. Main observational constraints extracted from the ARCTOMO data. Left: Mg IIλ2796 rest-frame EW as a function of R′ (Eq. (1)). Arctomographic data are shown in blue, and data from Huang et al. (2021) in orange, with 2σ upper limits denoted by downwards triangles. The black curve shows the best fit of the Huang et al. (2021) sample (Eq. (2)) and the vertical dashed line marks the region where our data tend to diverge from this curve (these data are excluded from our analysis). Right: distribution of the fractional difference of EW in the arctomo data as a function of the spaxel separation. The blue points show the values of all the spaxel pairs, while the blue curve and the gray band show the median values – with the 1σ uncertainties – of a binning of the same distribution (each bin contains 150 values). The black curves show the results for distributions obtained by choosing different orientation angles for the three MUSE fields (see Sect. 2.4 for more details).

where a₀ = 1.42, a₁ = −1.07. The ARCTOMO data lie on top of those obtained from the quasar survey and are well described by the quasar-based best-fit relation. The main difference can be seen in the region where R′≲10 kpc, where however there are no data available from the Huang et al. (2021) sample and the best-fit relation represents only an extrapolation. All the data with R′≲10 kpc are therefore excluded from our analysis (see Appendix A).

The second main observable that we use to constrain the coherence length of the cool CGM is given by the distribution of the EW fractional differences, denoted ΔEW, as a function of spaxel separation. We define the fractional EW difference between two spaxels as (e.g., Ellison et al. 2004; Rubin et al. 2018a; Okoshi et al. 2021)

$$\begin{array}{c}\hfill \mathrm{\Delta }\mathrm{EW}=\frac{|{\mathrm{EW}}_1-{\mathrm{EW}}_2|}{max({\mathrm{EW}}_1,{\mathrm{EW}}_2)}·\end{array}$$(3)

We calculate this quantity for all the pairs of spaxels with detections and combining together the three MUSE fields. The error on ΔEW is propagated from the uncertainties on the individual EWs (assuming that the errors of EW₁ and EW₂ are uncorrelated). We then create a distribution of ΔEW as a function of spaxel separation. To calculate the distances between spaxels located in different MUSE fields, we assume random orientation angles between the fields. A Monte Carlo simulation shows that these relative orientations do not affect the distribution of ΔEW versus separations. One of these realizations is shown in the right-hand panel of Fig. 2 (uncertainties are not shown for clarity). We bin separations by equal number of pairs (150) and compute the median ΔEW (blue solid line) and the 32nd and 68th percentiles (gray band) in each bin. Here, we can see that ΔEW is on average already greater than zero for contiguous spaxels and tends to increase with spaxel separation. We note that the peculiar shape of this distribution for separations of between 30 and 60 kpc is likely nonphysical, and is more likely due to the specific spaxel configuration we are using (a similar shape is indeed visible also in the synthetic observations of our models). In this work, we focus on the small spaxel separations in any case.

Finally, the three additional black curves on the right-hand side of Fig. 2 show the median values for three other distributions with different choices of field orientation: these lines are all within the gray band, showing that, as mentioned above, the orientation does not affect this distribution. This suggests that the EW distribution is isotropic, such that it does not show any increment or decrement of EW in a particular direction with respect to the host galaxy.

3. Models

3.1. Basic approach

To interpret our observables, we built empirical models of the cool CGM distribution. We created 2D maps of Mg II EW, where a crucial parameter is given by the gas coherence length, C_length. As in Rubin et al. (2018b), we define the coherence length (or scale) as the size over which the EW of Mg II remains constant. In practice, our maps are divided in cells, the size of which is equal to C_length: the EW cannot vary within a cell, while different cells can have different EW values.

The two main assumptions when creating our models are (i) the “intrinsic” distribution of EWs across the halo follows the results of the quasar survey from Huang et al. (2021), justified by the fact that our data are consistent with their findings; (ii) the EW distribution is isotropic across the plane, justified by the distribution of fractional EW differences in our data (see Sect. 2.4 and Fig. 2). We then performed synthetic observations, tracing MUSE-like spaxels that resemble the observations outlined in Sect. 2 and assuming different values for the coherence scale. The outputs of such artificial observations can be directly compared with the observational constraints outlined in Sect. 2.4, allowing us to estimate the choice of C_length that best reproduces the ARCTOMO data.

3.2. Creation of EW maps

We populate the rescaled (Eq. (1)) projected plane (x′,y′) surrounding the central galaxy (located at x′ = 0, y′ = 0) using the EW–R′ relation described by Eq. (2) and a declining covering fraction profile based on the constraints found by Huang et al. (2021)³. More specifically, we first divided the plane into a grid, with a cell size equal to the assumed coherence length, and then into concentric rings with a width equal to the coherence length itself. We then populated each ring with a number of “systems” n_ring equal to the total number of cells in the ring multiplied for the predicted covering fraction. This is done by assigning to the selected n_ring cells a value of EW taken from a Gaussian distribution centered on the mean value of EW predicted by Eq. (2) for those distances and with a width equal to the intrinsic scatter σ of the data from Huang et al. (2021). This scatter is equal to 0.278 for R′< 100 kpc, which is the region analyzed in this paper. Finally, the “empty” cells have EW = 0 Å⁴.

The outputs of the method described above are shown in Fig. 3, with their dependence on the gas coherence length. In particular, C_length is equal to 1, 4, and 10 kpc going from top to bottom. We can see that, even though the general trend is similar (by construction), the maps appear quite different from one another because of the different scale over which the EW is allowed to vary. We note how all maps show “holes”, illustrating the effect of the assumption that the cool CGM is made of clouds whose total covering fraction is lower than one. Both C_length and covering fractions have an effect on the ARCTOMO Mg II detections.

Fig. 3. EW maps obtained using the method described in Sect. 3.2, for coherence lengths of 1 kpc (top), 4 kpc (center), and 10 kpc (bottom). The spaxel distribution used to create the synthetic observations (Sect. 3.3) is superimposed in black.

3.3. Synthetic observations

To use the observational constraints outlined in Sect. 2.4, we need to perform synthetic observations. To this aim, we mimicked the MUSE observations by creating distributions of spaxels that follow the configuration of our three fields combined, correcting the position of each spaxel vertex for the galaxy mass (see Eq. (1)) and choosing a random orientation angle. An example of these distributions is shown in black superimposed on the maps of Fig. 3. We selected only those spaxels for which we have detections in the data.

From each one of these spaxels, we then extract the resulting EW, which is given by the average of all the EW values (both “systems” and “holes”) intercepted by the spaxel, each of them weighted by the fraction of the cell area that overlaps with the spaxel. We exclude all the spaxels within 10 kpc of the center: at these distances, Eq. (2) is indeed only an extrapolation and we see in Fig. 2 that our arctomographic data depart from this prediction. This region therefore does not satisfy the assumption that our data are consistent with the results from Huang et al. (2021). In any case, we verified that including these data does not affect the results reported in Sect. 4. The final EW distribution represents the model output and we can compare this directly with observations in order to constrain C_length, as we describe in detail in the following section.

3.4. Comparison with data and likelihood definition

The likelihood is composed of two different terms, which are related to the two main observational constraints reported in Sect. 2.4, which we call ℒ_scatter and ℒ_fractional; these are explained in detail below. The total likelihood is given by:

$$\begin{array}{c}\hfill ln{\mathcal{L}}_{\mathrm{tot}}=ln{\mathcal{L}}_{\mathrm{scatter}}+ln{\mathcal{L}}_{\mathrm{fractional}}.\end{array}$$(4)

ℒ_scatter is used to compare the variation of EW as a function of the impact parameter, focusing in particular on the intrinsic scatter of this distribution with respect to the best-fit relation, as this represents a proxy for the gas coherence length. This can be seen in Fig. 4, where we show the results of our synthetic observations on the three different model EW maps shown in Fig. 3. In the left panels, we can see the EW distributions extracted from our models as a function of R′, for coherence lengths of 1 (top), 4 (center), and 10 kpc (bottom). The downward triangles depict the “model nondetections” (excluded from our analysis), which we define as EWs lower than the observational detection limit in the arctomographic data, which is equal to 0.2 Å (see Table 1). We note how the predictions of the three models follow the best-fit relation of Eq. (2), as expected by construction. However, the outputs are very different from each other in terms of the intrinsic scatter of the distributions, which is much smaller for smaller C_length. Indeed, if the scale of variation of the EW is smaller than the extension of the source, different EW values will be averaged along the source extension, resulting in a smaller scatter overall than if the coherence scale was larger than the source. This was hypothesized and then demonstrated by Rubin et al. (2018b), where they performed this analysis for the first time using galaxies as extended background sources; here our extended sources are the MUSE spaxels.

Fig. 4. Outputs of the three models presented in Fig. 3, used to infer the likelihood term ℒscatter. From top to bottom, the models have coherence lengths of 1, 4, and 10 kpc. Left: model EW distributions as a function of R′, with the black curve showing the relation of Eq. (2). Downward arrows show spaxels whose predicted EWs are below the detection limit of 0.2 Å and are therefore classified as “model nondetections”. Right: distributions of the differences between the EW values and the corresponding values predicted by Eq. (2), for data (blue) and models (red).

In the right-hand column of Fig. 4, we show in red the binned differences (logEW − logEW_exp) between our synthetic data and the expectations from Eq. (2) (i.e., the black curve in the left panels), for the same models as those shown in the corresponding left panels. The blue histogram shows the distribution of the same quantity for the observational data. Both in models and data, we exclude the values with R′< 10 kpc for the reasons discussed above. We can see how the modeled distributions in the right panels of Fig. 4 differ from each other, as models with larger coherence lengths yield wider distributions. By comparing models and data, we can therefore estimate which coherence length is in better agreement with the observations. We finally define the likelihood term ℒ_scatter as the p-value obtained by comparing the logEW − logEW_exp distributions of models and data with a Kolmogorov–Smirnov (KS) test⁵ (Massey 1951).

ℒ_fractional is instead related to the comparison of the fractional differences. In Fig. 5 we show the outputs of our synthetic observations for the three maps shown in Fig. 3, for three different coherence lengths of 1, 4, and 10 kpc. The left panels show the distribution of ΔEW as a function of the spaxel separation; this distribution is calculated in the same way as for the data (see Sect. 2.4 and the right panel of Fig. 2). We can see how the shape of the distribution varies when changing C_length. To compare with the data, we focus on the small spaxel separations. Indeed, while above we show that our isotropic model is a good representation of our observations, any deviation in the data from this isotropic assumption (see Sect. 5.4) will affect the shape of the ΔEW distribution at large separations more significantly, which will therefore be inconsistent with the predictions of our models. By focusing on small separations, we are more confident that any difference between model and data is entirely due to the cool gas coherence scale.

Fig. 5. Outputs of the three models presented in Fig. 3 and used to infer the likelihood term ℒfractional. From top to bottom, the models have Clength = 1, 4, 10 kpc. Left: distributions of the fractional differences in EW as a function of the spaxel separation (same as right panel of Fig. 2, but for the models). Right: distributions for all the pairs with spaxel separations smaller than 5 kpc, for data (blue) and models (red).

To this end, we selected all the pairs with spaxel separations of less than 5 kpc and the resulting distributions for the three models are shown in red in the right panels of Fig. 5, with the distribution obtained from the data superimposed in blue. We can see how the width of the model distribution increases when increasing the coherence length from 1 to 4 kpc. Indeed, if C_length is significantly smaller than the length of the spaxel, different EW values will be averaged with each other, with the overall effect being a reduction of the differences between the EW values of different nearby spaxels with respect to models where C_length is more similar to the spaxel size. However, models with C_length = 10 kpc have a large amount of low ΔEW, because if the coherence scale is considerably larger than the length of the spaxels, adjacent spaxels will be probing the same EW values, significantly reducing ΔEW. Therefore, this observational constraint seems to represent a good diagnostic to infer the coherence scale of the cool CGM. The likelihood term ℒ_fractional is given by the p-value of a KS test performed between the model and data distributions shown in the right panels of Fig. 5.

3.5. Bayesian analysis

Having entirely determined the likelihood of Eq. (4), we can use it to constrain C_length. We performed a Bayesian analysis – using the nested sampling method, adopting the Python package DYNESTY (Skilling 2004) – over our only free parameter C_length. We used the likelihood function ℒ_tot defined in Sect. 3.4 (Eq. (4)) and a uniform prior on a range that varies from 0.5 to 12 kpc. We note that simply by looking at the observed EW maps of Fig. 1, one can already see large EW variations on adjacent spaxels, which implies that C_length must be of the order of a spaxel in linear size (i.e., a few kiloparsecs; see Table 1). If C_length ≪ 1 kpc or C_length ≫ 10 kpc, we would indeed expect a much smoother EW variation across the halo. Therefore, with this prior we are able to explore the full range of possible cool gas coherence scales that might be in agreement with the data. Once the value of C_length is chosen, to estimate the likelihood in each call, we extract a random realization of the model and the data. The latter is obtained by randomly drawing the values of the observed EWs from a Gaussian distribution centered on the EW value and with a standard deviation equal to its error (see also Rubin et al. 2018b) and by assuming orientation angles for the three fields equal to those used to generate the synthetic observations from the models.

4. Results

4.1. Coherence scale of the cool CGM

Before reporting the results of the Bayesian analysis, we look at how our models compare with the two observational constraints as a function of the assumed coherence length by analyzing the two parts of our likelihood, ℒ_scatter and ℒ_fractional, separately. Indeed, while we use the Bayesian analysis to accurately constrain the range of C_length that better reproduces the observations, it is important to know whether or not the predictions of our models are effectively in agreement with the data.

We explore values of the coherence length varying from 0.5 to 12 kpc and, for each choice of C_length, we create 1000 model and data realizations. For each realization, we finally compare models and data using a KS test, as in Sect. 3.4. The top panel of Fig. 6 shows the results related to the scatter in the EW–R′ relation (ℒ_scatter) and the bottom panel shows those related to the distributions of fractional difference of EW (ℒ_fractional). The black curves show the median values of the KS-test p-value, while the gray bands show the extent of the 32nd and 68th percentiles of the distributions. We consider models with p-values higher than 0.05 (marked by the horizontal dashed line) in agreement with the data.

Fig. 6. Comparison between models with different coherence lengths and the ARCTOMO data. The two plots show the p-values of the KS test between model and data as a function of the coherence length adopted in the model. The black curve shows the median value of the 1000 realizations created for each Clength value (see Sect. 4.1), with the gray band showing the 32nd and 68th percentiles. The horizontal dotted line depicts pvalue = 0.05, meaning that models with a probability higher than this threshold can be considered in agreement with the data. Top: comparison of the scatter in the EW–R′ relation. Bottom: comparison of the EW fractional differences for spaxel separations smaller than 5 kpc. This figure shows that models with 1.5 ≲ Clength/kpc ≲ 7 are consistent with the data.

We can see how the consistency with the data depends on the adopted value of C_length. In particular, the constraint given by the scatter in the EW–R′ relation (top panel of Fig. 6) can be used to put a lower limit on the value of the coherence length, at about 1.3 kpc (where the median p-value is equal to 0.05), which is in line with the findings of Rubin et al. (2018b), who used the same constraint to evaluate the coherence scale of the cool CGM. Using this constraint, the upper limit on C_length is at about 10 kpc or more. The comparison with the fractional differences of EW (bottom panel) allows us to again put a lower limit at about 1.5 kpc and a much tighter constraint on the upper limit of the coherence length, at about 7 kpc, which means a more stringent constraint on C_length. Together, the two constraints define a range of models – with 1.5 ≲ C_length/kpc ≲ 7 – that are in agreement with the ARCTOMO data.

We now discuss the more quantitative results of the Bayesian analysis, where we adopt the full likelihood described by Eq. (4). The main result of this work is shown in Fig. 7, where we report the posterior distribution of the coherence length, the free parameter in our analysis. The vertical dashed lines show the value of the median of the distribution, with the 2σ uncertainties (0.025 and 0.975 quantiles; see also Table 2). We can see that C_length is well defined between 1.4 and 7.8 kpc, with a median value equal to 4.2 kpc. The 2σ limits of the posterior distribution capture the whole range of models that are in agreement with the ARCTOMO data, as shown in the two panels of Fig. 6. These are the tightest bounds on C_length to date. We stress that this is possible thanks to the unique nature of the ARCTOMO data, which allow us to have multiple and extended sources (the MUSE spaxels) around individual galaxies.

Fig. 7. Posterior distribution of the coherence length as found with our Bayesian analysis. The vertical lines show the 2.5, 50, and 97.5 percentiles of the distribution, corresponding to the median value and its 2σ uncertainties, whose values are shown on top. The ARCTOMO data are best reproduced by coherence lengths varying from 1.4 to 7.8 kpc.

4.2. How universal is our result?

Given the complexity of the gaseous galactic halos, one could expect the coherence length of the cool CGM to depend for example on its distance from the central galaxy or to vary across different galaxy halos. In the following sections, we discuss these aspects in more detail, and present our analysis of various subsamples of the original set of data that was used to obtain the results reported above.

4.2.1. Variation of C_length with impact parameter

Above, we show that at least some of the properties of the cool CGM, such as EW and covering fraction (e.g., Huang et al. 2021), seem to be related to its projected distance from the galaxy host (or impact parameter). Although it is important to keep in mind that these are all projected quantities, such trends might be related to the origin and dynamics of the cool CGM and to its connection with the central galaxy. Using our dataset, we investigated whether a similar dependency is also present for the gas coherence length.

We divide our original sample in three bins, for spaxels with R′< 30 kpc, 30 < R′/kpc < 50, and R′> 50 kpc. For each subsample, we then perform the same Bayesian analysis as that performed on the whole sample. The results of this procedure are shown in Fig. 8 and in Table 2. We find a slight trend of increasing coherence length at larger projected distances from the central galaxies, with the median values of the posterior distributions going from 3.9 kpc in the innermost bin to almost 7 kpc in the outermost one. However, it is hard to say whether this trend is physical or simply due to the low-number statistics. The three posterior distributions indeed share the same 2.5 percentile value (1.3 kpc) and the main difference between them is in the distribution width. The large “tails” at large C_length in the middle and right panels are likely due to the fact that most of the spaxels in our sample are within 30 kpc of the galaxy (47 detections for R′< 30 kpc, 21 detections for 30 < R′/kpc < 50, and 7 detections for R′> 50 kpc) and therefore the constraints on C_length at large R′ are very loose⁶.

Fig. 8. Same as Fig. 7, but for the three subsamples obtained by selecting the spaxels with R′< 30 kpc (left), 30 < R′/kpc < 50 (center), and R′> 50 kpc (right).

Given the considerations above, we argue that the dependence of the coherence length on the impact parameter, if any, is only minor. This could imply that the cool CGM originates from similar physical processes across the galactic halo, which tend to create similar gas structures independently of the distance from the central galaxies. More data are needed to study this in greater detail.

4.2.2. Variation of C_length between single galaxies

We performed the same analysis as described in Sect. 3.5 but this time selecting the three halos separately, selecting the spaxels from only one of the three MUSE fields each time. The results of these tests are shown in Fig. 9 and Table 2. We can note how the posterior distributions of the gas coherence length are slightly different in the three cases, indicating that this property can vary across different galaxy halos. While the distributions of the PSZ1 G311.65−18.48 and the SGAS J1226+2152 fields are very similar, with the main difference being the smaller width of the latter (due to the larger number of detections in SGAS J1226+2152), the posterior distribution for the RCS2 032727−132623 field seems to point to a larger gas coherence length in the cool CGM of this galaxy with respect to the other two. However, this difference does not seem to be significant, given that in all cases the posterior distributions are relatively extended and the three coherence lengths are consistent with each other when considering the uncertainties. The slight difference in the C_length of RCS2 032727−132623 with respect to the other fields might be due to the fact that the spaxels in this field are on average at larger impact parameters, possibly reinforcing the idea of a connection between C_length and the impact parameter. A larger number of spaxels for each halo would be needed to more accurately constrain the coherence length of the cool CGM of single galaxies.

Fig. 9. Same as Fig. 7, but selecting the spaxels of the three MUSE fields separately: PSZ1 G311.65−18.48 (left), RCS2 032727−132623 (center), and SGAS J1226+2152 (right).

It is particularly important to note that the three median values of the distributions are all within the range of coherence lengths that we find when combining the three fields (Fig. 7). This suggests that our finding of a range of coherence lengths of between about 1.5 and 8 kpc, which we present in Sect. 4.1, is not driven by one particular galaxy, but likely represents a common property of the cool CGM of these three objects (as we initially assumed). This result seems to confirm that similar galaxies are surrounded by cool CGM with a similar structure. A larger sample of galaxies is needed to confirm this finding.

5. Discussion

5.1. What is the coherence scale?

We have constrained the typical scale over which the absorption strength of the cool CGM around three galaxies at z ≲ 1 does not vary. From an observational point of view, a coherence scale, or length, of a few kiloparsecs means that, when observing the cool CGM around a single galaxy with multiple probes, we do not expect to observe strong variations in Mg II EW for (projected) separations smaller than this scale, while we expect little correlation for larger separations. This can be seen in the right-hand panel of Fig. 2 and has previously been found by other studies (see Sect. 5.2). The fact that we find an entire range of coherence lengths that seem to reproduce our data well (1.4 ≲ C_length/kpc ≲ 7.8), means that it is likely that the exact value of this scale varies across the halo, although we find that there is little (if any) dependence of C_length on the projected distance from the central galaxy. More data are needed to better understand the variation of the coherence length across the halos of galaxies.

The physical nature of this scale is not trivial. As already discussed extensively in Rubin et al. (2018b), a variation in EW may be due not only to differences in the number of clouds intercepted by the line of sight, but also to differences in their kinematics (and partly also in their gas metallicity). Observations of higher spectral resolution (e.g., Lopez et al. 1999; Rauch et al. 2002; Zahedy et al. 2016; Krogager et al. 2018) and/or more sophisticated 3D modeling (see e.g., Li et al. 2023 and Sect. 5.2) are needed to distinguish between different kinematic components along the line of sight and to better understand what property is driving the change in EW. The coherence length that we derived using our empirical models of EW is instead a projected quantity and does not allow us to form conclusions as to the physical properties of the individual CGM clouds. However, we speculate that, regardless of the exact property that drives the change in EW, the coherence scale should represent the size of the cool gas structures over which such properties change. Each one of these structures could be composed of a relatively large cloud with a size comparable to C_length, or of multiple smaller clouds whose properties are related to each other and that cluster over a coherence scale. This length is a direct observable that can be tested with upcoming observational data and with theoretical models (see Sect. 5.3), and it is independent of any physical assumption about the cool CGM.

5.2. Comparison with previous works

The most relevant comparison to our study is the work of Rubin et al. (2018b), which indeed inspired the present exploration. In a similar fashion to the later developments of the present study, these authors used a fiducial model for the cool CGM based on the quasar surveys of Chen et al. (2010b) and Werk et al. (2013), and compared the predictions of this model with their data (a sample of 27 galaxy–galaxy pairs) in order to estimate the gas coherence length. The main constraint was provided by the comparison of the scatter in the EW–R′ distribution of models and data. This coincides with the first part of the likelihood that we used in this work (ℒ_scatter).

In our analysis, we applied two main refinements to the models of Rubin et al. (2018b): (i) we created two-dimensional distributions of absorbers across the halo, also considering the covering fraction as a function of the projected distance from the center, while Rubin et al. (2018b) explored “slices” of the CGM at fixed impact parameters, with a unity covering fraction; and (ii) we used the actual distribution of MUSE spaxels and their shape to extract the synthetic observations, while Rubin et al. (2018b) used more idealized distributions of circular beams to account for the spatial extension of their background galaxies.

Regardless of the above differences, and of the use of a model based on a different survey (Huang et al. 2021), our results are in perfect agreement with those of Rubin et al. (2018b), who estimated a lower limit for the coherence length of ∼1.9 kpc. Our estimate (1.4 < C_length/kpc < 7.8) is consistent with this value, and, in particular, the plot in the top panel of Fig. 6 shows that, using the same observational constraint as Rubin et al. (2018b), we would obtain a very similar result, with very loose constraints on the coherence length upper limit. By using the additional observable given by the EW differences across the spaxel pairs (which was evidently not possible in Rubin et al. 2018b, given the different observational setup) we were also able to put a firm upper limit on the coherence length of the gas, at about 8 kpc. The similarity between our results and those of Rubin et al. (2018b), who used around 30 absorbing galaxies for their analysis, further points to the conclusion that the cool CGM of galaxies with similar properties has a similar “universal” structure (see Sect. 4.2.2).

Various other observational studies have tried to estimate the typical scale of the cool CGM, either directly through the use of multiply lensed quasars, or indirectly, using photo-ionization models. Here, we summarize some of the main results of these works and we refer to Rubin et al. (2018b) for a more complete and extensive review (see their Sect. 3.2). The use of lensed quasars to spectroscopically study the CGM began more than 40 years ago (e.g., Young et al. 1981; Smette et al. 1995) and has been extensively adopted in recent decades in order to investigate the structure of this cool gas (e.g., Rauch et al. 1999, 2002; Lopez et al. 1999, 2005; Ellison et al. 2004; Chen et al. 2014; Zahedy et al. 2016; Rubin et al. 2018a; Okoshi et al. 2021). It is not trivial to compare our findings with those of these works, as different studies use different gas properties, such as column densities of various ions and/or the gas velocity structure (Chen et al. 2014), to derive a coherence scale of the CGM. Even when looking at Mg II absorption, the coherence scale is often based only on the number of coincidences and “anticoincidences” (detection or nondetection of the absorption) in the QSO lines of sight (e.g., Ellison et al. 2004), a method that has been shown to have limitations (Martin et al. 2010). Here instead, in addition to taking full advantage of the arc-tomography spatial sampling, we use the variations of the absorption strength across the halo. The general picture that emerges when consolidating the results from previous studies is that weak Mg II absorbers (EW ≲ 0.4 Å) tend to exhibit a smaller coherence scale (less than 1 kpc) than stronger systems (like those present in the ARCTOMO data), which show a coherence length of a few kiloparsecs or longer, which is in agreement with our findings.

Additional constraints on the structure of the cool CGM have been derived using photo-ionization models (e.g., CLOUDY Ferland et al. 1998, 2013, 2017), which allow estimation of the thickness of the absorbers⁷. Overall, cloud sizes of less than 1 kpc are found at redshifts of z ≳ 2 (e.g., Crighton et al. 2015), while at z ≲ 1, the scales seem to be larger in general, that is, of the order of several up to even hundreds of kiloparsecs (e.g., Werk et al. 2014; Keeney et al. 2017, but see Stern et al. 2016, who assuming that the cool gas spans a variety of densities, found that the cool clouds traced by Mg II in the COS-Halos sample have sizes of only a few tens of parsecs). More recently, Zahedy et al. (2021) applied photo-ionization models to four Lyman-limit systems at redshifts lower than 1 taken from the CUBS survey, and found that the clouds have thicknesses that can vary from 0.01 to 10 kpc, in agreement with our findings (although their mode is around 0.1 kpc). Nevertheless, it is important to note that, due to the various assumptions of ionization models, such as the strength and shape of the adopted UV background (see Acharya & Khaire 2022; Gibson et al. 2022), the cloud sizes inferred from these works remain relatively uncertain. Our estimate of the coherence length is one of the most accurate constraints of the structure of the cool CGM around galaxies at z ≲ 1 to date.

In recent years, more complex and physically motivated models have also been used to interpret observations of the cool CGM and to infer, among other properties, the size of the single clouds. Faerman & Werk (2023) built a model of the cool gas assuming pressure equilibrium between the T ∼ 10⁴ K medium and a hot corona (see Faerman et al. 2020) and allowing for nonthermal pressure support. These authors then compared their models with the COS-Halos data (e.g., Werk et al. 2013) and found that the observed column densities of neutral hydrogen and low and intermediate ions (among which Mg II) are best reproduced by filling factors of the cool gas of the order of 1%. This, together with the observed number of kinematic components (or clouds) along the line of sight, provides an upper limit on the cloud sizes of R_cl ≲ 0.5 kpc. This finding is not in disagreement with our results because, as discussed in Sect. 5.1, we are not modeling single clouds, and each absorbing system in our maps could be made of multiple clouds, which could therefore be smaller than our estimates for the coherence scale. As also mentioned in Faerman & Werk (2023), multiple clouds could be segregated in larger complexes and structures, which would be consistent with our work and with complementary implications from high-resolution data of the high-velocity clouds in the halo of the Milky Way (e.g., Tripp 2022). Finally, using dynamical semi-analytical models of infalling clouds accreted from the intergalactic medium, Afruni et al. (2022) reproduced the observed covering fraction, silicon column density, and kinematics of the cool absorbers detected by the AMIGA project (Lehner et al. 2020) in the halo of the Andromeda galaxy. In their best-fit models, typical cool clouds have radii of about 5 kpc at the virial radius and of less than 1 kpc at distances of a few tens of kiloparsecs from the central galaxy. Although we have no indication in our empirical models of the intrinsic location of the single absorbing clouds, our coherence length going from about 1.5 to 8 kpc seems to be in perfect agreement with the results of these models.

5.3. Implications for hydrodynamical simulations of the CGM

While the coherence scale that we estimate in this work does not directly reflect the physical 3D distribution or the size of the single clouds, it gives an important constraint on the typical scales over which the cool gas is expected to cluster in separate structures (see Sect. 5.1). Our finding therefore represents a very useful observable for comparison with current (and future) theoretical models of the circumgalactic medium.

From the point of view of high-resolution hydrodynamical simulations, single clouds can have sizes of smaller than a parsec because of processes such as fragmentation, or shattering, which might lead to small cloudlets with lengths similar to the gas cooling length (see McCourt et al. 2018). As part of the hypothesis that these small cloudlets are long-lived, one possible picture is that the cool CGM is composed of a fog of these cloudlets, which permeate the whole galaxy halo: this has been proposed especially for the cool CGM at redshifts of z ≳ 2 to explain a small filling factor and simultaneously a covering fraction of close to unity seen in observations (see references in McCourt et al. 2018). However, in order to fully understand the structure of the cool gas,many other physical effects must be considered: small droplets of gas could quickly evaporate into the hot gas that surrounds them because of thermal conduction (Armillotta et al. 2017; Afruni et al. 2023) or they could coagulate to form larger structures (e.g., Gronke & Oh 2020b, 2023), which could later even grow via the induced condensation of the hot gas (e.g., Marinacci et al. 2010; Gronke & Oh 2018; Kooij et al. 2021). The final picture emerging from this type of simulation is therefore that the cool CGM is likely composed of a spectrum of structures, from a fog of cloudlets to more “monolithic” clouds, and the predominance of one or the other depends on the physical properties of the gas (Gronke & Oh 2020b). Our results suggest that, at least in the galaxies probed in this work, the cool CGM presents substantial differences (in Mg II EW) on scales of between 1.4 and 7.8 kpc. We speculate that this finding can be used by hydrodynamical simulations to put upper limits on the larger scales of the cascade of structures that is believed to form the cool gas around galaxies.

Given their parsec-scale resolution, the studies mentioned above can only be focused on small parts of the galactic halo. To have a more extended view, one needs to resort to cosmological and zoom-in simulations, which have been extensively used in the last decade to study the cool CGM (see Faucher-Giguère & Oh 2023, for a recent review). These studies can be very useful for obtaining a full picture of the distribution of the CGM across the halos of galaxies, but have necessarily much lower resolutions with respect to the idealized simulations. Even in the best zoom-in cases (e.g., Peeples et al. 2019), the resolutions can reach a maximum of about 500 pc, which is not sufficient to accurately predict the structure of the cool gas. van de Voort et al. (2019) indeed showed that increasing the resolution in this type of simulation drastically changes the distribution and the covering fraction of the cool gas, indicating no sign of convergence. Therefore, while considerable effort has been made to compare the predictions of this type of simulation with observations of the cool CGM (see e.g., Appleby et al. 2021; DeFelippis et al. 2021), the results have not led to any firm conclusions. With this caveat in mind, it is worth mentioning that, recently, cool clouds with sizes of a few kiloparsecs, possibly in agreement with our findings, have started to be identified even in cosmological simulations such as TNG-50 (Nelson et al. 2020; Ramesh et al. 2023). A CGM coherence scale could be directly extracted from these simulations and compared with the findings of this study.

Future simulations and, in particular, works that aim to link small- and large-scale simulations (e.g., Huang et al. 2022; Fielding & Bryan 2022; Weinberger & Hernquist 2023) will hold the key to a more accurate comparison with the available observational data. Our result regarding the coherence length is one of the fundamental observational constraints that such models should aim to reproduce.

5.4. Limitations and assumptions of this study

An important assumption we made in this study is that the gas is distributed isotropically across the halo, meaning that the average strength of the absorption is the same along rings centered on the host galaxy. This is likely a simplistic assumption. Previous works (e.g., Schroetter et al. 2019) found a segregation of absorbers along the galaxy minor and major axis, suggesting that the cool CGM could be described by a wind plus disk model. In particular, Fernandez-Figueroa et al. (2022) found a dependence of the EW strength on the azimuthal angle for one of the three MUSE fields analyzed in this work, SGAS J1226+2152, although the effect is diluted when considering the three fields together, as we see in Sect. 2.4. Moreover, the presence of satellites in the halo of the main absorbing galaxy might affect the cool gas distribution: part of the absorbing cool gas might indeed come from stripping from these satellite galaxies, as extensively observed in the Local Universe (e.g., Poggianti et al. 2017; Putman et al. 2021), which would contradict our assumption of isotropy.

Nevertheless, by considering only the EW fractional differences at small spaxel separations in our analysis (see Sect. 3.4), we are minimizing the effects due to possible anisotropic distributions, which occur at larger scales. We verified that using larger spaxel separations does not significantly affect our main finding regarding the coherence length range; however, for any given choice of C_length, it reduces the agreement between models and data, indicating a slight level of anisotropy in our data. We therefore acknowledge the importance of also exploring anisotropic models in more detail, but we leave this for future studies, as it is outside the scope of this work. In future works, we will also address more realistic shapes for the gas absorbers, as the squared approximation adopted here is clearly simplistic. However, a more accurate absorber morphology (which is also essentially unknown to date) will likely not affect our findings on the gas coherence scale.

Another main assumption of this work is the creation of our empirical models based on the quasar survey from Huang et al. (2021), assuming that the cool CGM around our galaxies and theirs has similar properties. Our galaxies are at a median redshift of z ∼ 0.8, while those analyzed in Huang et al. (2021) are at a significantly lower redshift (z ∼ 0.21). During this cosmic time, the properties of the cool CGM might change, and while we show that our data are consistent with the EW–R′ relation found by Huang et al. (2021), the intrinsic scatter of this relation could be different, affecting our results on the gas coherence length. However, there remains no evidence of any dependence of this scatter across this range of redshifts. Nielsen et al. (2013) explored the relation between the Mg II absorption strength and the projected distance from the galaxy in the MAGIICAT sample, which is composed of almost 200 galaxies whose redshifts vary from 0.072 to 1.120. These authors found that the scatter of this relation has no dependence on redshift, which implies that our results, although based on a sample of galaxies at lower redshift than our galaxy sample, are most likely robust. To further test this, we ran an additional Bayesian test on C_length, leaving the intrinsic scatter σ used to create our models (see Sect. 3.2) free to vary between 0.05 and 0.5. We find that the posterior distribution of σ has a median at $0.{282}_{-0.124}^{+0.175}$ (2σ uncertainties), which, interestingly, is very similar to the original value of 0.278 (based on the data from Huang et al. 2021), and that the findings on C_length (see Table 2) do not change with respect to our fiducial model.

A limitation of the present study is that our sample is composed of only three galaxies, whose CGM is not necessarily in accordance with the average best-fit relation of Huang et al. (2021), which is based on hundreds of objects. The fact that our EWs lie on top of the Huang et al. (2021) relation seems to suggest that this is a good representation of the cool CGM of our galaxies, at least for R′> 10 kpc; however, including a larger number of halos in our sample might increase the scatter in our data, potentially affecting our findings. The tests described in Sect. 4.2.2 indicate that adding more galaxies to our sample (as long as they have similar properties to those analyzed here) might not influence our results, but a larger and more statistically significant sample of galaxies will be needed to confirm this.

We chose to work with binned spaxels of 0.6″ in order to maximize the number of probes, while still reducing the cross-talk between adjacent spaxels. The (small) amount of cross-talk has the effect of producing smaller EW differences between adjacent spaxels and this might affect our findings. To evaluate the significance of this effect, we ran the same Bayesian analysis, but this time using a 4 × 4 binning for the data (corresponding to 0.8″, see also Lopez et al. 2018; Tejos et al. 2021), which further reduces the seeing-induced cross-correlation. The results, reported in Table 2, are consistent with our fiducial findings, with the median value of the posterior distribution well within the fiducial range of coherence lengths. The main difference resides in the fact that the posterior distribution is much wider (the 97.5 percentile is equal to almost 11 kpc), meaning that the constraints on the coherence scale are looser with respect to our fiducial case, which is due to the smaller number of spaxels. We therefore conclude that, as long as the binned spaxels have dimensions comparable to the observational seeing, our results are relatively insensitive to the exact choice of spatial binning.

Finally, an inherent limitation of our dataset is its incompleteness at EW ≲ 1 Å. This can be seen in Fig. 10, where we show the EW distribution obtained by combining ten different bootstrap realizations of the full data sample (in blue), drawing values from Gaussian distributions centered on the observed EWs and with a standard deviation equal to the observed uncertainties (as in the Bayesian analysis; see Sect. 3.5). We can see how, for EWs weaker than about 1 Å, the data do not exhibit the exponential growth in the distribution expected for Mg II toward quasars (see Zhu & Ménard 2013). On the other hand, the synthetic distribution appears to also be incomplete below 1 Å (red distribution of Fig. 10). This is not surprising as we are excluding all the spaxels with nondetections in the data, which would yield a large number of low EW values in the model (most of them being at large impact parameters). The synthetic distribution was obtained by combining ten different model realizations with a coherence length of 4 kpc, which is consistent with our estimated best-fit value. Even though the shape of this distribution slightly changes depending on the chosen C_length, in all cases, models and data are consistent with each other in the region of weak EWs. This feature also holds when looking at the three fields separately.

Fig. 10. Distribution of Mg II EW in the data (blue) and models with Clength = 4 kpc (red). Both distributions are built by combining ten different realizations of models and data (see Sect. 5.4). Both data and models are incomplete for EW ≲ 1 Å.

Would a higher S/N affect our results? Having spectra with a higher S/N would allow us to increase the number of detections in our data, especially at R′≳50 kpc, where most of our current data show only nondetections. While it is not trivial to predict how a more complete sample would affect our current findings, such a dataset would be crucial in order to more accurately study the possible variation of the coherence length with the distance from the central galaxy, because we would have a much better sampling at large impact parameters. In conclusion, even though a higher S/N is desirable to conduct a more detailed analysis, S/N incompleteness does not affect the results of the current study.

6. Summary and conclusions

In this paper, we investigated the coherence length (C_length) of the cool CGM traced by Mg II absorption. We used arc-tomographic MUSE data of the cool CGM of three different star-forming disk galaxies at redshift ≲1, which allowed us to obtain a sample of Mg II EWs coming from almost 100 spaxels, spanning the galaxy halos up to several tens of kiloparsecs from the center. We developed 2D empirical models that describe the EW distribution and gas-covering fraction based on a recent QSO survey of the cool CGM of low-redshift galaxies. We then compared models and data through a Bayesian analysis in order to find the best coherence length that reproduces both the observed scatter in the distribution of EWs as a function of the impact parameter and the observed distribution of the fractional difference of EWs. Our main findings are as follows:

1. Our data are best reproduced by models with a cool gas coherence length in the range between 1.4 and 7.8 kpc (corresponding to the 2σ limits of our posterior distribution), meaning that the Mg II absorption strength does not vary significantly below these scales.

2. We find that there is no significant trend in the relationship between coherence length and the impact parameter, possibly indicating similar formation mechanisms for the cool gas at different distances from the central galaxies. There is a slight tendency to see larger coherence lengths at larger impact parameters, but the variation is likely due to the low-number statistics, and more data are needed to properly probe this trend.

3. We observe a small variation in the coherence lengths of the three individual galaxy halos, but this variation is well within our fiducial range. This suggests that similar galaxies are likely to be surrounded by CGM with a similar structure.

In conclusion, with this study we are able to put the most stringent constraints on the coherence scale of the cool CGM to date thanks to the use of unique data allowing us to spatially resolve this medium around individual galaxies at z ∼ 1. While here we made use of only three galaxy halos, a future larger sample of such data could help us to confirm and refine the findings of this paper. Moreover, while this analysis is entirely based on the strength of the Mg II absorption lines and on simple empirical models, future studies involving the gas kinematics, the absorption from other ions, and the use of more sophisticated and physically motivated models will help in building a more complete picture of the CGM structure. Such results would inform and constrain our current and future theoretical models of the gaseous halos of galaxies, and would help us to eventually understand how galaxies form and evolve in the Universe.

Acknowledgments

The authors would like to thank the referee for a constructive report. A.A. acknowledges the financial support of the Joint Committee ESO-Chile grant. S.L. acknowledges support by FONDECYT grant 1231187. E.J.J. acknowledges support from FONDECYT Iniciación en investigación 2020 Project 11200263 and the ANID BASAL project FB210003. L.F.B. acknowledges the support of FONDECYT project 23050313.

References

Appendix A: Absorption line fits

In Tables A.1, A.2, and A.3 and in Figs. A.1, A.2, and A.3, we report the results of the automated Gaussian fits (Sect. 2.2) that were used to derive the EWs used in this work, for respectively the PSZ1 G311.65−18.48, RCS2 032727−132623, and SGAS J1226+2152 fields. We show only the results concerning the successful fits, or detections (significance in the EW larger than 2 in both doublet lines), as the nondetections were excluded from our analysis. In the three figures, the normalized MUSE spectra are shown in green, with their respective uncertainties in blue, while the fits of the Mg IIλλ2796, 2803 doublet are reported in red, with the red ticks showing the positions of the line velocity centroids. We note that, for the SGAS J1226+2152 field, the Mg II lines are partially blended with Si IIλ1260 absorption, and so this is modeled by adding an additional Gaussian fit in our analysis for this field (see Tejos et al. 2021).

Fig. A.1. Results of the fitting analysis. Each panel shows a different spaxel (numbered arbitrarily) with Mg II detections in the PSZ1 G311.65−18.48 field. The normalized MUSE spectra are shown in green (uncertainties in light blue). The Gaussian fits of the Mg II doublet are reported in red.

For the reasons explained in Sect. 2.4, we do not consider all the spaxels with R′≲10 kpc⁸, and therefore the spectra from number 21 to 28 in the PSZ1 G311.65−18.48 field and those from 6 to 10 and from 12 to 14 in the SGAS J1226+2152 field are excluded from our analysis. The results of Sect. 4 are based on all the other spectra shown in Figs. A.1, A.2, and A.3. We note that despite fulfilling the criteria of our automated fitting analysis, some of the Gaussian fits appear potentially untrustworthy upon visual inspection. In particular, we consider those of the spectra 12 and 26 in the RCS2 032727−132623 field and that of spectrum 36 in the SGAS J1226+2152 field to be poor fits. To assess the extent to which these untrustworthy spectra influence our results, we ran an additional Bayesyan analysis (Sect. 3.5), excluding the EWs coming from the three spectra mentioned above. The resulting posterior distribution has a median C_length equal to 4.1 kpc and the 2σ limits are equal to 1.3 and 7.7 kpc, which are perfectly consistent with our fiducial findings. We conclude that this (small) set of uncertain fits does not affect the results of this paper.

Fig. A.2. Same as Fig. A.1, but for the RCS2 032727−132623 field.

Fig. A.3. Same as Fig. A.1, but for the SGAS J1226+2152 field. An additional Gaussian that accounts for the Si IIλ1260 absorption is fitted.

All Tables

All Figures

Fig. 1. Maps of EWs of the 2976 Mg II doublet line for all the detections found in the three fields analyzed in this work: PSZ1 G311.65−18.48 (top), RCS2 032727−132623 (center), and SGAS J1226+2152 (bottom). Both the image (observed) plane and the reconstructed absorber plane (dashed-line rectangle) are shown, with the crosses showing the positions of the three identified absorbing galaxies, referred to collectively as G1. In SGAS J1226+2152, the absorption on top of G1 is shown for completeness but not used in the analysis.

In the text

Fig. 2. Main observational constraints extracted from the ARCTOMO data. Left: Mg IIλ2796 rest-frame EW as a function of R′ (Eq. (1)). Arctomographic data are shown in blue, and data from Huang et al. (2021) in orange, with 2σ upper limits denoted by downwards triangles. The black curve shows the best fit of the Huang et al. (2021) sample (Eq. (2)) and the vertical dashed line marks the region where our data tend to diverge from this curve (these data are excluded from our analysis). Right: distribution of the fractional difference of EW in the arctomo data as a function of the spaxel separation. The blue points show the values of all the spaxel pairs, while the blue curve and the gray band show the median values – with the 1σ uncertainties – of a binning of the same distribution (each bin contains 150 values). The black curves show the results for distributions obtained by choosing different orientation angles for the three MUSE fields (see Sect. 2.4 for more details).

In the text

	Fig. 3. EW maps obtained using the method described in Sect. 3.2, for coherence lengths of 1 kpc (top), 4 kpc (center), and 10 kpc (bottom). The spaxel distribution used to create the synthetic observations (Sect. 3.3) is superimposed in black.
In the text

Fig. 4. Outputs of the three models presented in Fig. 3, used to infer the likelihood term ℒscatter. From top to bottom, the models have coherence lengths of 1, 4, and 10 kpc. Left: model EW distributions as a function of R′, with the black curve showing the relation of Eq. (2). Downward arrows show spaxels whose predicted EWs are below the detection limit of 0.2 Å and are therefore classified as “model nondetections”. Right: distributions of the differences between the EW values and the corresponding values predicted by Eq. (2), for data (blue) and models (red).

In the text

Fig. 5. Outputs of the three models presented in Fig. 3 and used to infer the likelihood term ℒfractional. From top to bottom, the models have Clength = 1, 4, 10 kpc. Left: distributions of the fractional differences in EW as a function of the spaxel separation (same as right panel of Fig. 2, but for the models). Right: distributions for all the pairs with spaxel separations smaller than 5 kpc, for data (blue) and models (red).

In the text

Fig. 6. Comparison between models with different coherence lengths and the ARCTOMO data. The two plots show the p-values of the KS test between model and data as a function of the coherence length adopted in the model. The black curve shows the median value of the 1000 realizations created for each Clength value (see Sect. 4.1), with the gray band showing the 32nd and 68th percentiles. The horizontal dotted line depicts pvalue = 0.05, meaning that models with a probability higher than this threshold can be considered in agreement with the data. Top: comparison of the scatter in the EW–R′ relation. Bottom: comparison of the EW fractional differences for spaxel separations smaller than 5 kpc. This figure shows that models with 1.5 ≲ Clength/kpc ≲ 7 are consistent with the data.

In the text

	Fig. 7. Posterior distribution of the coherence length as found with our Bayesian analysis. The vertical lines show the 2.5, 50, and 97.5 percentiles of the distribution, corresponding to the median value and its 2σ uncertainties, whose values are shown on top. The ARCTOMO data are best reproduced by coherence lengths varying from 1.4 to 7.8 kpc.
In the text

	Fig. 8. Same as Fig. 7, but for the three subsamples obtained by selecting the spaxels with R′< 30 kpc (left), 30 < R′/kpc < 50 (center), and R′> 50 kpc (right).
In the text

	Fig. 9. Same as Fig. 7, but selecting the spaxels of the three MUSE fields separately: PSZ1 G311.65−18.48 (left), RCS2 032727−132623 (center), and SGAS J1226+2152 (right).
In the text

	Fig. 10. Distribution of Mg II EW in the data (blue) and models with Clength = 4 kpc (red). Both distributions are built by combining ten different realizations of models and data (see Sect. 5.4). Both data and models are incomplete for EW ≲ 1 Å.
In the text

	Fig. A.1. Results of the fitting analysis. Each panel shows a different spaxel (numbered arbitrarily) with Mg II detections in the PSZ1 G311.65−18.48 field. The normalized MUSE spectra are shown in green (uncertainties in light blue). The Gaussian fits of the Mg II doublet are reported in red.
In the text

	Fig. A.2. Same as Fig. A.1, but for the RCS2 032727−132623 field.
In the text

	Fig. A.3. Same as Fig. A.1, but for the SGAS J1226+2152 field. An additional Gaussian that accounts for the Si IIλ1260 absorption is fitted.
In the text