An Ensemble Study of Turbulence in Extended QSO Nebulae at z ≈ 0.5–1

Turbulent motions in the circumgalactic medium play a critical role in regulating the evolution of galaxies, yet their detailed characterization remains elusive. Using two-dimensional velocity maps constructed from spatially extended [O ii] and [O iii] emission, Chen et al. measured the velocity structure functions (VSFs) of four quasar nebulae at z ≈ 0.5–1.1. One of these exhibits a spectacular Kolmogorov relation. Here, we carry out an ensemble study using an expanded sample incorporating four new nebulae from three additional quasi-stellar object (QSO) fields. The VSFs measured for all eight nebulae are best explained by subsonic turbulence revealed by the line-emitting gas, which in turn strongly suggests that the cool gas (T ∼ 104 K) is dynamically coupled to the hot ambient medium. Previous work demonstrates that the largest nebulae in our sample reside in group environments with clear signs of tidal interactions, suggesting that environmental effects are vital in seeding and enhancing the turbulence within the gaseous halos, ultimately promoting the formation of the extended nebulae. No discernible differences are observed in the VSF properties between radio-loud and radio-quiet QSO fields. We estimate the turbulent heating rate per unit volume, Q turb, in the QSO nebulae to be ∼10−26–10−22 erg cm−3 s−1 for the cool phase and ∼10−28–10−25 erg cm−3 s−1 for the hot phase. This range aligns with measurements in the intracluster medium and star-forming molecular clouds but is ∼103 times higher than the Q turb observed inside cool gas clumps on scales ≲1 kpc using absorption-line techniques. We discuss the prospect of bridging the gap between emission and absorption studies by pushing the emission-based VSF measurements to below ≈10 kpc.


INTRODUCTION
The circumgalactic medium (CGM) is the outermost, gaseous envelope of a galaxy, extending beyond the visible stellar disk and containing the majority of the baryons in the galaxy.This main gas reservoir records critical information about a galaxy's past and ongoing interactions with the surrounding environment.Due to the tenuous nature of the CGM, absorption spectroscopy using bright background sources -predominantly quasi-stellar objects (QSOs) -has been the main probe of gaseous halos, yielding sensitive constraints on gas density, temperature, metallicity, and ionization state (e.g., Chen 2017;Tumlinson et al. 2017;Rudie et al. 2019;Péroux & Howk 2020;Donahue & Voit 2022;Faucher-Giguère & Oh 2023).
Over the past decade, the advent of wide-field, highthroughput integral field spectrographs (IFSs) has provided a spatial resolving power that complements the pencil-beam probe from QSO absorption spectroscopy, greatly aiding in the investigation of the CGM.Various dynamical processes in the CGM, such as infalls, outflows, and tidal interactions, can now be spatially and spectrally mapped by IFSs via strong nebular emission lines (e.g., Epinat et al. 2018;Johnson et al. 2018;Rupke et al. 2019;Chen et al. 2019).One particularly exciting prospect with these resolved kinematic measurements is the robust constraint of turbulent motions in low-density gas.
With a high Reynolds number, ionized, diffuse plasma such as the CGM is expected to be turbulent (see e.g., Burkhart 2021, for a recent review), which can manifest as large density Chen et al. fluctuations commonly observed in extended emission at tens of kpc scales in gaseous halos (e.g., Travascio et al, in prep).Turbulence plays a critical role in several key processes in the CGM, such as mixing/transporting metals (e.g, Pan & Scannapieco 2010), facilitating multiphase structure formation (e.g., Gaspari et al. 2018;Fielding et al. 2020), and offsetting radiative cooling (e.g., Zhuravleva et al. 2014).Until recently, observing turbulence in circumgalactic/intergalactic gas has had to rely on two approaches employing high-resolution absorption line spectra of background QSOs.One approach is to observe line widths of ions with different masses and isolate the turbulent contribution to the velocity profile along the line of sight (e.g., Rauch et al. 1996;Rudie et al. 2019;Qu et al. 2022;Chen et al. 2023a).Alternatively, if multiple lines of sight (e.g., to gravitationally lensed QSO images) are available, turbulence can be measured as a function of transverse separation between the lines of sight to form the structure functions for the line-of-sight velocities (Rauch et al. 2001).With the advent of IFS, spatially-resolved velocity maps of entire gaseous galactic halos can now be obtained in one shot, enabling the simultaneous measurement of the turbulent power spectrum over a wide range of scales, thus providing multiple independent constraints on the nature of turbulence and the turbulent energy transfer in the gas.
Recently, Chen et al. (2023b) (hereafter Paper I) obtained two-dimensional line-of-sight velocity maps of line-emitting gas around four QSOs up to scales of ∼ 100 kpc using IFS data.Taking advantage of the spatially-resolved velocity maps from IFS observations, these authors constructed velocity structure functions (VSFs),   , defined as where x and r represent, respectively, a position in the velocity map and the distance vector between two positions separated by r.The exponent  is generally referred to as the order of the VSFs, and ⟨⟩ denotes the mean value averaged over all available velocity pairs separated by .As can be seen from the definition,   quantifies the scale-dependent variance of a velocity field (e.g., Frisch 1995), and has been commonly used to probe the dynamical state of the interstellar medium (ISM) in local H II regions (e.g., Wen & O'dell 1993;Ossenkopf & Mac Low 2002;Federrath 2013;Arthur et al. 2016;Padoan et al. 2016;Melnick et al. 2021;García-Vázquez et al. 2023) as well as in the intracluster medium (ICM) in nearby cool-core clusters (Li et al. 2020;Ganguly et al. 2023).
While the uncertainties remained large for three QSO nebulae, Paper I found that in one particular nebula, the gas dynamics can be unambiguously characterized by the Kolmogorov relation, expected for subsonic, isotropic and homogeneous turbulent flows.Building upon the sample studied in Paper I, in this follow-up paper, we include results from four nebulae discovered in three new QSO fields.Combining this new sample with the previous one establishes a sample of eight QSO nebulae that allows us to carry out an ensemble study of the empirical properties of CGM turbulence in distant QSO host halos.The QSOs are all luminous, with a bolometric luminosity of ∼ 10 47 erg s −1 , and span a range in redshift from  ≈ 0.5 to  ≈ 1.1.The nebulae are revealed in [O II]  3727, 3729 and/or [O III]  5008 line emission (see Figure 1) and are selected to have an extended, contiguous emission area ≳ 1500 kpc 2 .Table 1 summarises the properties of the QSOs in the sample.Out of the seven QSOs, four are radio-loud, and three are radio-quiet.
This paper is organized as follows.In Section 2, we describe the observations of the ensemble sample and the subsequent velocity measurements using the emission line features.Based on the spatially-resolved velocity maps, we present the VSFs for all eight nebulae in Section 3. We discuss the implications of the results in Section 4 and conclude in Section 5. Throughout this paper, we adopt a flat ΛCDM cosmology with  0 = 70 km s −1 Mpc −1 , Ω M = 0.3 and Ω Λ = 0.7 when deriving distances, masses and luminosities.All distances quoted are in physical/proper units.

OBSERVATIONS AND DATA ANALYSIS
To constrain the turbulent energy spectrum, we follow the approach described in Paper I to construct the VSFs of four nebulae found in three new QSO fields, PKS 0405−123, HE 0238−1904, and PKS 0552−640.In this section, we briefly summarize the IFS observations and the steps we took to construct a spatially-resolved velocity map based on a line profile analysis of [O II]3727, 3729 and [O III]5008 emission lines in these QSO fields.

IFS observations
To measure the spatially-resolved kinematics in the plane of the sky for the QSO nebulae in our sample, we use the IFS observations obtained using the Multi-Unit Spectroscopic Explorer (MUSE; Bacon et al. 2010) on the VLT UT4.The Wide-Field-Mode (WFM) was used to observe all seven fields, offering a field-of-view (FOV) of 1 ′ × 1 ′ for a single pointing and a spatial sampling of 0. ′′ 2 per pixel.MUSE covers a wavelength range of 4750-9350 Å and has a spectral resolving power of  ≈ 2000-4000, with a higher resolution at longer wavelengths.Note- Atmospheric seeing FWHM measured using the QSO at 7000 Å.
To improve the quality of line fitting, each combined data cube was convolved with a Gaussian kernel of FWHM= 0. ′′ 7.This yielded a total PSF FWHM of ≈ 0. ′′ 9-1.′′ 0, corresponding to a projected separation of 6-9 kpc at the redshifts of these QSOs.

Construction of velocity maps
As described in Paper I, the main steps to construct a two-dimensional velocity map include removing the contamination from the QSO point spread function (PSF), subtracting continuum flux across the whole MUSE FOV, constructing optimally extracted narrow-band images for [O II]  3727, 3729 and [O III]  5008 lines using threedimensional masks, and finally fitting Gaussian components to the emission signals and optimizing the parameters via an MCMC analysis.Readers can find the detailed descriptions and associated technical considerations of each step in Paper I. Note that to increase the signal-to-noise ratio for faint spaxels in the outskirts of a nebula, we smooth the data cubes in the spatial dimension with a two-dimensional Gaussian kernel of full-width-at-half-maximum of FWHM = 0. ′′ 7, leading to a total PSF FWHM of ≈ 0. ′′ 9-1.′′ 0 (see Table 2), corresponding to ≈ 6-9 kpc at the QSO redshifts.
A subset (≈10-20%) of spaxels in the nebulae (mostly towards the inner region in the vicinity of the QSOs) exhibit multiple velocity components, which can be identified clearly with the [O III]  5008 line.With MUSE spectral resolution and due to the doublet nature of the [O II]  3727, 3729 line, multiple velocity components are only obvious for narrow features with a velocity dispersion ≲ 50 km/s.In Paper I, we demonstrated that different ways of handling the multicomponent spaxels (e.g., adopting the flux-weighted mean velocity versus using the velocity of the strongest component) do not lead to significant differences in the VSF measurements.The insensitivity of the VSFs to the treatment of multi-component spaxels can be attributed to the relatively small proportion of spaxels requiring a multi-component fit, and that the majority of such spaxels exhibit a single prominent component that dominates the kinematics.Therefore, we opt to take the simple approach of using a single Gaussian function when fitting the lines.
We also treat [O II]  3727, 3729 and [O III]  5008 from the same spaxels separately when conducting the line fitting, allowing the two lines to have different velocities and line widths.This decision is motivated by the observation that for spaxels requiring multiple velocity components, there exists spatial variation in the [O III]/[O II] ratio across different components, resulting in a different flux-weighted mean velocity for the two lines.In addition, the two lines have different footprints within the same nebula due to different signal-to-noise ratios and emission strengths.Therefore, to keep the analyses simple without sacrificing the accuracy of the velocity measurements, we opt to measure [O II] and [O III] separately.

VSF measurements
For the three new QSO fields presented in this paper, we show the continuum-and QSO-subtracted narrow-band images in Figure 1.The narrow-band images for PKS0454−22, J0454−6116, J2135−5316 and TXS0206−048 have already been presented in Figure 1 of Paper I.
As described in Section 3.5 of Paper I, to ensure the robustness of the VSF measurements, we exclude spaxels with a velocity uncertainty larger than 45 km/s.We also examine the velocity map for each nebula in tandem with the broadband images from either MUSE or HST to identify spaxels that are likely to originate from continuum sources.If such spaxels exhibit distinctly different velocities and line widths from the rest of the nebula, we exclude them because such continuum sources are likely to be separate from the rest of the nebula, and are simply projected to be within the nebula footprint.Finally, we exclude spaxels that are outliers (≈ 2 per cent tail on both the blue and red ends) in the probability density distribution of the velocities in each field.After the above-mentioned steps, all spaxels left in the velocity maps are included in the subsequent VSF calculation, as shown in the top left panels of Figures 9-16.Summing over all spaxels included in the VSF analyses, the total luminosity and area for each nebula are listed in Table 3 1 and 9-12; also see Figure 2 of Johnson et al. 2018).For the purpose of this paper, we analyze the southern and eastern nebulae of PKS0405−123 separately and refer to them as PKS0405−123 S and PKS0405−123 E, and we do not include the nebula immediately surrounding the QSO in this field due to its relatively small size.
We measure the VSFs up to order  = 6 for all eight nebulae following the definition of Equation 1. VSFs with  > 6 become too noisy to provide meaningful constraints.Due to the spatial correlation between spaxels that are separated by distances less than the size of the total PSF, not all velocity pairs in each distance bin are independent.Therefore, to obtain a more robust estimate of the uncertainty in the VSF measurements, we adopt the modified bootstrap method described in Paper I. In addition, as shown in Paper I, the spatial correlation due to atmospheric seeing and the additional Gaussian smoothing applied to the datacubes preferentially removes power from small scales and steepens the VSFs.This smoothing effect can be explicitly accounted for by employing a Gaussian-convolved 2nd-order VSF,  ′ 2 , where Γ ′ is a Gaussian-convolved velocity autocorrelation function, Here Γ() and Γ  () are the autocorrelation function of the velocity field and the smoothing kernel, respectively.A more detailed derivation for Equation 3 can be found in Equations 2-7 of Paper I.
To quantify the slopes of the 2nd-order VSFs, we adopt a single power-law model: When fitting the observed  ′ 2 with a power-law model, we conduct the convolution in Equation 3 numerically, and find the best-fitting  2 with the Scipy curve fit routine for each of the 1000 modified bootstrap samples described above to obtain the mean and dispersion of  2 .Note that we only Neb.S.
Neb.E. consider non-negative slopes of  2 , which is motivated by data and avoids divergent values at  = 0.With the IFS data, observations are confined to projected quantities both in velocity and spatial separations.Therefore, we report the VSF measurements using the line-of-sight velocities and the projected spatial separation  proj in the plane of the sky.The potential limitations due to the projection effect will be discussed in further detail in Section 4.5.

