# Transverse clues on the kiloparsec-scale structure of the circumgalactic medium as traced by C IV absorption

S. Lopez, A. Afruni, D. Zamora, N. Tejos, C. Ledoux, J. Hernandez, T. A. M. Berg, H. Cortes, F. Urbina, E. J. Johnston, et al.

► To cite this version:

S. Lopez, A. Afruni, D. Zamora, N. Tejos, C. Ledoux, et al. Transverse clues on the kiloparsec-scale structure of the circumgalactic medium as traced by C IV absorption. Astronomy & Astrophysics - A&A, 2024, 691, 10.1051/0004-6361/202451200. insu-04838840

HAL Id: insu-04838840

https://insu.hal.science/insu-04838840v1

Submitted on 15 Dec 2024

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# Transverse clues on the kiloparsec-scale structure of the circumgalactic medium as traced by C IV absorption

S. Lopez¹,*, D. A. Afruni¹,²,³, D. Zamora¹, N. Tejos⁴, C. Ledoux⁵, J. Hernandez⁶, T. A. M. Berg⁷,
H. Cortes¹, F. Urbina¹, E. J. Johnston⁸, L. F. Barrientos⁶, M. B. Bayliss⁹, R. Cuellar¹, J. K. Krogager¹⁰,¹¹
P. Noterdaeme¹⁰,¹², and M. Solimano⁸

¹ Departamento de Astronomía, Universidad de Chile, Casilla 36-D, Santiago, Chile

² Kapteyn Astronomical Institute, University of Groningen, Landleven 12, 9747 AD Groningen, The Netherlands

³ Dipartimento di Fisica e Astronomia, Università di Firenze, Via G. Sansone 1, 50019 Sesto Fiorentino, Firenze, Italy

⁴ Instituto de Física, Pontificia Universidad Católica de Valparaíso, Casilla, 4059 Valparaíso, Chile

⁵ European Southern Observatory, Alonso de Córdova, 3107 Vitacura, Casilla 19001, Santiago, Chile

⁶ Instituto de Astrofísica, Pontificia Universidad Católica de Chile, Av. Vicuña Mackenna 4860, 7820436, Macul Santiago, Chile

⁷ NRC Herzberg Astronomy and Astrophysics Research Centre, 5071 West Saanich Road, Victoria, B.C. V9E 2E7, Canada

⁸ Instituto de Estudios Astrofísicos, Facultad de Ingeniería y Ciencias, Universidad Diego Portales, Av. Ejército Libertador 441, Santiago, Chile

⁹ Department of Physics, University of Cincinnati, Cincinnati, OH 45221, USA

¹⁰ French-Chilean Laboratory for Astronomy, IRL 3386, CNRS and U. de Chile, Casilla 36-D, Santiago, Chile

¹¹ Centre de Recherche Astrophysique de Lyon, Université de Lyon 1, ENS-Lyon, UMR5574, 9 Av Charles André, 69230 Saint-Genis-Laval, France

¹² Institut d'Astrophysique de Paris, CNRS-SU, UMR 7095, 98bis bd Arago, 75014 Paris, France

Received 20 June 2024 / Accepted 3 October 2024

## ABSTRACT

The kiloparsec-scale kinematics and density structure of the circumgalactic medium (CGM) is still poorly constrained observationally, which poses a problem for understanding the role of the baryon cycle in galaxy evolution. Here we present VLT/MUSE integral-field spectroscopy ($R \approx 1800$) of four giant gravitational arcs exhibiting $W_0 \gtrsim 0.2$ Å C IV absorption at eight intervening redshifts, $z_{abs} \approx 2.0-2.5$. We detected C IV absorption in a total of 222 adjacent and seeing-uncorrelated sight lines whose spectra sample beams of ("de-lensed") linear size $\approx 1$ kpc. Our data show that (1) absorption velocities cluster at all probed transverse scales, $\Delta r_{\perp} \approx 0-15$ kpc, depending on system; (2) the (transverse) velocity dispersion never exceeds the mean (line-of-sight) absorption spread; and (3) the (transverse) velocity autocorrelation function does not resolve kinematic patterns at the above spatial scales, but its velocity projection, $\xi^{arc}(\Delta v)$, exhibits a similar shape to the known two-point correlation function toward quasars, $\xi^{QSO}(\Delta v)$. An empirical kinematic model suggests that these results are a natural consequence of wide-beam observations of an unresolved clumpy medium. Our model recovers both the underlying velocity dispersion of the clumps ($70-170$ km s$^{-1}$) and the mean number of clumps per unit area ($2-13$ kpc$^{-2}$). The latter constrains the projected mean inter-clump distance to within $\approx 0.3-0.8$ kpc, which we argue is a measure of clump size for a near-unity covering fraction. The model is also able to predict $\xi^{arc}(\Delta v)$ from $\xi^{QSO}(\Delta v)$, suggesting that the strong systems that shape $\xi^{arc}(\Delta v)$ and the line-of-sight velocity components that define $\xi^{QSO}(\Delta v)$ trace the same kinematic population. Consequently, the clumps must possess an internal density structure that generates both weak and strong components. We discuss how our interpretation is consistent with previous observations using background galaxies and multiple quasars as well as its implications for the connection between the small-scale kinematic structure of the CGM and galactic-scale accretion and feedback processes.

**Key words.** galaxies: evolution – galaxies: halos – galaxies: high-redshift – intergalactic medium – quasars: absorption lines

## 1. Introduction

The widespread presence of metals in the diffuse inter-galactic and circumgalactic media (IGM and CGM, respectively; e.g., Cowie et al. 1995; Ellison et al. 2000; Simcoe et al. 2004; Songaila 2005; Ryan-Weber et al. 2006; Becker et al. 2019; Cooper et al. 2019; D'Odorico et al. 2023; Bordoloi et al. 2024) suggests a continuous metal enrichment since $z > 4$ that remained constant until $z = 2$ (McQuinn 2016; D'Odorico et al. 2022; Galbiati et al. 2023). Produced in galaxies (e.g., Pettini et al. 2001; Adelberger et al. 2005; Lofthouse et al. 2023; Banerjee et al. 2023) and their environments

(Shen et al. 2012), CGM metals enter a baryon cycle powered by galactic-scale feedback and re-accretion mechanisms that are believed to regulate star formation and ultimately galaxy evolution (e.g., Kereš et al. 2005; Oppenheimer & Davé 2008; Faucher-Giguère et al. 2011; Christensen et al. 2016). Both the outflowing (Rupke et al. 2005; Oppenheimer & Davé 2006) and inflowing (Nelson et al. 2015; Fielding et al. 2017; Faucher-Giguère & Oh 2023) gas may be clumpy, resulting in poor mixing of the metals on "small scales" (Schaye et al. 2007), defined here as $\le 1$ kpc. Thus, the CGM small-scale kinematics and spatial structure must be intimately connected to the baryon cycle of galaxies. Constraining the former observationally has become a cornerstone for understanding the latter and validating hydro simulations of increasingly higher

<footnote>* Corresponding author; slopez@das.uchile.cl</footnote>
resolution (e.g., Oppenheimer & Davé 2008; Wiersma et al. 2010; Cen & Chisari 2011; Rahmati et al. 2016; Bird et al. 2016; Finlator et al. 2020; Hummels et al. 2019; Peeples et al. 2019; Marra et al. 2024).

Measuring the clumpiness of the CGM is difficult because this medium is diffuse, and the rare bright background sources (e.g., quasars) required to detect it in absorption do not provide transverse sampling. The only option to directly probe the transverse dimension seems to rely on the even scarcer multiple background sources. In fact, resolved spectroscopy of lensed quasars and galaxies already have provided evidence that the $T = 10^4$ K enriched gas is clumpy on kiloparsec scales (Rauch et al. 1999, 2001a,b; Ellison et al. 2004; Lopez et al. 2005, 2007, 2018, 2020; Chen et al. 2014; Zahedy et al. 2016; Rubin et al. 2018a,b; Péroux et al. 2018; Krogager et al. 2018; Kulkarni et al. 2019; Zahedy et al. 2019; Mortensen et al. 2021; Bordoloi et al. 2022; Afruni et al. 2023).

An alternative way to probe the CGM structure on various scales is through measuring line-of-sight and transverse velocity clustering. This technique has been applied on different data sets, for instance, (i) the auto-correlation function of absorbers either along the line-of-sight at high spectral resolution (Steidel 1990; Pichon et al. 2003; Scannapieco et al. 2006; Boksenberg & Sargent 2015; Rauch et al. 1996; Fathivavsari et al. 2013) or transversely using lensed or multiple quasars in general (Rauch et al. 2001a; Coppolani et al. 2006; Tytler et al. 2009; Martin et al. 2010; Mintz et al. 2022; Maitra et al. 2019; Gontcho A Gontcho et al. 2018; Dutta et al. 2024; Hennawi et al. 2006); (ii) the absorber-galaxy cross-correlation either along the line-of-sight using background galaxies (Steidel et al. 2010; Turner et al. 2017) or transversely (Adelberger et al. 2005; Lofthouse et al. 2023; Banerjee et al. 2023; Galbiati et al. 2023); and (iii) other cross-correlations such as quasar-absorber (Hennawi & Prochaska 2007; Vikas et al. 2013) or outflow-absorber (Rauch et al. 2001b; Theuns et al. 2002). The broad picture that has emerged is that metals and galaxies trace the same overdensities. However, none of the observations can really disentangle the absorbing galaxy from the absorption system. Only a handful of them address the line-of-sight kinematics (e.g. Turner et al. 2017) and how this is supposed to be entangled with the spatial structure (Stern et al. 2016), and even fewer have really been able to measure the level of transverse structure (Rauch et al. 1999, 2001a,b).

In this article, we take advantage of multiplexed spectroscopy of giant gravitational arcs (hereafter “ARCTOMO<fnref n="1" /> data”) to measure the velocity clustering of intervening triply ionized carbon (C IV) across kiloparsec scales. The strong and easy-to-identify C IV λλ1548, 1550 doublet is arguably the most sensitive metal tracer of both the cool and warm CGM at high redshifts (Chen et al. 2001; Bordoloi et al. 2014; Rudie et al. 2019). On the other hand, ARCTOMO data have demonstrated the potential to add a wealth of unique and new spatial information on CGM scales (Lopez et al. 2018, 2020; Tejos et al. 2021; Fernandez-Figueroa et al. 2022; Afruni et al. 2023), hence the timely combination. We build on a blind survey of $z \approx 2$ C IV in all available ARCTOMO fields (originally targeted for $z \approx 1$ Mg II) that resulted in eight intervening C IV systems toward four arcs.

The paper is organized as follows. In Sect. 2, we describe the observations, the data reduction, and the spectrum extraction. Then in Sect. 3 we describe the automated C IV identification and line-profile fitting used, paying special attention to survey

completeness, and in Sect. 4 we address the subtleties of dealing with lensed fields. We used these data for two kinds of analysis on C IV kinematics: First, in Sect. 5, we carry out a direct assessment of line-of-sight and transverse kinematic properties, and secondly, in Sect. 6, we measure transverse velocity clustering and compare it with quasar line-of-sight observations. Thereupon, in Sect. 7, we present a gas-kinematics model that explains the results, reproduces the arc signal out of the quasar kinematics, and predicts independent observations. Finally, in Sect. 8, we discuss the implications of our findings, concluding with a summary of the results in Sect. 9.