RESULTS
Using the velocity maps constructed for individual nebulae, we proceed with the VSF analysis using the full sample of eight extended nebulae.Recall that while it is relatively straightforward to measure the VSFs using spatially-resolved velocity maps, a primary systematic uncertainty is possible contributions to the observed signal from coherent bulk motions projected in the plane of the sky.To account for this uncertainty, we follow the approach adopted by Paper I and consider a simple, unidirectional velocity gradient model parameterized as (, ) =  +  + , where  and  are coordinates of individual spaxels, and , , and  are the free parameters.For each nebula, we measure the VSFs with and without the best-fitting two-dimensional bulk-flow model subtracted.The amplitudes of the best-fitting gradient for the [O II]  3727, 3729 and [O III]  5008 emission lines in each field are listed in Table 4.We estimate the uncertainty of this velocity gradient by repeating the fitting with 1000 randomly perturbed velocity maps based on the MCMC chains for each spaxel and find that the uncertainties are small (≲ 0.1 km s −1 ) for all nebulae.Therefore, we do not list the uncertainties in Table 4.To identify possible coherent motions dominant along the radial or tangential directions (for example in the case of strong outflows or inflows), Paper I also calculated the VSFs using radial and tangential velocity pairs separately and found the VSF measurements to be comparable along these two directions.For the newly analyzed nebulae in this pa- Separation r proj (kpc) [OII] grad.removed 10 1 10 2 Separation r proj (kpc) [OIII] 10 1 10 2 Separation r proj (kpc) [OIII] grad.removed  The observed 2nd-order VSF  ′ 2 () for the nebulae of radio loud QSOs.Vertical dashed lines mark the size of the total PSF FWHM for each field.The data points and the error bars show the mean and the standard deviation for the 1000 measurements obtained through the modified bootstrap method (see Section 2).The dashed and solid gray lines show, respectively, the expected  2 for Kolmogorov turbulence before and after convolving with an appropriate PSF.As different fields have slightly different PSF sizes, we use the mean value of the PSF FWHM for all radio-loud fields included in the panels when constructing the expected Kolmogorov VSF.The dotted gray lines show a power-law with a slope of 1 (i.e., Burger's turbulence).The four panels from left to right show the results using the [O II] and [O III] lines, both from direct measurements and after removing a uni-directional gradient (see text).Bottom row: Same as the top row but for radio-quiet fields.All nebulae exhibit VSFs with an increasing amplitude in velocity variance as a function of separation distance.PSF smoothing significantly steepens the apparent slopes of the VSFs and we will explicitly take into account this smoothing effect when fitting the VSFs with a power-law (see Section 3.2).per, we find a similar trend where radial and tangential VSFs show no clear differences and are therefore not included in the presentation here.
In this section, we first examine the general trend displayed in the second-order VSF across all eight QSO nebulae.Then we quantify and compare the best-fitting VSF slope over a finite range of spatial scale where the measurements can be characterized by a power-law function.Finally, we explore the presence or absence of extended self-similarity (ESS; see, e.g., Benzi et al. 1993) in turbulent flows in QSO host halos by measuring the higher-order VSFs.

The overall shape of VSFs
Figure 2 shows the observed 2nd-order VSFs,  ′ 2 , for the eight nebulae in our sample.Radio-loud and radio-quiet fields are shown in the top and bottom rows, respectively.The vertical dashed lines mark the FWHM of the total PSF for each field (see Table 2).To guide the visual comparison, we overplot the expected  2 for Kolmogorov turbulence, with the dashed gray line showing the intrinsic 2/3 slope and the solid gray line showing the observed shape of  ′ 2 after convolving with an appropriate PSF.Because different fields have slightly different PSF sizes, we use the mean value of the PSF FWHM for radio-loud (-quiet) fields when constructing the expected Kolmogorov  ′ 2 for the top (bottom) row.We also show the power-law with a slope of 1 (e.g., Burger's turbulence), without convolving with a PSF, in dotted gray lines.The comparison between the data and the model  2 with slopes 2/3 and 1 underlines the importance of including the PSF effect when quantifying the observed VSF slopes.In particular, if the probed distance separation,  proj , is ≲ 10-20 times the PSF FWHM, the PSF smoothing effect can significantly steepen the apparent slope of the VSFs and a naive visual inspection will lead to the wrong conclusion that the VSF slopes are steeper than their intrinsic values.The VSFs obtained using the gradient-removed velocity maps are also included in Figure 2 for comparison.As shown in Figure 2, all nebulae in our sample exhibit an overall increasing trend of velocity fluctuations with increasing spatial scale.The values of ⟨Δ 2 ⟩ range from ≈ 5000-10,000 km 2 /s 2 at  proj ≈ 10 kpc to ≈ 10, 000-80,000 km 2 /s 2 at  proj ≈ 50 kpc.The results based on the [O II]  3727, 3729 and [O III]  5008 lines are consistent within the uncertainty for fields with both lines.In general, we do not expect the VSFs constructed from [O II]  3727, 3729 and [O III]  5008 lines to be identical, because the footprints of the two emission lines in the nebulae do not overlap completely due to the different signal-to-noise ratios of the two lines at different locations.For regions with overlapping footprints from both [O II]  3727, 3729 and [O III]  5008 emission, the line-of-sight velocities can also differ for spaxels with multiple velocity components and varying [O III]/[O II] line ratios between components, as discussed in Section 2.2.We will show below that the VSFs from [O II] and [O III] lines lead to consistent constraints on the dynamical state of the gas.In addition, the removal of a large-scale, unidirectional velocity gradient generally flattens the VSFs via preferentially reducing the power at larger distance separations.Nonetheless, the constrained slopes for a single power-law fit are consistent before and after the removal of the gradient, as we will discuss in the following section.

2nd-order VSF slopes
As shown in Figure 2, all VSFs exhibit structures that deviate from a single power-law.In particular, at larger separations, the VSFs can show an overall decreased amplitude (e.g., TXS0206−048 at  proj ≳ 60 kpc), an overall enhanced power (e.g., J0454−6116 at  proj ≳ 30 kpc), or an oscillatory behavior (e.g., HE0238−1904 at  proj ≳ 30 kpc).Such deviations may reflect different levels of velocity fluctuations in the central regions of the nebulae versus the outskirts, as velocity pairs at larger separations are predominantly constructed from spaxels in the outskirts.In addition, large-scale periodic oscillations in the velocity fields can manifest as oscillations in the VSFs at large separations (e.g., García-Vázquez et al. 2023).The VSF measurements at larger separations are also more uncertain due to a combined effect of fewer pair counts and uncertain velocity centroids as a result of fainter signals in the outskirts of a nebula.
Taking into account the above-mentioned factors, we restrict the fitting to be within a finite range of spatial scales, [ 1 ,  2 ], when employing a single power-law model to quantify the slopes of the VSFs.The lower limit  1 is chosen to be the FWHM of the total PSF for each field (see Table 2), while the upper limit  2 is chosen through a series of trial and error such that we obtain the lowest reduced  2 for the best-fitting model within this range.We refer to  2 as the VSF turnover scale and will discuss its correlation with the nebula size later in Section 4.4.When constraining the VSF slopes, we explicitly incorporate the smoothing effect in the 2nd-order VSF models before comparing them with the data, as described in Section 2.3.The [ 1 ,  2 ] values as well as the best-fitting slopes for the 2nd-order VSF,  2 , are listed in Table 4 using both the directly measured line-of-sight velocity maps and the gradient-removed velocity maps.As mentioned above, removing a large-scale, unidirectional gradient tends to flatten the VSF, leading to a smaller  2 and weaker constraints on the VSF slopes.The comparisons between best-fitting powerlaw models and the data for PKS0454−22, J0454−6116, .PKS 0405−123 S exhibits a steeper slope but this is also a system that shows a large-scale velocity gradient across the nebula.After removing a unidirectional velocity gradient, the VSF is consistent with the Kolmogorov expectation.Below we discuss these three categories individually. Nebulae Note- The best-fitting slopes are derived from 1000 modified bootstrap samples, as discussed in Section 2. These slopes correspond to the intrinsic power-law slopes for  2 , with our fitting process explicitly addressing the PSF smoothing effect in the measured  ′ 2 .The reported values are medians along with the 16 th and 84 th percentiles.The 3 th and 97 th percentiles are approximately double the uncertainty estimates listed here for all fields.For the unconstrained results, we present 95% upper limits for the slope, assuming the observed pair separations fall within the inertial range.If the available pair separations are close to injection scales, then no robust constraints can be obtained.We exclusively consider non-negative power-law slopes, in line with the discussion in Section 2.  Measurements obtained after removing a 2D velocity gradient (see § 3). Lower and upper bounds in the projected distance separation,  proj , in the unit of kpc, within which the power-law slopes of the VSFs are constrained (see § 3.2). Best-fitting 2D velocity gradient, in the unit of km/s/kpc.the steepening of the VSF in this nebula is a strong effect of projection smoothing if the depth of the nebula is larger than the projected distance scales in the plane of the sky (see discussions in Section 4.5).Moreover, the line-of-sight velocity maps for both [O II] and [O III] show a possible velocity shear along the NW-SE direction (see Figures 9 and 10).The bestfitting direction and amplitude for the velocity gradient are consistent between both emission lines, suggesting that the bulk flow can plausibly contribute to the VSF measurements, leading to steeper VSF slopes.Indeed, the VSFs become flatter after we remove a unidirectional velocity gradient, resulting in slope upper limits consistent with the Kolmogorov expectation albeit with larger uncertainties.
In summary, using the direct measurements of the line-ofsight velocity fields based on the [O II]  3727, 3729 and/or [O III]  5008 emission lines, five out of eight nebulae exhibit a 2nd-order VSF slope that is consistent with the expected value of 2/3 for Kolmogorov turbulence while three exhibit a flatter VSF.Incidentally, the three nebulae with a flatter VSF are also the smallest in the sample (see Table 3).It is possible that the observations do not have a sufficiently large dynamic range for securing a robust constraint on the shape of the VSF (see, e.g., Federrath et al. 2021).

Extended self-similarity (ESS) in turbulent flows
In addition to measuring the 2nd-order VSF slope  2 , Paper I also explored the presence of ESS, in which a simple power-law function holds between VSFs of different orders on spatial scales that are outside of the inertial range where the Kolmogorov relation applies (see, e.g., Benzi et al. 1993).This ESS is particularly useful for inferring the energy cascade rate when the inertial range is not well established.Compared with the slopes of VSFs of individual orders, the ESS slope ratios are often better constrained with a higher statistical significance thanks to the tight correlation between different orders.In addition, an enhanced level of intermittency in a velocity field will suppress the VSF slopes at higher orders compared with the slopes of lower orders (e.g., Frisch 1995), making the ESS slope ratios a valuable diagnostic for the underlying gas dynamics.Here we explore the presence or absence of ESS in the QSO nebulae by measuring the VSFs up to order  = 6.We obtain the slope ratios   / 3 for  = 1-6 by fitting a single power-law model to the  ′  vs.  ′ 3 measurements.As discussed in Paper I, the smoothing effect due to the data PSF does not change the ESS slope ratios.The results are displayed in Figure 3, where the data points represent the median values obtained from fitting the   5008 measurements, as well as their corresponding velocity residual maps after removing a coherent, unidirectional gradient.  / 3 is the best-fitting power-law slope for the relation between the observed th-order VSF,  ′  , and  ′ 3 .The data points represent the median values obtained from fitting the 1000 modified bootstrap samples (see Section 2.3), and the error bars indicate the 16 th and 84 th percentiles.The solid curves represent the expected ratio of   / 3 for subsonic Kolmogorov turbulence, taking into account the intermittency correction presented in She & Leveque (1994).The dashed curves represent the expected ratio for supersonic magnetohydrodynamic turbulence, as presented in Boldyrev (2002).The dash-dotted (dotted) curves indicate the   / 3 ratio derived from numerical hydrodynamic turbulent simulations for a Mach number of 0.9 (6.1), as presented in Pan & Scannapieco (2011).Finally, the blue dash-dotted curves represent the expected   / 3 ratio for Kolmogorov turbulence without the intermittency correction, scaling simply as /3.The top row shows the results for radio-loud fields, while the measurements for radio-quiet fields are shown at the bottom.Except for the field of HE0238−1904, all nebulae exhibit ESS slope ratios in agreement with expectations for subsonic motions using directly measured velocity maps and/or velocity residual maps after removing a coherent unidirectional gradient.None of the nebulae show signatures of supersonic motions with a Mach number ≳ 6. 1000 modified bootstrap samples (see Section 2.3), and the error bars indicate the 16 th and 84 th percentiles.The correlation between 2nd-and 3rd-order VSFs for each nebula are displayed in the right-most panels of Figures 9-16.
Specifically, we measure   / 3 using the [O II]  3727, 3729 and [O III]  5008 velocity maps as well as their corresponding residual maps after removing a unidirectional velocity gradient.Figure 3 shows the ESS slope ratios, with radio-loud fields in the top row and radioquiet fields at the bottom.We also overplot the expected   / 3 ratios from different theoretical considerations and numerical simulations, including the Kolmogorov expectation of   / 3 = /3 (blue dashed curve), the Kolmogorov turbulence with intermittency correction (solid curve; She & Leveque 1994), the expectation for supersonic magnetohydrodynamic turbulence (dashed curve; Boldyrev 2002), and numerical predictions for hydrodynamic turbulence with Mach numbers of M = 0.9 and 6.1 (dash-dotted and dotted curves; Pan & Scannapieco 2011).In general, the ratio   / 3 is expected to be suppressed significantly at larger 's in supersonic flows with a high Mach number.This can be seen in Figure 3 where the numerical simulations predict that for gas motions with M = 6.1,   / 3 does not increase significantly for  > 3, showing a plateau in the   / 3 curve (dotted lines).
While the strongest distinguishing power for different scenarios comes in at higher orders, the measurements are also more uncertain.In addition, removing a large-scale gradient from the velocity field can change the   / 3 ratios to be more consistent with predictions for lower Mach numbers (e.g., see the trend for PKS0405−123 S).Within the 16 th and 84 th measurement percentile range and considering the results both before and after removing the large-scale velocity gradient, seven out of eight nebulae in our sample show ESS slope ratios consistent with expectations from subsonic turbulence (black solid curve, blue dash-dotted curve, and dash-dotted curve in Figure 3).For the nebula surrounding HE0238−1904, the   / 3 ratios are consistent with the predictions for supersonic magnetohydrodynamic turbulence as presented in Boldyrev (2002), suggesting that the Mach number of gas motions in this field may be higher than that in other nebulae.Given that this field has a constrained  2 value that is consistent with the Kolmogorov expectation as discussed above, additional effects (e.g., the presence of a dynamically important magnetic field) might contribute to a relatively small  2 in tandem with suppressed   / 3 ratios.A more detailed investigation into the properties of this nebula (e.g., ionization state, interactions with group member galaxies) is needed to further shed light on the possible physical causes for this difference in ESS slope ratios, and a larger sample is required to examine whether the HE0238−1904 nebula is a special case.Overall, no system in our sample exhibits ESS slope ratios that indicate gas motions with M ≳ 6.

DISCUSSION
We have shown that the 2nd-order VSF measured for eight QSO nebulae in our sample exhibits a range of slopes.While five of the nebulae in our sample are consistent with the expected slope of 2/3 for Kolmogorov turbulence, the remaining three exhibit a shallower slope.Despite a range of 2nd-order VSF observed in these QSO nebulae, the measurements suggest that turbulent flows in the [O II]  3727, 3729 and [O III]  5008 line-emitting clouds are subsonic.The subsonic dynamical state of the gas is further corroborated by the ESS slope ratios   / 3 , which are consistent with theoretical or numerical expectations for subsonic systems with M ≲ 1 in seven out of eight nebulae.None of the systems shows   / 3 measurements that are indicative of highly supersonic flows with M ≳ 6.
In addition, we do not observe significant differences between radio-loud and radio-quiet QSO fields in terms of nebula size, line emission luminosity, VSF slopes, VSF amplitude, and turbulent energy heating rate.Recall that five of the nebulae in our sample occur near radio-loud QSOs, while the remaining three reside in radio-quiet halos.The main distinguishing characteristic between radio-loud and radio-quiet QSOs is the presence of powerful jets in radio-loud sources that can result in large-scale structures like radio lobes spanning from tens to thousands of kpc in size (e.g., Mullin et al. 2008).The mechanical energy contained in the collimated jets and the associated inflated bubbles is estimated to be ∼ 10 41 -10 46 erg/s (e.g., Heckman & Best 2014).If a significant portion of this energy can be deposited into the CGM as kinematic energy, we may expect the VSFs from radio-loud and radio-quiet fields to exhibit different properties.While previous studies have found that radio jets are the dominant mechanism for driving fast outflows in the inner ≲ 10 kpc regions in radio galaxies (e.g., Nesvadba et al. 2017), a lack of correlation between the observed VSFs and the radio power suggests that the effect of radio jets may be limited to the inner regions and have little influence on the gas kinematics on scales ≳tens of kpc.This is in agreement with simulation predictions for the ICM in cool-core clusters (e.g., Yang & Reynolds 2016).A larger sample with both radio-loud and radio-quiet sources will be helpful to draw robust conclu-sions regarding the difference (or lack thereof) in the CGM dynamics between these two populations.
In this section, we first discuss the implications for the dynamical state of the gas in the multiphase CGM and infer the energy transfer rate in these QSO host nebulae.We then discuss potential caveats associated with observational limitations, including projection effects, finite nebula sizes, and the small number of systems in the current sample.

Implications for the multiphase CGM dynamics
Based on the velocity dispersion of member galaxies in the QSO host group environment, the halo mass of the QSO hosts in our sample is estimated to be ≈ 10 13 -10 14 M ⊙ (see e.g., Johnson et al. 2018;Helton et al. 2021;Johnson et al. 2022;Liu et al. 2023).This mass range suggests a viral temperature of  ≈ 10 6 -10 7 K for the underlying hot halo (e.g., Mo et al. 2010).Meanwhile, the sound speed of the gas can be calculated by  s = √︁  B / p , where  = 5/3 is the adiabatic index for an ideal monatomic gas,  B is the Boltzmann constant,  is the mean atomic weight (which is 0.588 for fully ionized gas), and  p is the proton mass.For the cool gas of  ≈ 10 4 K,  s,cool ≈ 15 km/s, while for the hot medium of  ≈ 10 6 -10 7 K,  s,hot ≈ 150-500 km/s.Therefore, for the nebulae in our sample, the Mach number calculated using the sound speed of the cool gas is M cool = √ 3 pos / s,cool ≈ 7-18, and M hot = √ 3 pos / s,hot ≈ 0.2-1.8 using  s,hot for the hot gas.Here  pos is the velocity dispersion in the plane of the sky.As we will discuss below in Section 4.3,  pos is typically smaller than the velocity dispersion along the line of sight, and the Mach numbers will be larger (M cool ≈ 9-20 and M hot ≈ 0.3-2.0)when estimated using the line-of-sight velocity dispersion.
Given the contrast between the two Mach numbers, M cool and M hot , the subsonic motions revealed by the VSFs of the nebulae suggest that the [O II] and [O III] emission originates from cool gas clumps embedded in the ambient hot medium.If these cool clumps are in pressure equilibrium with the hot halo, then they can serve as tracers for the kinematics of the volume-filling plasma.The scenario of a dynamicallycoupled multiphase gaseous system is supported by absorption line studies on CGM kinematics of  ∼ 2 star-forming galaxies (e.g., Rudie et al. 2019) as well as by recent measurements in the core regions of nearby galaxy groups and clusters (e.g., Li et al. 2020;Olivares et al. 2022).There has also been an increasing number of theoretical and numerical predictions arguing for a shared dynamical state across different gas phases (e.g., Gaspari et al. 2018;Gronke & Oh 2018;Schneider et al. 2020;Mohapatra et al. 2022) The dynamical coupling likely happens due to a combination of physical processes involving cooling, the exchange of mass and momentum between cool and hot phases, and the competition between cool clump formation and cloud crush-ing at different mass/length scales.Turbulence is expected to facilitate these processes, which in turn further feed into the development of turbulence in the gaseous halo.In the absence of turbulence, the condensed cool clumps tend to settle in more organized structures such as a disk.The extended morphological features of the nebulae in our sample suggest that turbulence is significant in these gaseous halos.Phenomenologically, Gaspari et al. (2018) proposed an empirical criterion of  cool / eddy ≲ 1 for the condensation and survival of cool gas in clusters and groups, where  cool is the gas cooling time and  eddy is the eddy turnover time.Based on the VSF measurements, we can calculate the eddy turnover time via  eddy ≈  −1/3  2/3 , where  is the energy transfer rate per unit mass at the spatial scale  (for more discussion on  see Section 4.2 below).For the nebulae in our sample, we estimate  eddy ≈ 60-150 Myr at  ≈ 10 kpc and  eddy ≈ 150-300 Myr at  ≈ 50 kpc.While we cannot obtain an estimation for  cool due to the absence of temperature and metallicity measurements of the hot phase, our measured  eddy is in agreement with the estimated values ( eddy ≈ 100-200 Myr for galaxy groups) that fulfill the gas condensation criterion in Gaspari et al. (2018) (see their Figure 5).
In addition, turbulence in the CGM can also be produced by Kelvin-Helmholtz instability during the accretion of cool gas streams (e.g., Vossberg et al. 2019;Mandelker et al. 2019;Li et al. 2023), and motions of fragmented cool gas clumps in disrupted, turbulent mixing zones near the accreting streams are predicted to be subsonic in numerical simulations (e.g., Aung et al. 2019).Among our sample, the nebula in the field of TXS0206−048 exhibits compelling signs of cool, filamentary gas accretion from large scales (Johnson et al. 2022), suggesting that the observed subsonic turbulence may be in part produced through the accreting streams.
Finally, previous studies have identified a correlation between the presence of close companions around the QSOs and the presence of strong, extended nebular line emission (see e.g., a narrow-band imaging survey by Stockton & MacKenty 1987).In our sample, the morphokinematics of some nebulae (e.g., PKS0405−123, HE0238−1904, TXS0206−048) reveal that part of the line-emitting gas originates from stripped ISM of group member galaxies as indicated by consistent line-ofsight velocities between the galaxies and extended nebulae (see e.g., Johnson et al. 2018;Helton et al. 2021;Liu et al. 2023).It is natural to assume in these cases that the tidal interactions between group member galaxies disturb the gas and enhance turbulence and thermal instabilities in the hot halo, leading to more efficient cooling and cool clump condensation.The stripped ISM can also serve as massive cool gas seeds that facilitate the coagulation of smaller clumps, aiding in subsequent stochastic mass growth in the cool phase (e.g., Gronke et al. 2022).The significance of this environmental effect on the formation of extended nebulae is supported by  et al. (2023a).The turbulent heating rates in the QSO nebulae, the cool-core cluster ICM, and the star-forming molecular clouds are on average ∼ 1000 times higher than that in Complex C and cool gas clumps probed in absorption.
the fact that the nebulae in PKS0405−123, HE0238−1904, and TXS0206−048 are much larger in area than the nebulae in fields such as J0454−6116 and J2135−5316 where no massive close companions with consistent line-of-sight velocities were found in the nebulae footprint.

Energy transfer rate over seven decades in spatial scale
As described in Paper I, the energy transfer rate per unit mass  can be calculated via the "four-fifths law" (Kolmogorov 1941;Frisch 1995): For Kolmogorov turbulence,  is a constant at all scales within the inertial range.For VSFs flatter (steeper) than the Kolmogorov expectation, the energy transfer rate would be higher (lower) on smaller spatial scales.Across different nebulae in our sample and on different scales between 10-60 kpc, the estimated  shows a range of values between ≈ 0.02 cm 2 s −3 and ≈ 0.2 cm 2 s −3 .For nebulae with both [O II]  3727, 3729 and [O III]  5008 measurements, the values obtained using these two lines are consistent within uncertainty.This estimated range for  with our sample is comparable to the measurements for H filaments in core regions of nearby cool-core clusters (Li et al. 2020;Ganguly et al. 2023) and molecular clouds in nearby H II regions (e.g., Hennebelle & Falgarone 2012).Much lower estimates of  ≈ 10 −7 -10 −3 cm 2 s −3 were obtained for CGM cool clumps probed through absorption line spectroscopy (Rauch et al. 2001;Chen et al. 2023a), and a Milky Way high-velocity cloud (HVC; Marchal et al. 2021).
To gain further insights into the differences between these dynamical systems, we convert the estimated  to a turbulent heating rate per unit volume via  turb = , where  is the gas density and can span a wide range for gas in different phases.For the QSO nebulae in our sample, the [O II]  3727, 3729 doublet line ratios suggest a median upper limit of gas density for the  ∼ 10 4 K cool phase of ≲ 40 cm −3 (Liu et al. 2023), while an estimate of ≈ 1-5 cm −3 is obtained assuming pressure equilibrium between typical AGN-illuminated [O II]-emitting gas and the hot halo (Johnson et al. 2022).Based on the [S II]6716, 6731 doublet ratio, observations of spatially extended nebula illuminated by the active galactic nucleus (AGN) in the Teacup galaxy at  ∼ 0.1 show that the gas density at distances of a few kpc away from the galaxy center is ≲ 10 cm −3 (Venturi et al. 2023).Therefore, we adopt a range of 1-40 cm −3 for the cool phase gas when calculating  turb to account for this wide range of uncertainty.For the hot phase with  ≈ 10 6 -10 7 K, we adopt a density range of 0.01-1 cm −3 (e.g., Li et al. 2018).We obtain an estimated  turb of ≈ 10 −26 -10 −22 erg cm −3 s −1 for the cool gas and ≈ 10 −28 -10 −25 erg cm −3 s −1 for the hot gas, as shown by the blue and red shaded regions in Figure 4. Ganguly et al. (2023) constrained the  turb of the ICM in the core regions of nearby cool-core clusters to be ≈ 10 −28 -10 −25 erg cm −3 s −1 (the gray shaded region in Figure 4), in agreement with our result for the hot phase.For star-forming molecular clouds, measurements across a wide range of spatial scales of ≈ 0.01-100 pc led to an estimate of  turb ≈ 10 −27 -10 −24 erg cm −3 s −1 as presented in Hennebelle & Falgarone (2012) and shown by the brown shaded region in Figure 4. Marchal et al. (2021) measured the density and kinematics of a bright concentration region near the edge of Complex C, an HVC in the Milky Way, which resulted in an estimated  turb ≈ 10 −30 -10 −28 erg cm −3 s −1 as shown by the green shaded region in Figure 4. Using non-thermal velocity widths of resolved absorption profiles and clump sizes inferred from photoionization models, Chen et al. (2023a) constrained  turb to be ≈ 10 −30 -10 −27 erg cm −3 s −1 for spectrally resolved cool clumps with a size scale of ≈ 10 pc -1 kpc in the CGM.These are shown by the gray data points in Figure 4.
It can be seen that the turbulent heating rates in the QSO nebulae, the cool-core cluster ICM, and the star-forming molecular clouds are on average ∼ 1000 times higher than that in the MW HVC and cool gas clumps probed in absorption.Given that both Complex C and cool absorption clumps are expected to be in relatively quiescent, undisturbed environments (Chen et al. 2023a), a possible explanation for this difference is that feedback due to star formation and AGN activities can significantly elevate the turbulent energy in the gaseous halos.However, caveats remain in this interpretation.As discussed in the previous section, the galaxy environments of the largest extended nebulae hint towards the scenario where tidal/merger interactions play a key role in stirring up the gas and facilitating the formation of multiphase structures, and the presence of a large amount of cool gas near the QSOs can lead to more efficient black hole accretion (e.g., Prasad et al. 2015;Voit et al. 2017).In this case, the elevated turbulent energy might be a precursor for fueling these luminous QSOs instead of a consequence of QSO feedback.
For the first time, we are able to determine turbulent energy transfer rate in the diffuse cosmic gas over seven decades in spatial scale from ∼ 0.01 pc to ∼ 100 kpc, but the measurements rely on two distinct approaches at different spatial scales.In particular, in the circumgalactic space, where we see three orders of magnitude difference in  turb from large to small scales, such distinction is also accompanied by differences in the way turbulence energy is determined.The gas turbulence probed in emission likely reflects the relative motions between different line-emitting clumps that trace the hot gas dynamics (as discussed in Section 4.1), while highresolution absorption line studies likely probe turbulence internal to individual clouds.Therefore, the lack of overlapping spatial scales probed by emission and absorption prevents us from forming a consistent picture of turbulent energy cascade in galaxy halos, while systematic uncertainties remain when comparing turbulent flows based on VSF measurements and those from absorption-line analyses.In Paper I, we discussed uncertainties associated with VSF measurements due to either projection effects (see also e.g., von Hoerner 1951; Xu 2020) or PSF smoothing (see further discussion in § 4.5).While the smallest area accessible in emission measures is limited to the PSF size of the data, the absorption line technique averages cloud properties over the beam size that is dictated by the black hole accretion disk size (i.e., on the order of ≪ 1 pc).At the same time, absorption-line analyses are subject to uncertainties in the photo-ionizing background radiation field.Future observations using AO-assisted ground-based IFSs and/or space-based IFSs can extend the small scales probed in the VSFs to ≲ 10 kpc for the line-emitting gas, bridging the gap in spatial scales accessible between emission and absorption studies.A sample of systems with both extended line emission and high-resolution absorption line data will also greatly aid in the investigation of this discrepancy in  turb .Finally, we note that these measurements of turbulent motions in QSO nebulae imply that turbulence is insufficient in providing the required energy to offset cooling at tens of kpc scales in the QSO environments.In Section 5.1 of Paper I, we compared the turbulent heating rate and the radiative cooling rate, utilizing measurements from TXS0206−048.The calculations considered the gas mass within a 50 kpc radius of a 5×10 13 M ⊙ halo, assuming an NFW mass profile and that gas of all phases is perfectly coupled dynamically.This approximate evaluation shows that the turbulent heating rate is on par with the luminosity of [O II] or [O III], yet it constitutes only approximately 0.05% of the QSO bolometric luminosity.

Velocity dispersion along the line of sight versus in the plane of the sky
In Figure 5, we show the velocity dispersion in the plane of the sky,  pos , versus the mean velocity dispersion along the line of sight, ⟨ los ⟩.  pos is quantified as the standard deviation of the line-of-sight velocity from spaxels included in the VSF measurements (see discussion in Section 2.3 and the velocity maps in Figures 9-16), and ⟨ los ⟩ is the mean line width (obtained through a single-component Gaussian fit) for the same set of spaxels.We show results for both [O II]  3727, 3729 and [O III]  5008 emission as they can differ in  pos and ⟨ los ⟩ due to the different footprints of the two lines.The statistical uncertainties of both velocity dispersions estimated through Monte Carlo resampling are small and are not shown in Figure 5.The measurements for  pos before and after removing the uni-directional velocity gradient in the plane of the sky are consistent with each other to within ≈ 20 km/s.Therefore, for clarity, we only show the values obtained using the directly measured [O II] and [O III] velocity maps.
It can be seen that for all nebulae in our sample,  pos ≲ ⟨ los ⟩.This observation agrees with the general trend seen in spatially-resolved data for H II regions where the velocity dispersion along the line of sight exceeds the velocity dispersion in the plane of the sky (e.g., Lagrois & Joncas 2011;Arthur et al. 2016;García-Vázquez et al. 2023).One possible explanation for this trend is the smoothing effect due to multiple line-emitting clouds along the line of sight contributing to the observed velocity centroid, leading to reduced velocity dispersion in the plane of the sky.In addition, the contribution from bulk/coherent motions along the line of sight will also result in larger  los .To investigate this possibility, we adopt a simple assumption that ⟨ los ⟩ 2 = [ 2 pos + ( grad,los ×  los ) 2 ], where  grad,los is the velocity gradient along the line of sight.We approximate the depth of the nebula  los to be the square root of the nebula size (see Table 3), and derive a velocity gradient of  grad,los ≈ 0.5-3 km/s/kpc for different nebulae.The range of this derived  grad,los is in qualitative agreement with the best-fitting velocity gradient in the plane of the sky (see Table 4), suggesting that bulk flows along the line of sight may be non-negligible.In contrast, the velocity dispersion across the plane of the sky provides a robust tracer of the underlying velocity variance at scales ≳ 10 kpc, particularly when a credible model for the coherent shear in the plane of the sky can be obtained with the spatially-resolved velocity measurements, as pointed out by previous studies (e.g., Stewart & Federrath 2022;García-Vázquez et al. 2023).

Power-law turnover scale for the VSFs
As discussed in Section 3.2, the shapes of the VSFs generally do not follow a single power-law across the entire range of scales probed.While additional structures in the VSFs may provide hints for different physical processes present in the nebulae, we caution that the limited nebula size and signalto-noise can hinder a robust interpretation of these structures.
In particular, we note that there is a moderate correlation (with a Spearman's  coefficient of 0.7) between the VSF turnover scale  2 (see Section 3.2) and the size of the nebula for both the [O II]  3727, 3729 and [O III]  5008 emission, as shown in Figure 6.This correlation indicates that the deviation of the VSF from a single power-law at larger scales is in part due to the limited nebula size probed by the data given the detection limit.Previous studies have also shown that boundaries of clouds/nebulae can artificially flatten the VSFs at large scales that mimic the signature of energy injection and affect the interpretation of the data (e.g., Ganguly et al. 2023;García-Vázquez et al. 2023).In addition, the smooth transition between the inertial range and the energy injection scale can cause the VSF slopes to taper off at a scale as small as half of the true energy injection scale (Federrath et al. 2021) and further complicate the interpretation of a flattening signal in the VSFs.
Given the abovementioned caveats, we refrain from interpreting  2 or VSF flattening scales in our sample as indicative of energy injection scales.However, Figure 7 indicates no discernible correlations between the constrained 2nd-order power-law slopes ( 2 ) and VSF turnover scale ( 2 ) or nebula size, underscoring the robustness of  2 measurements.Measurements from local H II regions reported by García-Vázquez et al. (2023) result in larger  2 values on average (shown in the blue shaded region in Figure 7), suggesting elevated Mach numbers in local H II regions and/or increased susceptibility to projection smoothing in their observations (for more discussion of projection effects see Section 4.5 below).

Limitations and caveats
A notable limitation in the present study arises from the projection effect inherent in the data.Several studies have investigated how VSFs are affected by the use of projected measurements.Analytically, von Hoerner (1951) derived that for volume-filling gas, the projection effect depends on the spatial scales probed: VSFs are steepened when measuring separation scales smaller than the depth of the cloud along the line of sight, while the VSF slopes recover to the intrinsic value at scales exceeding the cloud depth.This result is sometimes referred to as the "projection smoothing" effect and was independently confirmed by O'dell & Castaneda (1987) and Xu (2020).On the other hand, Zhang et al.  2023), which are on average higher than the slopes constrained for QSO nebulae.No discernible correlations are found between the 2nd-order power-law slopes ( 2 ) and VSF turnover scale ( 2 ) or nebulae size.
(2022) used numerical simulations to show that for spatially confined structures (e.g., isolated filaments), the projection effect flattens the VSFs.As we have discussed in Section 4.1, the dynamical state of the nebulae examined in this work indicates that the cool line-emitting gas is embedded in the hot ambient medium and traces the turbulent motions of the hot, volume-filling gas.Therefore, our measurements are more likely affected by the "projection smoothing" effect, suggesting that the intrinsic VSF slopes may be flatter than the values reported in Table 4, which still supports our interpretation of the subsonic/transonic gas motions.In addition to whether the gas is volume-filling or spatially confined, in reality, the projection effect will also depend on detailed properties of the system such as density/emissivity fluctuations and the three-dimensional geometry of the gas structure.Detailed investigations using high-resolution numerical simulations are needed to robustly quantify and calibrate the projection effect in more realistic environments.
Another main limitation of the current study is the restricted dynamic range in the VSF measurements, which is confined to approximately one decade or less in projected distance separation.This restriction prevented us from obtaining robust constraints on the VSFs slopes for several systems in our sample.While the largest separation is determined by the nebula size given the detection threshold, the smallest separation accessible in the data is dictated by the spatial sampling (i.e., angular size per spatial pixel) as well as the PSF size.As ground-based observations without adaptive optics (AO) are fundamentally limited by atmospheric seeing, improving the dynamic range towards small scales requires conducting AO-assisted observations on the ground (e.g., with VLT/ERIS in the infrared and using the Narrow-Field-Mode on VLT/MUSE in the optical) with longer exposure times to reach sufficient signal-to-noise.Alternatively, space-based IFSs such as JWST/NIRSpec with unprecedented spatial resolution have also started delivering an increasing sample of spatially-resolved observations of the CGM (e.g., Wylezalek et al. 2022;Veilleux et al. 2023).Finally, with a fixed PSF size, targeting systems at lower redshifts with a higher angular-to-physical size ratio can also help increase the VSF dynamic range.However, few extended (≳ 50 kpc) nebulae have been discovered at  < 0.5 (e.g., Rupke et al. 2019;Chen et al. 2019;Venturi et al. 2023) and additional effort is required to expand the sample size of low-redshift extended nebulae.

CONCLUSION
This paper presents an ensemble study of the turbulent motions in eight extended nebulae surrounding seven QSOs at  ≈ 0.5-1.1.Using the [O II]  3727, 3729 and/or [O III]  5008 emission lines, we measure the line-of-sight velocity fields and construct the velocity structure functions (VSFs).We probed the dynamical state of the gas illuminated by the QSO radiation field at scales ≈ 10-100 kpc.Our main conclusions are: • Five out of the eight nebulae in our sample have a constrained power-law slope of the 2nd-order VSFs,  2 , between ≈ 0.3-1.1, while the other three nebulae have loose constraints corresponding to 95% upper limits of ≲ 0.5-1.5, as shown in Figures 2 and 7 and discussed in Section 3.2.To within the 2- measurement uncertainty, the slopes are either consistent with the expectation from Kolmogorov turbulence or flatter, suggesting that the gas motions are subsonic.
• Removing a best-fitting unidirectional velocity gradient from the line-of-sight velocity maps flattens the VSFs in general, but also leads to larger uncertainties due to a reduced dynamic range in the VSFs that can be used for a single power-law fit.The results before and after removing a velocity gradient are consistent within the range of the uncertainty, as shown in Figure 2.
• Complementing the measurements for the 2nd-order VSF slopes,  2 , the ESS slope ratios   / 3 for  = 1-6 are also in agreement with the expectation of subsonic turbulence, as shown in Figure 3 and discussed in Section 3.3.The only exception is the nebula surrounding the QSO field HE0238−1904, where the   / 3 ratios are consistent with the supersonic MHD turbulence prediction by Boldyrev (2002) both before and after removing a uni-directional gradient field.A more detailed investigation of this field and a larger sample size are required to shed light on whether this field is a special case.
• The subsonic motions in the QSO nebulae suggest that the line-emitting cool clouds with  ∼ 10 4 K are embedded within a hot ambient medium with  ∼ 10 6 -10 7 K. Adopting the sound speed of the hot medium of  s,hot ≈500 km/s, we estimate the Mach number of the cool clouds to be ≈ 0.2-0.5, consistent with the observed VSF properties.The subsonic nature of gas motions supports a scenario where the cool clumps condense out of the hot gas, carrying the turbulent memory of the hot halo and serving as tracers of hot phase dynamics (see Section 4.1).
• No discernible differences are seen in VSF properties between radio-loud and radio-quiet QSO fields, suggesting that the collimated jets and their inflated bubbles do not play a critical role in shaping the dynamical state of the gas on ∼tens of kpc scales.
• Comparing the mean velocity dispersion along the line of sight, ⟨ los ⟩, and the velocity dispersion observed in the plane of the sky,  pos , we find that ⟨ los ⟩ ≳  pos for all fields (Figure 5).We discuss that projection effects and bulk motion along the line of sight are possible sources for the larger dispersion (see Section 4.3).
• The turbulent heating rate per unit volume,  turb , in the QSO nebulae is estimated to be ∼ 10 −26 -10 −22 erg cm −3 s −1 for the cool phase and ∼ 10 −28 -10 −25 erg cm −3 s −1 for the hot phase at scales ≈ 10-60 kpc.This range is in agreement with the measurements for intracluster medium and star-forming molecular clouds but is ∼ 1000 times higher than that estimated for Milky Way Complex C and cool circumgalactic gas clumps probed in low-ion absorption lines, as shown in Figure 4 and discussed in Section 4.2.While the difference in  turb might be a signpost for AGN/stellar feedback, a robust investigation into the systematics of the different measurements is required to shed light on this discrepancy.
Future observations of extended nebulae using AO-assisted IFSs on the ground (e.g., MUSE Narrow-Field-Mode) and/or space-based IFSs (e.g., JWST/NIRSpec IFU) will help extend the small scales probed in VSFs to ≲ 10 kpc, improving the robustness of the VSF constraints and bridging the gap between  turb measured by emission and absorption techniques.The findings of this ensemble study align with the recent emerging picture of the multiphase CGM where different gas phases are intricately connected throughout their formation and evolution history.For local H II regions as well as H II galaxies (at both low and high redshifts), a - relation corresponding to the correlation between the luminosity of the region/galaxy in a certain line emission (such as H and H) and its velocity dispersion is commonly observed (e.g., Melnick et al. 1987;González-Morán et al. 2021).In our QSO nebulae sample, however, we do not observe such a correlation, as shown in Figure 8.The contrast here likely arises from the different emission mechanisms for recombination lines versus collisionally-excited lines, the former of which is more well coupled to the total mass and stellar feedback in the H II regions/galaxies.In addition, QSOs are variable and the number of ionization photons output by QSOs is subject to significant changes on timescales of ≲tens of Myr (e.g., Schawinski et al. 2015;Sun et al. 2017;Shen 2021), further weakening a correlation between the luminosity of the surrounding nebulae and the velocity dispersion of the gas.

Figure 1 .
Figure 1.Continuum-and QSO-subtracted narrow-band images of the [O II] and [O III] emission from the three fields studied in this paper, based on the MUSE-WFM observations.The fields are shown in order of increasing redshift from left to right.Contours are at surface brightness levels of [5, 10, 50, 100] × 10 −18 erg s −1 cm −2 arcsec −2 .The yellow cross in each panel marks the quasar position.For PKS0405−123, we indicate the eastern nebula (Neb.E.) and the southern nebula (Neb.S.) with dashed boxes (same for both [O II] and [O III] emission), as these two nebulae are treated as separate systems despite originating from the same QSO field (see text for details).The narrow-band images for PKS0454−22, J0454−6116, J2135−5316 and TXS0206−048 were presented in Paper I.

Figure 2 .
Figure 2. Top row:The observed 2nd-order VSF  ′ 2 () for the nebulae of radio loud QSOs.Vertical dashed lines mark the size of the total PSF FWHM for each field.The data points and the error bars show the mean and the standard deviation for the 1000 measurements obtained through the modified bootstrap method (see Section 2).The dashed and solid gray lines show, respectively, the expected  2 for Kolmogorov turbulence before and after convolving with an appropriate PSF.As different fields have slightly different PSF sizes, we use the mean value of the PSF FWHM for all radio-loud fields included in the panels when constructing the expected Kolmogorov VSF.The dotted gray lines show a power-law with a slope of 1 (i.e., Burger's turbulence).The four panels from left to right show the results using the [O II] and [O III] lines, both from direct measurements and after removing a uni-directional gradient (see text).Bottom row: Same as the top row but for radio-quiet fields.All nebulae exhibit VSFs with an increasing amplitude in velocity variance as a function of separation distance.PSF smoothing significantly steepens the apparent slopes of the VSFs and we will explicitly take into account this smoothing effect when fitting the VSFs with a power-law (see Section 3.2).

Figure 4 .
Figure 4. Turbulent heating rate  turb at different scales for different physical systems.The red and blue shaded regions show the estimated  turb for QSO nebulae at scales of ≈ 10-60 kpc based on our sample.The calculations for the hot and cool gas phases assume densities of 0.01-1 cm −3 and 1-40 cm −3 , respectively (see text).The lower and upper bounds indicate the 16 th -84 th percentile ranges for measurements across all eight nebulae.Measurements from Ganguly et al. (2023) for ICM at scales ≈ 0.3-10 kpc are shown by the gray shaded region.Results for star-forming molecular clouds at scales ≈ 0.01-100 pc presented in Hennebelle & Falgarone (2012) are shown by the brown shaded region.Marchal et al. (2021) measured  turb for a bright concentration location in the HVC Complex C at scales ≈ 6-28 pc, as shown by the green shaded region.The gray points show the results for CGM cool clumps at scales ≈ 10 pc -1 kpc based on absorption line measurements presented inChen et al. (2023a).The turbulent heating rates in the QSO nebulae, the cool-core cluster ICM, and the star-forming molecular clouds are on average ∼ 1000 times higher than that in Complex C and cool gas clumps probed in absorption.

Figure 5 .
Figure 5.The positional velocity dispersion in the plane of the sky,  pos , versus the mean velocity dispersion along the line of sight, ⟨ los ⟩.The left and right panels show measurements using the [O II]  3727, 3729 and [O III]  5008 emission lines, respectively.The dashed line shows the relation  pos = ⟨ los ⟩, and the dotted line indicates the relation where ⟨ los ⟩ is twice the value of  pos .For all nebulae in our sample,  pos ≲ ⟨ los ⟩.

Figure 6 .
Figure 6.Nebula size versus the turnover scale  2 for all eight nebulae in the sample.The left and right panels show the values using the [O II]  3727, 3729 and [O III]  5008 emission lines, respectively.There is a moderate correlation between nebula size and the VSF turnover scale.

Figure 7 .
Figure 7. Top row: Nebula size versus the 2nd-order VSF slope  2 .Bottom row: The VSF turnover scale  2 versus the 2ndorder VSF slope  2 .The results based on the [O II]  3727, 3729 measurements are shown in the left panels while the results from [O III]  5008 are shown in the right panels.The horizontal blue shaded regions mark the measurements for local H II regions presented in García-Vázquez et al. (2023), which are on average higher than the slopes constrained for QSO nebulae.No discernible correlations are found between the 2nd-order power-law slopes ( 2 ) and VSF turnover scale ( 2 ) or nebulae size.

LFigure 8 .
Figure 8.The nebula emission line luminosity (for both [O II]  3727, 3729 and [O III]  5008 ) versus the velocity dispersion, both along the line of sight and in the plane of the sky.We do not observe a  −  correlation using the [O II] and [O III] extended emission surrounding QSOs.
B. MEASUREMENTS FOR INDIVIDUAL NEBULAEHere we present the VSFs measurements for individual nebulae in PKS0405−123, PKS0552−640, and HE0238−1904.The measurements for PKS0454−22, J0454−6116, and J2135−5316 can be found in the Appendix of Paper I.

Table 2
lists the coordinates, exposure time, and atmospheric seeing conditions of our sample.Out of the seven QSO fields, the measurements for four fields (PKS0454−22, J0454−6116, J2135−5316, and TXS0206−048) were presented in Paper I. The three newly included fields

Table 1 .
Summary of the QSO properties.
Number of spectroscopically-identified group member galaxies, including the QSO host. Velocity dispersion of the group.

Table 3 .
Summary of emission properties in QSO nebulae  .
3727, 3729 andBrandenburg & Lazarian 2013;Grete et al. 2021;les are consistent with the Kolmogorov slope.For TXS0206−048, only measurements with [O II]  3727, 3729 are available and the result is consistent with  2 = 2/3.While the VSF slope for the nebula PKS0405−123 E based on [O III]  5008 is flatter than 2/3 within the 16 th -84 th percentiles, the values within the 3 th -97 th percentiles using both [O II]  3727, 3729 and [O III]  5008 emission are in agreement with the Kolmogorov slope and therefore we consider the VSFs of this nebula consistent with the Kolmogorov expectation.Nebulae with  2 < 2/3: for the three nebulae in PKS 0454−22, J0454−6116 and J2135−5316, only upper limits of  2 can be obtained and the 95% limits derived from [O II]  3727, 3729 measurements are below 2/3.While the  2 upper limits obtained from the [O III]  5008 measurements are larger than 2/3, the smaller  2 upper limits obtained using [O II]  3727, 3729 suggest that the VSF slopes for these three nebulae are likely flatter than the Kolmogorov expectation.As we discussed in Paper I, the flatter VSFs may indicate the presence of multiple energy injection scales (e.g, ZuHone et al. 2016) and/or the effect of a dynamically important magnetic field (e.g.,Boldyrev 2006;Brandenburg & Lazarian 2013;Grete et al. 2021; Mohapatra et al. 2022).Nebulae with  2 > 2/3: Based on the directly measured velocity fields, PKS0405−123 S exhibits VSFs that are steeper than the expectation of Kolmogorov turbulence.The constraints are consistent using the [O II]  3727, 3729 and [O III]  5008 measurements.One possible explanation for

Table 4 .
Summary of the power-law slopes of the VSFs constructed using [O II] and [O III] lines  .
Ratios   / 3 for all eight nebulae in the sample based on both [O II]  3727, 3729 and [O III] Turbulence plays a critical role in facilitating nonlinear interactions within the gaseous halos, which in turn promote further developments of turbulence.For shaping the dynamical properties of gas traced by [O II] and [O III] at scales ≳ 10 kpc, environmental effects (e.g., tidal interactions, galaxy mergers, gas accretion) may dominate over QSO feedback.These findings can be directly compared with high-resolution numerical simulations to shed light on detailed physical mechanisms that govern the driving and development of turbulence in the CGM.We thank Fausto Cattaneo and Jenny Greene for helpful discussions throughout this work.We also thank Yuan Li for constructive feedback on the discussions of this paper.HWC and MCC acknowledge partial support from NSF AST-1715692 grants.ZQ acknowledges partial support from NASA ADAP grant 80NSSC22K0481.JIL is supported by the Eric and Wendy Schmidt AI in Science Postdoctoral Fellowship, a Schmidt Futures program.SC gratefully acknowledges support from the European Research Council (ERC) under the European Union's Horizon 2020 Research and Innovation programme grant agreement No 864361.FSZ acknowledges the support of a Carnegie Fellowship from the Observatories of the Carnegie Institution for Science.EB acknowledges support by NASA under award number 80GSFC21M0002.This research has made use of the services of the ESO Science Archive Facility and the Astrophysics Data Service (ADS)1.The analysis in this work was greatly facilitated by the following python packages: Numpy (Oliphant 2015), Scipy (Virtanen et al. 2020), Astropy (Astropy Collaboration et al. 2013; Price-Whelan et al. 2018), Matplotlib (Hunter 2007), and MPDAF (Bacon et al. 2016).This work was completed with resources provided by the University of Chicago Research Computing Center.