Throughout the paper we use a ΛCDM cosmology with the following cosmological parameters: $H_0 = 70$ km s⁻<fnref n="1" /> Mpc⁻<fnref n="1" />, $\Omega_m = 0.3$, and $\Omega_\Lambda = 0.7$. We also use the standard notations $N(\mu, \sigma^2)$ for a normal distribution with mean $\mu$ and variance $\sigma^2$ and $U_{[a,b]}$ for a continuous uniform distribution with support [a,b].

## 2. ARCTOMO data

### 2.1. Observations

The four arcs were observed with the Multi Unit Spectroscopic Explorer (MUSE; Bacon et al. 2010) mounted on UT4 ( Yepun) at the Very Large Telescope (VLT) in Paranal, Chile, as part of programs 297.A-5012 (PI Aghanim), 098.A-0459 (PI Lopez), and 0103.A-048 (PI Lopez). The observations were conducted in service mode with the MUSE wide field mode (WFM), which provides a field of view of $1' \times 1'$ sampled at 0.′/2/spaxel. They used either non-adaptive optics and a nominal wavelength range (NOAO-N) or adaptive optics and an extended wavelength range (AO-E), depending on the program (Table 1). These setups provide spectral coverage of ≈4700–9300 Å and ≈4600–9300 Å respectively<fnref n="2" />, and a resolving power ranging from $R \approx 1770$ at 4800 Å to $R \approx 3590$ at 9300 Å. Exposure times varied between 0.8 and 4.2 hours. The requested image quality resulted in final point spread functions (PSFs) ranging from 0.′/74 to 0.′/84, depending on the targeted field.

### 2.2. Data reduction

The data reduction was carried out using the ESO MUSE pipeline (v2.6, Weilbacher et al. 2020) in the ESO Recipe Execution Tool (EsoRex) environment (ESO CPL Development Team 2015). The master bias, flat field, and wavelength calibrations for each CCD were created from the associated raw calibrations, and they were applied to the raw science and standard-star observations as part of the pre-processing steps. Flux calibration was carried out using the standard star observations from the same nights as the science data, and the sky continuum was measured directly from the science exposures and subtracted off. The reduced data for each exposure were then stacked to produce the combined datacube, with a wavelength solution calibrated to vacuum. Any residual sky contamination was removed using the Zurich Atmosphere Purge code (ZAP, Soto et al. 2016). Finally, the datacubes were matched to the WCS of the corresponding HST data (references in Table 1). White-image stamps of the arcs are shown in the left-hand column of Fig. 1.

<footnote n="1">¹ https://sites.google.com/view/arctomo/home</footnote>

<footnote n="2">² A gap between 5760–6010 Å is present due to the contamination produced by the sodium laser used in AO.</footnote>
<table><caption>**Table 1.** Summary of targets and VLT/MUSE observations.</caption><thead><tr><th>Field name</th><th>RA<br>[h:m:s]</th><th>Dec<br>[°:′:′′]</th><th>z<sub>em</sub></th><th>Total Exp.<br>time [s]</th><th>Instrument<br>Mode</th><th>Effective<br>PSF <sup>(a)</sup> [′′]</th><th>Program ID</th><th>Lens<br>Ref. <sup>(b)</sup></th></tr></thead><tbody><tr><td>SGAS J0033+0242</td><td>00:33:41.52</td><td>+02:42:27.0</td><td>2.39</td><td>11 × 700</td><td>WFM-NOAO-N</td><td>0.84</td><td>098.A-0459(A)</td><td>(1)</td></tr><tr><td>SGAS J1527+0652</td><td>15:27:44.74</td><td>+06:52:19.4</td><td>2.76</td><td>25 × 600</td><td>WFM-AO-E</td><td>0.78</td><td>0103.A-0485(B)</td><td>(2)</td></tr><tr><td>PSZ1 G311.65–18.48</td><td>15:50:07.46</td><td>−78:11:23.2</td><td>2.37</td><td>2 × 1483</td><td>WFM-NOAO-N</td><td>0.82</td><td>297.A-5012(A)</td><td>(3)</td></tr><tr><td>SGAS J2111–0114</td><td>21:11:19.01</td><td>−01:14:30.4</td><td>2.86</td><td>15 × 606</td><td>WFM-AO-E</td><td>0.74</td><td>0103.A-0485(B)</td><td>(4)</td></tr></tbody></table>

Notes. <sup>(a)</sup>Point spread function at ≈5000 Å measured from fitting 2D Gaussian profiles to stars in the final combined MUSE cubes. <sup>(b)</sup>Lens model references: (1) Fischer et al. (2019), (2) Sharon et al. (2020), (3) Lopez et al. (2020), and (4) Sharon et al. (2020).

## 2.3. Binned spectra extraction

We extracted and combined the cube spectra optimally using 3 × 3 spaxel apertures (so binned spaxels are non-overlapping squares of 0′′6 on a side), a size deemed sufficiently large to minimize cross-talk between adjacent spaxels due to seeing smearing, while maximizing spatial sampling (e.g., Lopez et al. 2018). We show in Appendix C that this choice does not bias our results. We considered an arc binned spectrum to sample an independent “arc sight line” or “arc beam.” By construction, these spectra have heterogeneous signal-to-noise ratios. Throughout the article, we apply appropriate completeness corrections to deal with each arc’s particular selection function.

## 3. Absorption line analysis and definitions

We defined a C IV “ARCTOMO system,” or simply a “C IV system,” as having significant C IV absorption detected in at least one arc sight line. We borrowed the concept of “system” from the quasar absorption lines technique despite the fact that here we deal with extended sources so that a system could be probed by several adjacent sight lines. In this article, we refer to intervening systems, $c|z_{abs} - z_{em}|/(1 + z_{em}) > 3000 \text{ km s}^{-1}$, where $z_{abs}$ and $z_{em}$ are the system absorption and the source emission redshifts, respectively. The C IV system identification and line profile fitting were performed separately and automatically. We outline these two steps below.

### 3.1. C IV system identification

To automatically find candidate systems, binned spectra were pre-selected within a manually selected mask that contained the arc and sufficient sky around it. We identified candidate C IV doublets using the template-matching algorithm described in Noterdaeme et al. (2010) and Ledoux et al. (2015). The search runs over all masked spectra along the available redshift path (see Appendix B for details). A total of eight candidate C IV systems were found toward four arcs at redshifts $z_{abs}$ between 2 and 2.5.

### 3.2. Line profile fitting

In a second round, an automated fit of a double Gaussian plus a local continuum is ran over all masked spectra having S/N > 1 at each candidate C IV redshift. The Gaussians have tied doublet separation; free common width (constrained by FWHM<sub>obs</sub> ≥ 2 pixels); free amplitudes $A^{1548}$ and $A^{1550}$, such that $1 \le A^{1548}/A^{1550} \le 2$; and free velocity within ±4000 km s<sup>-1</sup> of

$z_{abs}$<sup>3</sup>. Only one C IV doublet is considered in each velocity window. We note that at MUSE spectral resolution, the absorption profiles are dominated by the unresolved LOS kinematics and not by the doublet ratios; hence, $A^{1548}/A^{1550} \approx 1$ can occur even if the lines do not reach zero level. This has also been observed in quasar spectra of similar quality (Cooksey et al. 2013).

Best-fit parameters and their 1-sigma errors were obtained for rest-frame equivalent width ($W_0 \equiv W_0^{1548}$; $\delta_{W0}$), doublet velocity ($v$; $\delta_v$), and line width ($\sigma_{obs}$; $\delta_\sigma$). We adjusted $z_{abs}$ so that the median velocity per system $\langle v \rangle = 0 \text{ km s}^{-1}$. The errors $\delta_v$ and $\delta_\sigma$, both relevant here, typically range within 10–25 km s<sup>-1</sup> (68% level; Fig. B.1). For a fit to be considered successful (i.e., a “detection”), we required a significance $W_0/\delta_{W0} \ge 2$ on both doublet lines and $\delta_v < 35 \text{ km s}^{-1}$ (≈ 1/2 pixel). We note that the rather loose S/N pre-selection ensured scanning of all arc spaxels, but the final “decision” on detections was taken autonomously by the fitting algorithm. Unsuccessful fits became “non-detections,” for which $2-\sigma$ upper limits were computed using $W_0 = 2 \times \text{FWHM}/\langle S/N \rangle/(1+z)$.

### 3.3. C IV sample

A total of 533 binned spectra were processed, resulting in 222 detections grouped in eight C IV systems. We nicknamed systems with a short name for the arcs and vowels for their incidence. Various maps and properties of these systems are presented in Fig. 1 and Table 2, respectively. Fitted profiles for all systems are shown in Figs. A.1–A.8. This sample comprises a unique data set of spatially resolved and significant C IV detections at redshifts 2.0–2.5.

### 3.4. Equivalent width completeness

Heterogenous S/N in a given system naturally leads to an incomplete distribution of fitted parameters. To account for each system’s $W_0$-completeness, we estimated the false-negative rate (FNR) of the fitting procedure. To this end, synthetic spectra were created by replacing real absorption features with randomly selected flux (at the continuum) from neighbor pixels. A synthetic doublet was injected randomly inside the velocity window where the fit was initially performed, and $W_0$ was computed as in Sect. 3.2. This test was conducted iteratively for each binned spectrum where a feature was found and for a certain range of equivalent widths. Whenever the fitting algorithm failed to detect the inserted doublet, the instance was flagged as a false negative. We computed the FNR as the number of false negatives divided by the sum of the total number of false negatives

<footnote n="3">³ Velocity window chosen according to the correlation analysis explained in Sect. 6.</footnote>
<figure><img src="image_2.png" /><figcaption>**Fig. 1.** MUSE arc images and C IV system maps. Columns from left to right: Column (1): White-image stamps of the gravitational arcs. Column (2): 0.6' × 0.6'-binned spaxels with C IV detections, colored according to their significance (Sect. 3.2). Spaxel coordinates as in Figs. A.1–A.8. In both columns, the PSF is represented by the lower-right circle. Columns (3), (4), and (5): Absorber-plane reconstructed maps of equivalent width ($W_0^{1548}$), velocity, and absorption spread, respectively. In all panels, the scale is indicated. In Columns (2)–(5), the triangles indicate $W_0 \le 0.2$ Å (2-σ) non-detections.</figcaption></figure>
<table><caption>**Table 2.** Summary of C IV systems.</caption><thead><tr><th>System</th><th>z<sub>abs</sub><br>(1)</th><th># of<br>spec.<br>(2)</th><th># of<br>det.<br>(3)</th><th>⟨S/N⟩<br>(4)</th><th>σ<sub>⊥</sub><br>[km/s]<br>(5)</th><th>⟨σ<sub>‖</sub>⟩<br>[km/s]<br>(6)</th><th>W<sub>0</sub>-range<br>[Å]<br>(7)</th><th>Δr<sub>⊥</sub>-range<br>[kpc]<br>(8)</th><th>Area<br>[kpc<sup>2</sup>]<br>(9)</th></tr></thead><tbody><tr><td>1527 a</td><td>2.0553</td><td>45</td><td>32</td><td>6.9</td><td>40.6±4.0</td><td>131.5±2.9</td><td>1.2–2.2</td><td>0.2–13.9</td><td>117.6</td></tr><tr><td>1527 b</td><td>2.1624</td><td>49</td><td>14</td><td>8.8</td><td>65.7±5.7</td><td>132.8±6.8</td><td>0.4–1.7</td><td>0.7–12.3</td><td>52.1</td></tr><tr><td>1527 c</td><td>2.5430</td><td>55</td><td>23</td><td>10.2</td><td>22.4±3.4</td><td>61.7±4.4</td><td>0.2–1.4</td><td>0.2–14.4</td><td>65.5</td></tr><tr><td>G311 a</td><td>2.1894</td><td>181</td><td>79</td><td>9.5</td><td>30.4±2.1</td><td>62.0±3.1</td><td>0.3–1.5</td><td>0.0–10.5</td><td>39.5</td></tr><tr><td>0033 a</td><td>2.2814</td><td>50</td><td>19</td><td>7.4</td><td>74.6±5.2</td><td>134.7±5.5</td><td>0.7–1.9</td><td>0.1–8.7</td><td>9.1</td></tr><tr><td>0033 b</td><td>2.3356</td><td>50</td><td>18</td><td>7.6</td><td>57.7±4.9</td><td>127.4±5.6</td><td>1.0–2.6</td><td>0.1–6.4</td><td>9.3</td></tr><tr><td>2111 a</td><td>2.4367</td><td>50</td><td>12</td><td>5.1</td><td>64.9±7.2</td><td>86.4±8.4</td><td>0.4–1.3</td><td>0.2–4.4</td><td>5.9</td></tr><tr><td>2111 b</td><td>2.5370</td><td>53</td><td>25</td><td>3.6</td><td>54.4±5.7</td><td>83.3±6.3</td><td>0.7–1.7</td><td>0.1–9.6</td><td>15.4</td></tr></tbody><tfoot><tr><td>Total</td><td></td><td>533</td><td>222</td><td></td><td></td><td></td><td></td><td></td><td>314.4</td></tr></tfoot></table>

**Notes.** Table columns: (1) Absorption redshift, (2) number of binned spectra surveyed, (3) number of C IV detections, (4) median pixel S/N, (5) transverse velocity dispersion, (6) median absorption spread, (7) range of C IV rest-frame equivalent widths, (8) range of transverse distances between spaxel centers, and (9) total area of C IV detections.

<figure><img src="image_3.png" /><figcaption>**Fig. 2.** Rest-frame equivalent width distributions. Upper panels: Completeness of $W_0$ and exponential fit. $W_0 \ge 2$ Å produces blended doublet lines and are not used in the fit. Lower panels: Measured and completeness-corrected $W_0$ distributions (filled and unfilled histograms, respectively). We note that system G311 a has a different y-scale.</figcaption></figure>

atives and true positives. The FNR was computed as a function of $W_0$ and modeled with an exponential function of the form $\text{FNR}(W_0) = \beta \exp(\alpha W_0)$. The $W_0$-completeness level was defined as $1 - \text{FNR}(W_0)$ (Fig. 2, upper panels) and used in Sect. 6 to correct the observed $W_0$ distribution (ibid., lower panels). Effects from S/N on other parameters are addressed in Sect. 5 and Appendix D.

## 4. Spatial information and de-lensing

Clearly the most innovative – and we assert most powerful – feature of the present data set is its spatial coherence and resolution. In the so-called image plane (Grossman & Narayan 1988), that is, the geometry recorded by the instrument, the present binned spaxels are non-overlapping squares (of 0.6 on a side); however, arcs are a consequence of strong lensing, so the observed spaxels do not represent the real geometry of the so-called absorber

plane at $z_{\text{abs}}$ (Lopez et al. 2018). To reconstruct the absorber plane, we used parametric lens models built for each field using Lens stool Tool (Jullo et al. 2007) and based on available HST imaging (for details, we refer the reader to the respective publications indicated in Table 1). We used these models to calculate the deflection angles and reconstruct the observed spaxel grid at the absorber plane by applying the lens equation to the spaxel vertices. The center of a spaxel in the absorber plane was taken to be the center of a spaxel in the image plane de-lensed back to the absorber plane.

Columns (3)–(5) of Fig. 1 show $2'' \times 2''$ stamps of the reconstructed absorber plane at each system’s redshift. Displayed are maps of $W_0$, $v$, and absorption spread. As in column (2), these maps consist of detections only, except for the triangles, which indicate $W_0 \le 0.2$ ($2-\sigma$) upper limits. This threshold roughly corresponds to a 50% $W_0$-completeness on average. The $W_0$ maps suggest spatial coherence on ~10 kpc scales down to $W_0 \approx 0.3$ Å.

The figure highlights the great advantage of the present data to reach down to sub-kiloparsec resolution at $z \approx 2$. This is due basically to the absorber redshift being close to $z_{\text{source}}$, where beam separations approach zero. The drawback is that in some fields, the de-lensing produces heavily distorted grids, due to high magnification close to critical curves. There is also some level of overlapping spaxels, and in some cases their vertices even appear flipped. Both effects are nonphysical and likely reflect the spatial inaccuracy of the lens model at high magnifications. In this paper, we use average separations between spaxel centers and average spaxel areas only; we neglect the shape of the reconstructed spaxels.

The median spaxel area of C IV detections was found to be $= 0.65$ kpc$^2$, so the median linear scale is ~0.81 kpc. By system, the spaxel areas range from 0.5 to 3.6 kpc$^2$. The total areas per system surveyed are small compared with strong C IV “sizes” of ~100 kpc (linear) obtained with different methods (Rauch et al. 2001a; Steidel et al. 2010; Martin et al. 2010; Rudie et al. 2019; Hasan et al. 2022).

## 5. Line-of-sight versus transverse velocity dispersion

For each system, we defined two velocity dispersions: (1) a “transverse dispersion,” $\sigma_{\perp}$, computed as the standard deviation
<table><caption>**Table 3.** Ancillary data.</caption><thead><tr><th>System</th><th>λ<sup>1548</sup><br>[Å]<br>(1)</th><th>⟨δ<sub>v</sub>⟩<br>[km/s]<br>(2)</th><th>σ<sub>inst</sub><br>[km/s]<br>(3)</th><th>σ<sub>⊥</sub><sup>MC</sup><br>[km/s]<br>(4)</th><th>α<br>(5)</th></tr></thead><tbody><tr><td>1527 a</td><td>4730.2</td><td>17.5</td><td>67.2</td><td>17.2</td><td>-3.0</td></tr><tr><td>1527 b</td><td>4896.0</td><td>20.8</td><td>64.9</td><td>27.2</td><td>-3.9</td></tr><tr><td>1527 c</td><td>5485.3</td><td>12.2</td><td>57.9</td><td>18.3</td><td>-5.4</td></tr><tr><td>G311 a</td><td>4937.8</td><td>14.4</td><td>64.3</td><td>19.8</td><td>-5.1</td></tr><tr><td>0033 a</td><td>5080.2</td><td>21.0</td><td>62.5</td><td>24.9</td><td>-2.3</td></tr><tr><td>0033 b</td><td>5164.2</td><td>17.8</td><td>61.5</td><td>19.2</td><td>-2.0</td></tr><tr><td>2111 a</td><td>5320.7</td><td>22.8</td><td>59.7</td><td>28.8</td><td>-2.6</td></tr><tr><td>2111 b</td><td>5476.0</td><td>20.2</td><td>58.0</td><td>22.1</td><td>-1.6</td></tr></tbody></table>

**Notes.** Table columns: (1) Wavelength of C IV zero-point velocity; (2) median velocity uncertainty (the individual velocity errors δv are used in Sect. 5 to carry out the bootstrap analysis on σ⊥ and in Sect. 6 to produce v = 0 mock catalogs. These are listed in Table A.1 and their distributions displayed in Fig. B.1); (3) instrumental profile width at λ<sup>1548</sup> (σ<sub>inst</sub> is used to compute σ‖ in Sect. 5. σ<sub>inst</sub> was obtained at each system's wavelength by interpolating the values measured by Mentz et al. (2016) using sky lines); (4) spurious σ⊥ (Sect. 5 and Appendix D); (5) slope of exponential fit to W<sub>0</sub>-completeness (Sect. 3.4).

of v of all spaxels along the arc (i.e., analogous to the projected velocity dispersion of galaxy clusters); and (2) a “parallel dispersion” (also “absorption spread”), defined as σ‖ = $\sqrt{\sigma_{\text{obs}}^2 - \sigma_{\text{inst}}^2}$, where σ<sub>inst</sub> = σ<sub>inst</sub>(λ) is the instrument spectral resolution (Table 3).

Figure 3 shows the distributions of v and σ‖ (blue and yellow histograms respectively). The vertical dashed lines indicate the median v and median σ‖ (≡ ⟨σ‖⟩) Each system's σ⊥ is indicated by the length of the blue horizontal line. One feature that stands out is that ⟨σ‖⟩ > σ⊥ in all systems, despite the fact that both dispersions vary across systems. Also, the velocity histograms are roughly symmetric. This suggests some level of Gaussianity and little contamination by line-of-sight outliers.

The measured σ⊥ could be biased high due to velocity outliers. On the other hand, sample incompleteness affects the low-W<sub>0</sub> end, and since W<sub>0</sub> and σ‖ are correlated by the Gaussian line fitting, ⟨σ‖⟩ may also be overestimated. To carry out a more robust comparison (system-by-system) between transverse and line-of-sight dispersions, we bootstrapped each sample by creating 1000 synthetic realizations of velocities and dispersions. These are drawn randomly from N(v, δ<sub>v</sub><sup>2</sup>) and N(σ‖, δ<sub>σ</sub><sup>2</sup>), where v, σ‖, and their respective variances correspond to the measured values. From the bootstrapped distributions per system, we took the median and standard deviations as corrected values. A Monte Carlo (MC) analysis (Appendix D) showed that the S/N selection function introduces a spurious dispersion of ≈17–29 km s<sup>-1</sup> (system dependent) in the spatial direction. This value was subtracted in quadrature, although it has only a marginal effect on σ⊥. As for ⟨σ‖⟩, the same MC analysis showed that our survey misses ≲10% of the systems detectable with MUSE (up to ~30% in 2111 a and 2111 b). It is hard to assess the effect of this bias on ⟨σ‖⟩, but it would translate into a less than 5% only upward bias in the “worst-case” scenario that all missed systems have σ‖ ≈ 0 km s<sup>-1</sup>.

Figure 4 shows the median absorption spread by system, ⟨σ‖⟩ as a function of σ⊥ with corresponding 1σ errors. We observed that ⟨σ‖⟩ > σ⊥ remains, and we note that it cannot be induced by the corrections described above. The relative spaxel areas probed

by each system (circle sizes) do not seem to correlate with any of the above measurements (nor the total areas probed). Conversely, the median W<sub>0</sub> broadly correlates with ⟨σ‖⟩, but this is expected based on the former results from fitting the velocity spread.

# 6. Transverse auto-correlation of velocities

Whilst the simple comparison between σ‖ and σ⊥ can have profound implications on the origin of the enriched gas and its kiloparsec-scale substructure, it is sensitive to outliers and hard to interpret. In this section, we use a more sophisticated statistical tool by measuring coherence in pairwise velocity differences and separations. To this end, we defined a transverse velocity auto-correlation function, computed it on a system-by-system basis, and tested the possibility that this function (1) unveils the spatial structure and (2) is related to the line-of-sight velocity correlation measured in quasar spectra (e.g., Steidel 1990; Petitjean & Bergeron 1994; Pichon et al. 2003; Scannapieco et al. 2006; Boksenberg & Sargent 2015; Rauch et al. 1996; Fathivavsari et al. 2013).

## 6.1. Definition

We defined a transverse velocity auto-correlation function, ξ(Δr⊥, Δv), as the excess probability of finding a given pairwise C IV absolute velocity difference, |Δv|, and we binned the spaxel separation on the reconstructed absorber plane, Δr⊥. We used the so-called natural estimator (e.g., Kerscher et al. 2000):

$$ \xi(\Delta r_{\perp}, \Delta v) = \frac{\langle DD \rangle}{\langle RR \rangle} - 1, \qquad (1) $$

where ⟨DD⟩ and ⟨RR⟩ are data-data and random-random pair averages, respectively. The pairs were created from data and random catalogs described below. To account for errors in ξ(Δr⊥, Δv), measured velocities were randomized with standard deviation δv, and ⟨DD⟩ was computed from the median values in the (Δr⊥, Δv) bin. The 1-σ errors result from adding in quadrature the standard deviation of the median statistics (“measurement” noise) with the square root of the number of data pairs (“shot” noise). See Appendix F.1 for more details.

## 6.2. Data catalogs

Each system's data catalog consists of entries of (1) a C IV sight line de-lensed RA-DEC coordinate; (2) a rest-frame equivalent width W<sub>0</sub> and its 1σ error, δ<sub>W</sub>; (3) a velocity v and its 1σ error, δ<sub>v</sub>; and (4) the same continuum S/N level used to pre-select spaxels. Data catalogs were used to create a list of pairwise velocity differences, |Δv|, and separations on the sky, Δr⊥. We assumed that the absorption signal is spatially independent from one spaxel to another.

## 6.3. Random catalogs

Random catalogs must account for each system's S/N selection function (e.g., Tejos et al. 2014), both in the spectral and transverse directions. Here, we followed a similar procedure as in Martin et al. (2010) and Mintz et al. (2022), albeit with important differences due to our particular data type.

For each system, we first created a new catalog filled with five repeated copies of the data catalog (including non-detections). This has the advantage of preserving the survey
<figure><img src="image_4.png" /><figcaption>**Fig. 3.** Distribution of C IV doublet velocity ($v$) and absorption spread ($\sigma_{\parallel}$) by system. The blue dashed line marks the zero-point velocity, and the orange dashed line indicates the median absorption spread ($\langle\sigma_{\parallel}\rangle$). The blue horizontal solid line has a size equal to the transverse velocity dispersion ($\sigma_{\perp}$); its blue side is arbitrarily placed at zero velocity for a visual comparison between $\sigma_{\perp}$ and $\langle\sigma_{\parallel}\rangle$. The data represented in this figure has not been corrected for incompleteness.</figcaption></figure>

<figure><img src="image_5.png" /><figcaption>**Fig. 4.** System-by-system median absorption spread versus transverse velocity dispersion. Errors result from a bootstrapping analysis (Sect. 5). Circle areas are proportional to the absorber-plane mean spaxel area (in square kiloparsecs), and colors indicate the median rest-frame equivalent width per system (in Å). $\sigma_{\parallel}$ has been corrected for the instrumental profile width (Sect. 5), and $\sigma_{\perp}$ has been corrected for the spurious dispersion induced by each system's S/N selection function (Appendix D).</figcaption></figure>

geometry<fnref n="4" />. Each catalog entry was populated by a (RA-Dec<sup>ran</sup>, W<sub>0</sub><sup>ran</sup>, v<sup>ran</sup>)-triplet, where (1) coordinates are uniformly distributed within the sky patch defined by all pre-selected

spaxels; (2) $W_0^{\text{ran}}$ is drawn randomly from the completeness-corrected $W_0$ distribution, with replacement – we did not use a model of the $W_0$ distribution (e.g., Mintz et al. 2022) because its shape is unknown in arc data; and (3) $v^{\text{ran}}$ is drawn randomly from a uniform distribution of velocities, $U_{[-4000,4000]}$, of spatially unresolved velocity components in the search window, which we consider the sample boundary (Mo et al. 1992). (See Appendix F.2 for more details on the third condition.) Entries were rejected if $W_0^{\text{ran}} < 3 \times \delta_W$, where $\delta_W = \text{FWHM}/\langle S/N \rangle / (1+z)$ is the detection limit set by the average S/N at the wavelength corresponding to $v^{\text{ran}}$. These entries become random non-detections. Otherwise, if accepted, the entries become random detections. In practice, the S/N is quite flat over the small wavelength range under consideration, implying that the final selection is determined mostly by the spaxel-to-spaxel S/N variations. Finally, RR pairs were built in exactly the same fashion as DD pairs. The fidelity of the random catalogs is tested in Appendix F.3 using mock catalogs of unclustered data.

## 6.4. Results

Figure 5 displays $\xi(\Delta r_{\perp}, \Delta v)$ computed for all eight systems. The arbitrary sampling in $\Delta r_{\perp}$ and $\Delta v$ minimizes the number of no-data bins while keeping enough S/N in the less populated bins. The signal indicates evident velocity clustering below ~100 km s⁻¹ in most systems. For the systems 1527 b and 0033 a, the velocity differences are more uniformly distributed. In general, no spatial correlation was observed. If present, either our “resolution element” (the total area per system) does not resolve it or the de-lensed coordinates are too uncertain over these scales.

Figure 6 shows the velocity projection of $\xi(\Delta r_{\perp}, \Delta v)$, hereafter $\xi^{\text{arc}}(\Delta v)$, created through merging bins in the spatial direction. Errors include both the bootstrapping analysis and Poisson statistics but are dominated by the latter. Consistent with Fig. 5, a significant amount of power is seen below ~100 km s⁻¹. For some systems (1527 b, 0033 a and 0033 b), a similar amount of power is seen for velocity differences up to 200 km s⁻¹. We note that none of these signals can be mimicked

<footnote n="4">⁴ We did not bootstrap because replacement would bias the pairs toward zero separation.</footnote>
<figure><img src="image_6.png" /><figcaption>Fig. 5. C IV transverse velocity autocorrelation, $\xi(\Delta r_{\perp}, \Delta v)$. Gray colored bins indicate no data.</figcaption></figure>

<caption>Fig. 6. Velocity projection of $\xi(\Delta r_{\perp}, \Delta v)$. Blue squares represent the arc data. The 1$\sigma$ errors (in pink) include both the bootstrapping analysis and the Poisson statistics described in Sect. 6.3 but are dominated by the latter. Dashed curves indicate the signal produced by a velocity field distributed as $N(0, \sigma_{\nu}^2)$, where $\delta_{\nu}$ are the individual measurement errors. Green curves show the model prediction (Sect. 7.2.4).</caption>

by measurement errors only (dashed lines), except perhaps for 1527 c. Remarkably, $\xi^{\text{arc}}(\Delta v)$ exhibits a strong similarity – in shape and amplitude – to the velocity two-point correlation function measured along quasar sight lines (hereafter $\xi^{\text{QSO}}(\Delta v)$). This seems sensible because $\xi^{\text{QSO}}(\Delta v)$ is measured for C IV “velocity components” in high-resolution (HR; $R \sim 40\,000$) spectra (Rauch et al. 1996; Scannapieco et al. 2006; Fathivavsari et al. 2013; Boksenberg & Sargent 2015) that our data cannot resolve. This similarity is this article’s main object of study.

## 7. Kinematic model

In this section we present a simple kinematic model that attempts to explain the results reported in Sect. 5 and Sect. 6. Subsequently, we show that this model provides testable clues on the kiloparsec-scale kinematic and spatial structure of the C IV gas.

### 7.1. Model setup

We propose an underlying population of sub-spaxel absorbing gas clumps (or “clouds”) that the ARCTOMO observations cannot resolve spatially but that nevertheless induce detectable absorption at the cloud velocity. In line with the lack of strong spatial structure seen in $\xi(\Delta r_{\perp}, \Delta v)$, we consider that the clouds are uniformly distributed on the sky and follow a normal distribution of projected velocities $N(0, \sigma_0^2)$ on the spatial scales probed here. To statistically re-create the observations, we envi-

sioned spaxels that sample $N$ such clouds on average (Fig. 7). We did not make any assumption regarding the intrinsic position of the clouds along the line of sight.

With this simplification, we computed the statistics of the sampled data, namely, (1) the distribution of the spaxel mean velocities and their dispersion, $\sigma_{\perp}$, and (2) the distribution of the spaxel velocity dispersions, $\sigma_{\parallel}$. Examples of (1) and (2) using $\sigma_0 = 100 \text{ km s}^{-1}$ and two different values of $N$ are shown in Fig. 7. Computing $\sigma_{\perp}$ is exactly the same as measuring the velocity dispersion in an ARCTOMO system. Computing $\sigma_{\parallel}$, on the other hand, is analogous to measuring the absorption dispersion in an ARCTOMO spectrum (but not exactly the same, as we are not dealing with absorption profiles). These definitions are based on the assumption that $\sigma_{\parallel}$ conveys the kinematic information of individual velocity components in the real observations (i.e., “within” a spaxel). Hence, panels (a) and (b) in Fig. 7 are the model analog to the observations displayed in columns 4 and 5 of Fig. 1, respectively.

The histograms in Fig. 7 show the distributions of $v$ and $\sigma_{\parallel}$. They resemble the data histograms in Fig. 3. Indeed, a suit of model realizations shows that the larger the $N$, the narrower and more separated the two distributions become. For velocities (blue histograms), this is straightforward to see as the well-known consequence of spatial re-sampling, a process which ‘blurs’ the signal. In our case,

$$ \sigma_{\perp} = \sigma_0 \times \frac{1}{\sqrt{N}}. \qquad (2) $$
<figure><img src="image_7.png" /><figcaption>Fig. 7. Kinematic model with N = 4 (upper row) and N = 25 clouds sampled per spaxel. In both cases the underlying velocity field dis- tributes as N(0, 100²). Left-hand column: Sampled velocities v (i.e., mean of N cloud velocities in the spaxel). Their standard deviation is the “transverse dispersion,” σ⊥. Middle column: Sampled standard devia- tions of the N cloud velocities, or “parallel” dispersion, σ‖. Right-hand column: Sample distribution of v and σ‖. The Horizontal blue line rep- resents the same as in Fig. 3.</figcaption></figure>

On the other hand, for σ_‖ (orange histograms), a large N not only narrows the distribution but it also shifts their peak to larger values. This effect is because, as more clouds are sampled per spaxel, their line-of-sight dispersions approach the original values of σ₀. Thus, there must exist a relation also between the original dispersion, σ₀, and the median line-of-sight dispersion, ⟨σ_‖⟩. An analytic solution (Kenney & Keeping 1951) exists for the mean only, not the median of the σ_‖ distribution. A good (within 1%) approximation for the median was found to be (Appendix E):

$$ \langle \sigma_{\parallel} \rangle \approx \sigma_0 \times \frac{N-1}{N}. \qquad (3) $$

As expected, the larger the number of clouds in a spaxel, the more the recovered dispersion approaches the original one. We note that the two equations above hold true independent of the physical spaxel size and depend only on N and σ₀.

## 7.2. Model predictions

### 7.2.1. Transverse versus parallel velocity dispersion

Equations (2) and (3) show that the spatial sampling determines σ⊥ and ⟨σ_‖⟩ and that this is a purely observational effect. Since we can measure σ⊥ and ⟨σ_‖⟩ in ARCTOMO data, the equations can be used to derive two intrinsic properties of the cloud distribution, N and σ₀. Solving for σ₀, one obtains:

$$ \langle \sigma_{\parallel} \rangle = \sigma_{\perp} \times \frac{N-1}{\sqrt{N}}, \qquad (4) $$

from which we derive N, the per-spaxel mean number of clouds. In Table 4, N is listed for each system. Equation (4) can also be used to re-create Fig. 4 for particular values of N. This is shown in Fig. 8, where the white dashed lines display the ⟨σ_‖⟩ versus σ⊥ relation for selected values of N. Since N does not need to be an integer, we interpreted it as the average number of clouds per spaxel. Equation (4) also explains why, under the assumption of Gaussianity, one should generally measure ⟨σ_‖⟩ > σ⊥ as a

<figure><img src="image_8.png" /><figcaption>Fig. 8. Same as Fig. 4 but showing the two predictions of our kinematic model on (a) ⟨σ_‖⟩ versus σ⊥ for a given number of clouds per spaxel, N (white dashed lines; Eq. (4)), and (b) the underlying line-of-sight dispersion, σ₀ (background color; Eq. (5)).</figcaption></figure>

result of the particular ARCTOMO observations “flattening” the underlying velocity field<fnref n="5" />.

Likewise, arranging Eqs. (2) and (3), this time to get rid of N, one finds

$$ \sigma_0 = \frac{\langle \sigma_{\parallel} \rangle \pm \sqrt{\langle \sigma_{\parallel} \rangle^2 + 4\sigma_{\perp}^2}}{2}. \qquad (5) $$

Finally, an equation emerges that gives the underlying C IV velocity dispersion out of two ARCTOMO observables. The solutions for each system are listed in Table 4, and the background of Fig. 8 displays σ₀ computed for all (σ⊥, ⟨σ_‖⟩) combinations, as indicated by the color scale. The uncertainties in Table 4 were propagated from those in σ_‖ and σ⊥. As a safety check for Eq. (5), if σ⊥ = 0 (no transverse dispersion), then σ₀ = ⟨σ_‖⟩; that is, the underlying dispersion is recovered (and all spectra show the same centroid velocity). The second solution, σ₀ = 0 (no absorption), is physical only for σ_‖ = 0.

### 7.2.2. Inter-cloud distances

So far, our model does not involve physical scales, but we can define the (projected) mean number of clouds per unit area as ⟨N_c⟩ ≡ N/A_spxel, where A_spxel is the mean spaxel area per system (hence, ⟨N_c⟩ is equivalent to the “counts-in-cylinder” statistics on galaxy scales; Berrier et al. 2011). With this definition, the projected mean distance between clouds is

$$ \langle d_c^{2D} \rangle \sim \left( \frac{1}{\langle N_c \rangle} \right)^{1/2} [\text{kpc}]. \qquad (6) $$

Using the values for ⟨N_c⟩ listed in Table 4, the present systems occur in structures separated on the sky by ⟨d_c^2D⟩ = 0.3–0.8 kpc. We discuss the implications of ⟨d_c^2D⟩ on cloud “sizes” in Sect. 7.2.3.

<footnote n="5">⁵ Interestingly, ⟨σ_‖⟩ = σ⊥ for N = 1 + φ, where φ is the “Golden ratio.”</footnote>
<table><caption>**Table 4.** Model predictions.</caption><thead><tr><th>System</th><th><i>N</i></th><th>σ<sub>0</sub><br>[km/s]</th><th>⟨<i>N</i><sub>c</sub>⟩<br>[kpc<sup>-2</sup>]</th><th><i>S</i><sub>c</sub><br>[kpc]</th></tr><tr><th></th><th>(1)</th><th>(2)</th><th>(3)</th><th>(4)</th></tr></thead><tbody><tr><td>1527 a</td><td>12.4±3.4</td><td>143.0±4.9</td><td>3.4±1.1</td><td>0.5–1.9</td></tr><tr><td>1527 b</td><td>5.9±1.7</td><td>159.8±8.9</td><td>1.6±0.5</td><td>0.8–1.9</td></tr><tr><td>1527 c</td><td>9.5±4.7</td><td>69.0±5.6</td><td>3.3±1.8</td><td>0.5–1.7</td></tr><tr><td>G311 a</td><td>6.0±1.4</td><td>74.4±3.7</td><td>12.0±3.5</td><td>0.3–0.7</td></tr><tr><td>0033 a</td><td>5.1±1.2</td><td>167.9±7.6</td><td>10.6±3.0</td><td>0.3–0.7</td></tr><tr><td>0033 b</td><td>6.7±1.8</td><td>149.6±7.4</td><td>13.0±4.2</td><td>0.3–0.7</td></tr><tr><td>2111 a</td><td>3.5±1.5</td><td>121.2±11.1</td><td>7.1±3.3</td><td>0.4–0.7</td></tr><tr><td>2111 b</td><td>4.1±1.5</td><td>110.2±8.5</td><td>6.7±2.8</td><td>0.4–0.8</td></tr></tbody></table>

**Notes.** Table columns: (1) number of clouds per spaxel, (2) underlying velocity dispersion, (3) mean number of clouds per unit area, and (4) cloud sizes for covering fraction $κ = 1$ (Sect. 7.2.3). Errors are propagated from $σ_⊥$ and $σ_‖$ uncertainties.

Likewise, $⟨N_c⟩$ could constrain the number density of clouds $n_c = N_c/L$, where $L$ [kpc] is the (unknown) total absorption length. In this case, the mean 3D distance between clouds is

$$⟨d_c^{3D}⟩ \sim \frac{1}{n_c^{1/3}} \sim \left(\frac{L}{⟨N_c⟩}\right)^{1/3} \quad [\text{kpc}]. \qquad (7)$$

From Table 4, we observed that the present systems arise in structures separated in space by $⟨d_c^{3D}⟩ = 2-4$ (L/100)$^{1/3}$ kpc. Current estimates on the extension of $W_0 \ge 0.3$ Å C IV halos amount to $\approx 100-200$ kpc (Steidel et al. 2010; Hasan et al. 2022).

### 7.2.3. Cloud sizes and covering fraction

While our model does not constrain cloud “sizes” directly, our results suggest that cloud sizes cannot exceed the arc beams. In fact, if the clouds were much larger than the beam, we would measure (within observational uncertainties) the same velocity centroids across the entire arc and therefore also $σ_⊥ \approx 0$. As a result, $ξ^{arc}(Δv)$ would be much steeper (given that $Δv \approx 0$ across spaxels) than what we measure (Fig. 6). The only way to produce the observed signals at larger $Δv$ is by having a cloud size that is about the same size or smaller than the spatial scale probed by a spaxel.

On the other hand, our model predicts a projected inter-cloud distance, $⟨d_c^{2D}⟩$ (Eq. (6)). This parameter is related to the cloud characteristic size, $S_c$, and to the covering fraction, $κ$, defined as the fraction of the beam area covered by all clouds in projection along the line of sight. Our data cannot directly disentangle $S_c$ and $κ$, but the following scenarios seem plausible:

1. $κ \approx 1$, at the limit of no overlap: This implies cloud sizes of $S_c \approx ⟨d_c^{2D}⟩$.

2. $κ \approx 1$, with considerable overlap: This implies $S_c \ge ⟨d_c^{2D}⟩$ but below the beam size for the reasons outlined above.

3. $κ < 1$: In this case, sizes remain unconstrained but below $⟨d_c^{2D}⟩$.

Table 4 displays $S_c = ⟨d_c^{2D}⟩$ (first scenario, which we consider the most likely), using Eq. (6) and the tabulated $⟨N_c⟩$ values. We emphasize that these considerations apply to clouds producing detectable ARCTOMO signals (i.e., typically $W_0 \ge 0.3$ Å).

In conclusion, our kinematic model constrains strong C IV cloud sizes to be of the order of or smaller than the probed

beam size (i.e., $\le 1$ kpc). Future higher spatial and spectral resolution data should help discern between the scenarios described above.

### 7.2.4. Transverse auto-correlation

Our kinematic model can also predict $ξ^{arc}(Δv)$ (Sect. 6.4). Recalling that if a random variable $X$ is normally distributed with variance $σ^2$, then $ΔX$ is normally distributed, too (with variance $2σ^2$), Eq. (2) implies that

$$\sigma_{\xi}^{\text{arc}} = \sigma_{0} \times \sqrt{\frac{2}{N}}, \qquad (8)$$

where $σ_{\xi}^{\text{arc}}$ is the $1σ$ width of $ξ^{arc}(Δv)$. The green curves in Fig. 6 compare this prediction with the measured $ξ^{arc}(Δv)$. Since the model only predicts the width of $ξ^{arc}(Δv)$, we display a Gaussian with an amplitude given by the maximum value of $ξ^{arc}(Δv)$ and a width given by Eq. (8), using $σ_0$ and $N$ listed in Table 4 (a comparison between fitted and predicted widths is shown in Fig. F.2). The measurement and the prediction match reasonably well, which demonstrates the self-consistency of the model (since both quantities are based on the same data).

### 7.2.5. Prediction using the quasar line-of-sight correlation

As an independent test, our model can also predict $ξ^{arc}(Δv)$ from $ξ^{QSO}(Δv)$. To this end, we first stacked and averaged $ξ^{arc}(Δv)$ over all eight systems. Next, we assumed that the C IV clouds abide to a normal distribution that is captured by both $ξ^{arc}(Δv)$ and $ξ^{QSO}(Δv)$; in other words, the clouds are responsible for both the parallel signal toward quasars and the transverse signal toward arcs. This is a strong assumption that we discuss further below in Sect. 8.1. From Eq. (8), it follows that $σ_{\xi}^{\text{arc}} = σ_{\xi}^{\text{quasar}} / \sqrt{N}$. Boksenberg & Sargent (2015) fit $ξ^{QSO}(Δv)$ with the sum of two Gaussians having $σ_{\xi}^{\text{quasar}} = 80$ and $185$ km s$^{-1}$ (also with a narrower one that our data do not resolve). Scannapieco et al. (2006) and Fathivavsari et al. (2013) do not provide fits, but their correlation functions have widths consistent with these. In Fig. 9, we display a comparison between our transverse correlation and the line-of-sight “QSO prediction” using $N = 6.7$, the weighted average number of clouds per spaxel in the arc data. The match is quite good ($\le 2σ$ deviation), lending independent support to our kinematic model and to the inferred number of clouds per spaxel. It also suggests that the present (modest) number of C IV systems is only moderately affected by cosmic variance.

## 8. Discussion

The proposed kinematic model seems acceptable on two different fronts: a self-consistent comparison between ARCTOMO transverse and line-of-sight kinematics, and a prediction of the transverse velocity auto-correlation from the quasar line-of-sight kinematics. The first one is solely observational, as we have outlined above. But the second suggests that the weak C IV components (only detected in HR quasar spectra) and the strong components (our arc detections) may trace the same physical structures. We conclude this article by discussing how our results are consistent with even more quasar observations and how together they provide glimpses of the CGM on both the sub-kiloparsec and galactic scales.
<figure><img src="image_9.png" /><figcaption>**Fig. 9.** Stacked average $\xi^{\text{arc}}(\Delta v)$ of all eight arc systems and QSO prediction (details in Sect. 7.2.5). The dashed curves indicate the range of arc systems.</figcaption></figure>

## 8.1. Cloud density structure

A natural explanation for a single cloud population giving rise to both strong and weak components is that each cloud has an internal density structure. Consequently, as a second-order approximation of the model, in the following we refer to cloud “cores” (the central parts of the clouds responsible for the strong components) and “halos” (external parts of the clouds responsible for the weak components).

Figure 10 illustrates the proposed situation. The left-hand panel displays a column-density “map” of cloud cores and halos. The white dot represents a narrow beam (e.g., a quasar beam), and the squares represent the arc beams (e.g., ARCTOMO spaxels). In this simple example, the beams pierce a clumpy medium represented by the superposition of just four clouds. The right-hand panels display the respective spectra. These result from spatially averaging the flux within the beam areas (see the figure caption for more details on the implementation). By model construction, the narrow beam collects the kinematic information of many clouds through intercepting mostly their larger cross-section cloud halos. Each cloud produces a velocity component in an HR spectrum (lower-right panel). Conversely, a LR spectrum of the wide beam recovers only an average velocity and a line-of-sight dispersion ($\sigma_{\|}$), and it is dominated by a few cloud cores inside the beam (upper-right panels). As expected, $\langle\sigma_{\|}\rangle > \sigma_{\perp}$ (see the figure caption for how $\sigma_{\|}$ and $\sigma_{\perp}$ are represented).

A scenario of clouds producing both weak and strong components also explains why removing the high column-density components from quasar HR data does not affect $\xi^{\text{QSO}}(\Delta v)$ (Scannapieco et al. 2006). Those authors find little change in $\xi^{\text{QSO}}(\Delta v)$ when $N = 10^{13,14,15} \text{ cm}^{-2}$ components are excluded (up to 10% of their sample, according to the reported column density frequency). In the present interpretation, this effect is not simply due to the strong components being a small fraction of the total but more fundamentally to both strong and weak components sharing the same velocity field.

A cloud density structure equivalent to our model has already been proposed. For instance, Hummels et al. (2024) have introduced the concept of a cool CGM complex, which is composed of multiple cloudlets of various masses and velocities that lead to a column density structure, as in our case. However, we emphasize that our model aims to explain the observed kinematics only. Besides, the canonical one-to-one association between the cloud and the velocity component may also be too simplistic (e.g., Faerman & Werk 2023; Marra et al. 2024; Li et al. 2024).

## 8.2. Other multiple sight line observations

The most direct comparison between arc and quasar results can be made with the – unfortunately very limited – sample of multi-sight line quasar observations subtending approximately kiloparsec separations. The first direct estimate of CIV coherence length using lensed quasars (Rauch et al. 2001b) delivered ~0.3 kpc (50% transverse variations). This value is consistent with subsequent studies that found CIV transverse structure below 1 kpc (Tzanavaris & Carswell 2003; Lopez et al. 2007; Rubin et al. 2018b). Along with the ten times larger $\langle d_c^{3D} \rangle$ computed in Sect. 7.2.2, the above narrow-beam size constraints imply a small volume filling factor of ~$10^{-3}$, which is in turn comparable with theoretical predictions (McCourt et al. 2018; Gronke & Oh 2020; Liang & Remming 2020; Li et al. 2024). On the other hand, combining those sizes with our $\langle d_c^{2D} \rangle$ constraint (Eq. (6)) leads to a near-unity covering fraction. This is again consistent with theoretical predictions (Liang & Remming 2020). We conclude that our arc results are consistent with quasar observations and theoretical predictions of the CIV spatial domain.

Regarding the line-of-sight direction, we measured $\langle\sigma_{\|}\rangle \approx 60-130 \text{ km s}^{-1}$, a range which seems consistent with the velocity shear found toward lensed quasars. Rubin et al. (2018b) find a velocity structure $\Delta v \approx 100 \text{ km/s}$ over kiloparsec scales; Rauch et al. (2001b) found $\approx 60 \text{ km s}^{-1}$ over 10 kpc; Ellison et al. (2004) measured $\approx 15 \text{ km s}^{-1}$ on sub-kiloparsec scales; and Lopez et al. (2007) found some cases with $60 \text{ km s}^{-1}$ shear over 1 kpc. Thus, also in the parallel direction, arc and quasar observations seem consistent with each other.

## 8.3. Departures from Gaussianity and large-scale kinematics

Our model assumes an underlying Gaussian velocity field. In order to gauge departures from this idealization, we considered a spatial gradient in velocity, $\Delta v$, produced, for instance, by galaxy-scale motions such as orbiting, co-rotating, or outflowing gas. Assuming the spatial scale of $\Delta v$ is much larger than the spaxel size, it is straightforward to see that $\sigma_{\perp}$ will increase but $\sigma_{\|}$ will not. This is because the former is an inter-beam measurement, whereas the latter is an “intra-beam” one. In this case, Eq. (2) does not hold anymore, and in Fig. 4 all points become shifted to the right. A simple numerical test shows that $N$ becomes underpredicted as soon as $\Delta v$ is comparable to $\sigma_0$. Thus, departures from Gaussianity could be detected in ARCTOMO systems exhibiting $\sigma_{\perp} \gtrsim \langle\sigma_{\|}\rangle$ and $\xi(\Delta r_{\perp}, \Delta v)$ power both at large spatial and velocity scales. None of the present systems clearly qualifies in this group, but in principle our experiment could be used to identify large-scale kinematic motions of the absorbing CGM.

On the other hand, the large underlying velocity dispersions we derived are intriguing. The largest $\sigma_0$ values in Table 4 suggest velocity dispersions of galaxy groups, and one possibility is that the unresolved clumps are bound to different group
<figure><img src="image_10.png" /><figcaption>**Fig. 10.** Semi-schematic visualization of the kinematic model. Left-hand panel: Column density map, *N*(*x*). Here, *x* denotes the on-sky position. The map was created by superposing four clouds with *N*(*C IV*) = 10¹⁴ cm⁻² (cores) and 10¹² cm⁻² (halos) at random positions and line-of-sight velocities, *v*<sub>c</sub>. The clouds are colored according to *v*<sub>c</sub>. The white dot represents a narrow beam (e.g., a quasar pencil beam), and the four squares indicate wide beams (e.g., ARCTOMO spaxels). An optical-depth cube τ(*v*, *x*) was created by evaluating Voigt profiles (Liang & Kravtsov 2017) with *v* = *v*<sub>c</sub>, *N* = *N*(*x*), and *b* = 10 km s⁻¹, and spectra were generated by averaging e⁻⁻τ within the beam area. Upper right panels: Wide-beam spectra at infinite and MUSE spectral resolution (HR and LR, respectively). The green squared symbols with error bars indicate LR individual line velocity centroids and σ<sub>||</sub> dispersions, respectively, obtained from Gaussian fits. The shaded region indicates the overall σ<sub>⊥</sub> dispersion of the four LR velocities. Lower right panel: Narrow-beam spectrum at infinite spectral resolution. We note that *v*<sub>c</sub> is indicated with the same colors as in the left-hand panel.</figcaption></figure>

galaxies. Another one, of course, is that they feature the long-
sought manifestation of super winds (Voit 1996; Pettini et al.
2001). Associations between C IV in quasar spectra and star-
forming galaxies find kinematic separations at this level, and
these have been associated with clustering (Adelberger et al.
2005; Lofthouse et al. 2023; Banerjee et al. 2023; Galbiati et al.
2023), inflows (Turner et al. 2017), or filaments outside the
galaxy virial radius (Galbiati et al. 2023; Banerjee et al. 2023).
All of these tests suggest that C IV in quasar spectra trace
galaxies, although none of them really can disentangle the
absorbing galaxy from the absorption system. Our transverse
observations literally add a new dimension to the under-
standing of the origin of metal-rich gas in the z = 2–3
cool-warm CGM.

8.4. *The path forward*

It may seem surprising that the present ARCTOMO data,
although unable to resolve the absorbing clouds neither spatially
nor spectroscopically, still carry the 3D kinematic information
encoded in their absorption profiles. We have shown that these
data, along with simple assumptions about the kinematics, can
lead to a set of realistic and testable predictions.

But the current tomographic data do not yet allow for a com-
plete disentangling between different global dynamic models
of the C IV-bearing CGM. The definitive pieces of the puzzle
will be obtained through (1) the detection of the C IV galax-
ies responsible for the arc signals, a challenging objective that
requires space-based near-infrared observations, and (2) spa-

tially resolved HR observations of extended background sources,
where little has been done yet (e.g., Diamond-Stanic et al. 2016).
Along with more sophisticated models that include density struc-
ture, we assert that such observations shall dramatically improve
our understanding of the cool-warm CGM’s small-scale struc-
ture. We consider, for instance, the HR curves in the upper-right
panels of Fig. 10. If the kinematic model tested in this article
is reliable, a small number of clouds per spatial resolution ele-
ment of approximately kiloparsec size should unfold as resolved
velocity components.

9. Summary

We have presented spatially resolved VLT/MUSE observations
of four giant gravitational arcs that offer C IV coverage. We
focused on the kinematic properties of C IV absorption detected
at 〈z<sub>abs</sub>〉 ~ 2.3, a redshift that enables sub-kiloparsec resolution
of the absorbers thanks to lens magnification. Our experimen-
tal setup allowed us to analyze C IV absorption in 222 adjacent,
uncorrelated beams that pierce eight intervening C IV systems.
Our results are as follows:

1. Significant absorption is detected across almost all the arcs, implying C IV extensions (W₀ ≥ 0.3 Å) of at least ≈10 kpc in the reconstructed absorber plane (Fig. 1). We ran an automated search of C IV absorption and computed doublet velocities, absorption spreads, and equivalent widths. Each system shows evident velocity clustering in the transverse direction (Fig. 3), and we set out to measure such clustering and investigated its origin.
2. On average, the transverse velocity dispersion, $\sigma_{\perp}$, is found to be smaller than the per-system median line-of-sight dispersion, $\langle\sigma_{\parallel}\rangle$ (Fig. 4).

3. To measure spatial clustering, we computed a transverse auto-correlation function of CIV velocities ($\xi(\Delta r_{\perp}, \Delta v)$; Fig. 5) and its velocity projection ($\xi^{\text{arc}}(\Delta v)$; Fig. 6). $\xi(\Delta r_{\perp}, \Delta v)$ does not show evident spatial patterns, perhaps due to insufficient resolution. On the other hand, the average $\xi^{\text{arc}}(\Delta v)$ exhibits great similarity in shape and amplitude with the line-of-sight two-point correlation measured in HR quasar spectra, $\xi^{\text{QSO}}(\Delta v)$.

4. To aid a comparison between wide- (e.g., arc) and narrow- (e.g., quasar) beam observations, we introduced a simple kinematic model in which the absorption profiles result from groups of *N* clumps (clouds) sampled at the sub-spelax level and that produce a mean velocity and line-of-sight dispersion (Fig. 7). The model successfully explains the $\sigma_{\perp} < \sigma_{\parallel}$ inequality under the assumption of an underlying Gaussian field with dispersion $\sigma_0$. It also consistently predicts *N* (Eq. (4)) and $\sigma_0$ (Eq. (5)) out of the observables $\sigma_{\perp}$ and $\sigma_{\parallel}$ (Fig. 8). Combining this information with the average spaxel areas allowed us to constrain the number of clouds per unit area to within $N_c \approx 2-13 \text{ kpc}^{-2}$ and $\sigma_0 \approx 70-170 \text{ km s}^{-1}$, depending on the system. The model also constrains the projected mean inter-cloud distance to within $\approx 0.3-0.8 \text{ kpc}$. The covering fraction and kinematic considerations put strong CIV sizes between those values and the spaxel size, 0.7–1.9 kpc.

5. Our model also predicts the width of $\xi^{\text{arc}}(\Delta v)$ out of the independently computed $\xi^{\text{QSO}}(\Delta v)$ reasonably well. Since $\xi^{\text{QSO}}(\Delta v)$ is dominated by the bulk of weak velocity components that the arc data cannot resolve, this match led us to conclude that a single population of CIV clouds must be responsible for both signals. This in turn implies that the clouds must have a density structure, for instance, a radial density profile that produces a few strong components that shape $\xi^{\text{arc}}(\Delta v)$ and many weak components that shape $\xi^{\text{QSO}}(\Delta v)$ and only get detected in HR spectra (e.g., of quasars; Fig. 10). In such a scenario, $\xi^{\text{QSO}}(\Delta v)$ is dominated by the weak components, as has been noted elsewhere.

6. We have discussed how our model and observations are compatible with extant observations of multiple quasars, both in the transverse and the parallel dimensions.

Upcoming large optical facilities will routinely enable tomography of the CGM using not only lensed galaxies but also normal galaxies. Comparison of these wide-beam observations with the extant narrow-beam statistics, as shown here for a handful of systems, promises a better understanding of the connection between the small- and the galactic-scale structure of the high-redshift CGM.

## Data availability

Appendix A displaying the CIV spectra, corresponding line-profile fits, and best-fit parameters can be found in the Zenodo repository at https://zenodo.org/records/13891488

Acknowledgements. We would like to thank the anonymous referee for a thorough and critical review of our analysis. This work has benefited from discussions with Max Gronke, Joseph Hennawi, and Claudio Lopez. S. L. and N. T. acknowledge support by FONDECYT grant 1231187. M. S. was financially supported by Becas-ANID scholarship #21221511, and also acknowledges ANID BASAL project FB210003.

## References

Adelberger, K. L., Shapley, A. E., Steidel, C. C., et al. 2005, ApJ, 629, 636

Afruni, A., Lopez, S., Anshul, P., et al. 2023, A&A, 680, A112

Bacon, R., Accardo, M., Adjali, L., et al. 2010, SPIE Conf. Ser., 7735, 773508

Banerjee, E., Muzahid, S., Schaye, J., Johnson, S. D., & Cantalupo, S. 2023, MNRAS, 524, 5148

Becker, G. D., Pettini, M., Rafelski, M., et al. 2019, ApJ, 883, 163

Berrier, H. D., Barton, E. J., Berrier, J. C., et al. 2011, ApJ, 726, 1

Bird, S., Rubin, K. H. R., Suresh, J., & Hernquist, L. 2016, MNRAS, 462, 307

Boksenberg, A., & Sargent, W. L. W. 2015, ApJS, 218,

Bordoloi, R., Tumlinson, J., Werk, J. K., et al. 2014, ApJ, 796, 136

Bordoloi, R., O'Meara, J. M., Sharon, K., et al. 2022, Nature, 606, 59

Bordoloi, R., Simcoe, R. A., Matthee, J., et al. 2024, ApJ, 963, 28

Cen, R., & Chisari, N. E. 2011, ApJ, 731, 11

Chen, H.-W., Lanzetta, K. M., & Webb, J. K. 2001, ApJ, 556, 158

Chen, H.-W., Gauthier, J.-R., Sharon, K., et al. 2014, MNRAS, 438, 1435

Christensen, C. R., Davé, R., Governato, F., et al. 2016, ApJ, 824, 57

Cooksey, K. L., Kao, M. M., Simcoe, R. A., O'Meara, J. M., & Prochaska, J. X. 2013, ApJ, 763, 37

Cooper, T. J., Simcoe, R. A., Cooksey, K. L., et al. 2019, ApJ, 882, 77

Coppolani, F., Petitjean, P., Stoehr, F., et al. 2006, MNRAS, 370, 1804

Cowie, L. L., Songaila, A., Kim, T.-S., & Hu, E. M. 1995, AJ, 109, 1522

Diamond-Stanic, A. M., Coil, A. L., Moustakas, J., et al. 2016, ApJ, 824, 24

D'Odorico, V., Finlator, K., Cristiani, S., et al. 2022, MNRAS, 512, 2389

D'Odorico, V., Bañados, E., Becker, G. D., et al. 2023, MNRAS, 523, 1399

Dutta, R., Acebron, A., Fumagalli, M., et al. 2024, MNRAS, 528, 1895

Ellison, S. L., Songaila, A., Schaye, J., & Pettini, M. 2000, AJ, 120, 1175

Ellison, S. L., Ibata, R., Pettini, M., et al. 2004, A&A, 414, 79

ESO CPL Development Team 2015, Astrophysics Source Code Library [record ascl:1504.003]

Faerman, Y., & Werk, J. K. 2023, ApJ, 956, 92

Fathivavsari, H., Petitjean, P., Ledoux, C., et al. 2013, MNRAS, 435, 1727

Faucher-Giguère, C.-A., & Oh, S. P. 2023, ARA&A, 61, 131

Faucher-Giguère, C.-A., Kereš, D., & Ma, C.-P. 2011, MNRAS, 417, 2982

Fernandez-Figueroa, A., Lopez, S., Tejos, N., et al. 2022, MNRAS, 517, 2214

Fielding, D., Quataert, E., McCourt, M., & Thompson, T. A. 2017, MNRAS, 466, 3810

Finlator, K., Doughty, C., Cai, Z., & Díaz, G. 2020, MNRAS, 493, 3223

Fischer, T. C., Rigby, J. R., Mahler, G., et al. 2019, ApJ, 875, 102

Galbiati, M., Fumagalli, M., Fossati, M., et al. 2023, MNRAS, 524, 3474

Gontcho A Gontcho, S., Miralda-Escudé, J., Font-Ribera, A., et al. 2018, MNRAS, 480, 610

Gronke, M., & Oh, S. P. 2020, MNRAS, 494, L27

Grossman, S. A., & Narayan, R. 1988, ApJ, 324, L37

Hamilton, A. J. S., & Tegmark, M. 2004, MNRAS, 349, 115

Hasan, F., Churchill, C. W., Stemock, B., et al. 2022, ApJ, 924, 12

Hennawi, J. F., & Prochaska, J. X. 2007, ApJ, 655, 735

Hennawi, J. F., Strauss, M. A., Oguri, M., et al. 2006, AJ, 131, 1

Hummels, C. B., Smith, B. D., Hopkins, P. F., et al. 2019, ApJ, 882, 156

Hummels, C. B., Rubin, K. H. R., Schneider, E. E., & Fielding, D. B. 2024, ApJ, 972, 148

Jullo, E., Kneib, J.-P., Limousin, M., et al. 2007, New J. Phys., 9, 447

Kenney, J. F., & Keeping, E. S. 1951, in Mathematics of statistics, ed. D. Van Nostrand (Princeton, N.J.)

Kereš, D., Katz, N., Weinberg, D. H., & Davé, R. 2005, MNRAS, 363, 2

Kerscher, M., Szapudi, I., & Szalay, A. S. 2000, ApJ, 535, L13

Krogager, J. K., Noterdaeme, P., O'Meara, J. M., et al. 2018, A&A, 619, A142

Kulkarni, V. P., Cashman, F. H., Lopez, S., et al. 2019, ApJ, 886, 83

Ledoux, C., Noterdaeme, P., Petitjean, P., & Srianand, R. 2015, A&A, 580, A8

Li, Z., Gronke, M., & Steidel, C. C. 2024, MNRAS, 529, 444

Liang, C., & Kravtsov, A. 2017, ArXiv e-prints [arXiv:1710.09852]

Liang, C. J., & Remming, I. 2020, MNRAS, 491, 5056

Lofthouse, E. K., Fumagalli, M., Fossati, M., et al. 2023, MNRAS, 518, 305

Lopez, S., Reimers, D., Gregg, M. D., et al. 2005, ApJ, 626, 767

Lopez, S., Ellison, S., D'Odorico, S., & Kim, T.-S. 2007, A&A, 469, 61

Lopez, S., Tejos, N., Ledoux, C., et al. 2018, Nature, 554, 493

Lopez, S., Tejos, N., Barrientos, L. F., et al. 2020, MNRAS, 491, 4442

Maitra, S., Srianand, R., Petitjean, P., et al. 2019, MNRAS, 490, 3633

Marra, R., Churchill, C. W., Kacprzak, G. G., et al. 2024, MNRAS, 527, 10522

Martin, C. L., Scannapieco, E., Ellison, S. L., et al. 2010, ApJ, 721, 174

McCourt, M., Oh, S. P., O'Leary, R., & Madigan, A.-M. 2018, MNRAS, 473, 5407

McQuinn, M. 2016, ARA&A, 54, 313

Mentz, J. J., La Barbera, F., Peletier, R. F., et al. 2016, MNRAS, 463, 2819

Mintz, A., Rafelski, M., Jorgenson, R. A., et al. 2022, AJ, 164, 51

Mo, H. J., Jing, Y. P., & Boerner, G. 1992, ApJ, 392, 452

Mortensen, K., Keerthi Vasan, G. C., Jones, T., et al. 2021, ApJ, 914, 92
Nelson, D., Genel, S., Vogelsberger, M., et al. 2015, MNRAS, 448, 59

Noterdaeme, P., Srianand, R., & Mohan, V. 2010, MNRAS, 403, 906

Oppenheimer, B. D., & Davé, R. 2006, MNRAS, 373, 1265

Oppenheimer, B. D., & Davé, R. 2008, MNRAS, 387, 577

Peeples, M. S., Corlies, L., Tumlinson, J., et al. 2019, ApJ, 873, 129

Péroux, C., Rahmani, H., Arrigoni Battaia, F., & Augustin, R. 2018, MNRAS,
479, L50

Petitjean, P., & Bergeron, J. 1994, A&A, 283, 759

Pettini, M., Shapley, A. E., Steidel, C. C., et al. 2001, ApJ, 554, 981

Pichon, C., Scannapieco, E., Aracil, B., et al. 2003, ApJ, 597, L97

Rahmati, A., Schaye, J., Crain, R. A., et al. 2016, MNRAS, 459, 310

Rauch, M., Sargent, W. L. W., Womble, D. S., & Barlow, T. A. 1996, ApJ, 467, L5

Rauch, M., Sargent, W. L. W., & Barlow, T. A. 1999, ApJ, 515, 500

Rauch, M., Sargent, W. L. W., & Barlow, T. A. 2001a, ApJ, 554, 823

Rauch, M., Sargent, W. L. W., Barlow, T. A., & Carswell, R. F. 2001b, ApJ, 562, 76

Rubin, K. H. R., Diamond-Stanic, A. M., Coil, A. L., Crighton, N. H. M., & Stewart, K. R. 2018a, ApJ, 868, 142

Rubin, K. H. R., O'Meara, J. M., Cooksey, K. L., et al. 2018b, ApJ, 859, 146

Rudie, G. C., Steidel, C. C., Pettini, M., et al. 2019, ApJ, 885, 61

Rupke, D. S., Veilleux, S., & Sanders, D. B. 2005, ApJS, 160, 115

Ryan-Weber, E. V., Pettini, M., & Madau, P. 2006, MNRAS, 371, L78

Savitzky, A., & Golay, M. J. E. 1964, Anal. Chem., 36, 1627

Scannapieco, E., Pichon, C., Aracil, B., et al. 2006, MNRAS, 365, 615

Schaye, J., Carswell, R. F., & Kim, T.-S. 2007, MNRAS, 379, 1169

Sharon, K., Bayliss, M. B., Dahle, H., et al. 2020, ApJS, 247, 12

Shen, S., Madau, P., Aguirre, A., et al. 2012, ApJ, 760, 50

Simcoe, R. A., Sargent, W. L. W., & Rauch, M. 2004, ApJ, 606, 92

Songaila, A. 2005, AJ, 130, 1996

Soto, K. T., Lilly, S. J., Bacon, R., Richard, J., & Conseil, S. 2016, MNRAS,
458, 3210

Steidel, C. C. 1990, ApJS, 72, 1

Steidel, C. C., Erb, D. K., Shapley, A. E., et al. 2010, ApJ, 717, 289

Stern, J., Hennawi, J. F., Prochaska, J. X., & Werk, J. K. 2016, ApJ, 830, 87

Tejos, N., Morris, S. L., Finn, C. W., et al. 2014, MNRAS, 437, 2017

Tejos, N., López, S., Ledoux, C., et al. 2021, MNRAS, 507, 663

Theuns, T., Viel, M., Kay, S., et al. 2002, ApJ, 578, L5

Turner, M. L., Schaye, J., Crain, R. A., et al. 2017, MNRAS, 471, 690

Tytler, D., Gleed, M., Melis, C., et al. 2009, MNRAS, 392, 1539

Tzanavaris, P., & Carswell, R. F. 2003, MNRAS, 340, 937

Vikas, S., Wood-Vasey, W. M., Lundgren, B., et al. 2013, ApJ, 768, 38

Voit, G. M. 1996, ApJ, 465, 548

Weilbacher, P. M., Palsa, R., Streicher, O., et al. 2020, A&A, 641, A28

Wiersma, R. P. C., Schaye, J., Dalla Vecchia, C., et al. 2010, MNRAS, 409, 132

Zahedy, F. S., Chen, H.-W., Rauch, M., Wilson, M. L., & Zabludoff, A. 2016, MNRAS, 458, 2423

Zahedy, F. S., Chen, H.-W., Johnson, S. D., et al. 2019, MNRAS, 484, 2257
## Appendix B: On system identification and completeness

The C IV system identification runs over all masked spectra along the whole available redshift path and provides redshift candidates that feed the subsequent line profile fitting (described in Sect. 3.2). The algorithm employs Pearson correlation. To begin, whole-spectrum continua are estimated through iterative Savitzky-Golay filtering (Savitzky & Golay 1964). The identification is performed along every continuum-normalized spectral pixel having S/N> 4, wherein the observed flux is correlated with a template. The template is built as two inverted Gaussian profiles defined by the characteristic separation of the C IV doublet and the instrumental spectral resolution. High correlation (Pearson coefficient $r > 0.7$, typically) indicates a possible true doublet (Noterdaeme et al. 2010; Ledoux et al. 2015). The candidates are visually examined and classified as a C IV system if high correlation occurs in at least one binned spectrum. Eight systems are found (Table 2); notably, their $W_0$ distributions reach such high values that these systems would have been easily discovered by eye. This suggests we are not missing systems of this kind. To roughly estimate the system completeness at the low-$W_0$ end, we inject synthetic C IV doublets with properties $z_i$ and $W_0$ and run the search in the vicinity of the i-pixel of each binned spectrum. We obtain on average ≈ 50% recovery wherein all spectra in a system have a doublet with $W_0 \sim 0.7$ Å.

<figure><img src="image_11.png" /><figcaption>Fig. B.1. Distribution of errors in centroid velocity ($\delta_v$) and absorption spread ($\delta_\sigma$).</figcaption></figure>

## Appendix C: Possible spatial sampling biases

To quantify spaxel-to-spaxel cross-talk effects, we performed two tests:

1. On the comparison between line-of-sight and transverse velocity dispersions (Sect. 5): we re-extracted and re-fit spectra using a larger aperture of 0.8 × 0.8-binned spaxels. Consistently, this decreases the number of spectra (by ~ 40%) but should also counteract any possible cross-talk effects due to a FWHM ≈ 0.8 PSF. The corresponding version of Fig. 4 shows broad consistency with the original one, that is, $\sigma_\perp < \sigma_\parallel$ holds. Therefore, this test indicates that the inequality is not an artifact induced by the chosen 0.6 × 0.6 aperture.

2. On the transverse auto-correlation of velocities (Sect. 6): a 40% decrement in detections precludes the analysis of some of the systems due to their small number of pairs. Instead, to quantify the possible effect of 0.6 × 0.6 apertures on $\xi(\Delta r_\perp, \Delta v)$ we recompute it excluding all neighbor spaxels in the image plane; that is, pairs made of information from adjacent binned spaxels are not considered. This exercise removes ~ 10% of the sample pairs. Figure C.1 shows the results (for a better comparison we keep the mock and

<figure><img src="image_12.png" /><figcaption>Fig. C.1. Same as Fig. 6 but excluding pairs by neighbor spaxels in the image plane.</figcaption></figure>

model curves as in Fig. 6). In all systems the correlation signal, although noisier, is preserved, indicating it is not driven by pairs of adjacent spaxels. Therefore, velocities measured using the 0.6 × 0.6 aperture can be considered independent, that is, not biased by spaxel cross-talk effects. The exception is perhaps system 1527 c, whose correlation function is close to the signal expected from measurement errors only (indicated by the dashed lines in both figures). However, this system also shows the strongest overlapping effect in the reconstructed absorber plane, so $\xi(\Delta r_\perp, \Delta v)$ could be biased toward $\Delta v = 0$ km s$^{-1}$ due to limitations of the lens model.

## Appendix D: Spurious transverse dispersion

In Sect. 5 we compare $\sigma_\perp$ and $\langle\sigma_\parallel\rangle$. Naturally, we expect that low S/N spectra will lead to an increased scatter in the velocity uncertainties measured from Gaussian fitting, resulting in wider velocity distributions and thus higher $\sigma_\perp$. In order to assess the magnitude of this systematic effect, we perform a MC re-sampling test. Each detected C IV doublet is ‘displaced’ to the same redshift, simulating a single-velocity distribution. This is achieved by replacing the originally fitted absorption lines with randomly selected flux from the adjacent pixels. After the doublet is essentially erased a double Gaussian profile is multiplied with the flux. These profiles have the same properties as the absorption features found in the data, for the corresponding spectra, except for the redshift which is set to the defined velocity center of the system. $\sigma_\perp$ obtained from fitting these doublets (Table 3) should be purely due to each system’s S/N selection function.

## Appendix E: Sample dispersion approximation

Equation (3) is obtained by trial and error on data produced by realizations of the model described in Sect. 7.1. Figure E.1 shows realizations for $\sigma_0 = 100$ km s$^{-1}$ and six different number of clouds per spaxel, $N$. Each of these realizations deliver a sample dispersion ($\langle\sigma_\parallel\rangle$) mean and median (blue and orange colors, respectively). The curves show solutions for the mean and the median (same respective colors). The exact solution for the mean is (Kenney & Keeping 1951):

$$ \langle\sigma_\parallel\rangle_{\text{mean}} = \sqrt{\frac{2}{N} \frac{\Gamma(\frac{N}{2})}{\Gamma(\frac{N-1}{2})}} \sigma_0, \qquad (E.1) $$

where $\Gamma(n)$ is the complete gamma function. In this paper we use instead the approximation for the median given by Eq. (3).
<figure><img src="image_13.png" /><figcaption>Fig. E.1. Mean and median sample dispersion versus number of scores per spaxel for σ₀ = 100 km s⁻¹.</figcaption></figure>

## Appendix F: On the velocity auto-correlation function

### F.1. Error budget

To account for the variance of $\xi(\Delta r_{\perp}, \Delta v)$ due to uncertainties in the velocity measurements, $\sigma_{\xi}^2$(measurement), we ran 500 realizations of $\xi(\Delta r_{\perp}, \Delta v)$ by varying velocities according to $N(v, \delta_v^2)$ and computed the variance of the median statistics. On the other hand, the variance of $\xi(\Delta r_{\perp}, \Delta v)$ due to shot noise, $\sigma_{\xi}^2$(Poisson), was computed according to

$$ \sigma_{\xi}^{2}(\text{Poisson}) = (DD/N_{DD}^{2} + RR/N_{RR}^{2})/(RR/N_{RR})^{2}, \quad (\text{F.1}) $$

where DD (RR) is the number of data-data (random-random) pairs in a ($\Delta r_{\perp}$, $\Delta v$) bin and $N_{DD}$ ($N_{RR}$) is the total number of respective pairs.

The variance of $\xi(\Delta r_{\perp}, \Delta v)$ was taken to be the sum in quadrature of the above variances. This is dominated by $\sigma_{\xi}^2$(Poisson) due to both the small number statistics and the small velocity errors compared to the velocity bins (Figs. 5 and 6).

### F.2. Random catalogs

Sampling a uniform distribution $U_{[a,b]}$ (of variance $\sigma_u^2 = (b-a)^2/12$) results in a normal distribution of variance $\sigma_N^2 = \sigma_u^2/N$, where N is the number of averaged velocities. This implies that even a uniform velocity field will show velocity correlation if recorded with a beam that does not resolve the clouds spatially. Taking advantage of this property, we draw random ARCTOMO velocities from $N(0, \sigma_{RR}^2)$, where the random variance $\sigma_{RR}^2 = (8000^2/12)/N$, that is, $\sigma_{RR}^2$ is the binned variance of $U_{[-4000,4000]}$. For some systems we extend this range until convergence is achieved.

### F.3. Mock catalogs

Despite the safeguards to account for each system's complex selection function, some limitations remain: (1) random non-detections change the polygon edges (Hamilton & Tegmark 2004) and thus bias $\xi$ close to the edges in the spatial direction; (2) individual data errors (both in velocity and spatially) can add spurious signal in $\xi(\Delta r_{\perp}, \Delta v)$.

To address these issues we build realistic mock catalogs by creating re-sized replicas of the data samples, exactly like in Sect. 6.3, but having a random uniform $v$ in the range (-4000, +4000) km s⁻¹, random RA-DEC normally distributed around spaxel centers with standard deviation 2.4 kpc (two delensed spaxels), and random $W_0$ drawn from the quasar $W_0$ distribution (Cooksey et al. 2013). As in Sect. 6.3, detections and non-detections are defined by confronting $W_0$ with the detection limit set by the S/N. Each mock entry passes the same S/N criterion outlined above Fig. F.1 shows $\xi(\Delta r_{\perp}, \Delta v)$ measured for a mock data set, based on system 1527 a. As expected, there is no significant clustering power at any scale, which lends support to the technique and also to the fidelity of the random catalogs.

On the other hand, to re-create the spurious signal that a single velocity across spaxels would produce, mock catalogs are created as above but having velocities distributed as $N(0, \delta_v^2)$, where $\delta_v$ are the individual measurement errors. The dashed curves in Fig. 6 show the resulting $\xi(\Delta r_{\perp}, \Delta v)$.

<figure><img src="image_14.png" /><figcaption>Fig. F.2. Per-system 1-$\sigma$ width of $\xi^{\text{arc}}(\Delta v)$ displayed in Fig. 6. This is computed via a single Gaussian fit of $\xi^{\text{arc}}(\Delta v)$ (x-axis) or predicted by the kinematic model described by Eq. (8), using $\sigma_0$ and N in Table 4 (y-axis). Error bars represent fit and propagated errors, respectively.</figcaption></figure>
<figure><img src="image_15.png" /><figcaption>**Fig. F.1.** $\xi(\Delta r_{\perp}, \Delta v)$ measured for a realistic mock data set with no correlation. (a) In the scatter plot the blue (gray) dots correspond to the mock data (random) pairs, and the side panels show the respective projections in number counts. In the 2d-histogram the color scale indicates $\xi(\Delta r_{\perp}, \Delta v)$ between -1 and 1. (b) Velocity projection of $\xi(\Delta r_{\perp}, \Delta v)$. The shaded pink regions indicate $1\sigma$ errors.</figcaption></figure>