Effects of cloud geometry and metallicity on shattering and coagulation of cold gas, and implications for cold streams penetrating virial shocks Open Access

ABSTRACT

Theory and observations reveal that the circumgalactic medium (CGM) and the cosmic web at high redshifts are multiphase, with small clouds of cold gas embedded in a hot, diffuse medium. We study the ‘shattering’ of large, thermally unstable clouds into tiny cloudlets of size |$\ell _{\rm shatter}\sim {\rm min}(c_{\rm s}t_{\rm cool})$| using idealized numerical simulations. We expand upon previous works by exploring the effects of cloud geometry (spheres, streams, and sheets), metallicity, and an ionizing ultraviolet background. We find that ‘shattering’ is mainly triggered by clouds losing sonic contact and rapidly imploding, leading to a reflected shock that causes the cloud to re-expand and induces Richtmyer–Meshkov instabilities at its interface. The fragmented cloudlets experience a drag force from the surrounding hot gas, leading to recoagulation into larger clouds. We distinguish between ‘fast’ and ‘slow’ coagulation regimes. Sheets are always in the ‘fast’ coagulation regime, while streams and spheres transition to ‘slow’ coagulation above a critical overdensity, which is smallest for spheres. Surprisingly, |$\ell _\mathrm{shatter}$| does not appear to be a characteristic clump size even if it is well resolved. Rather, fragmentation continues until the grid scale with a mass distribution of |$N(\gt m)\propto m^{-1}$|⁠. We apply our results to cold streams feeding massive (⁠|$M_{\rm v}\lower.5ex\rm{\,\, \buildrel\gt \over \sim \,\,}10^{12}\, {\rm M}_\odot$|⁠) galaxies at |$z\lower.5ex\rm{\,\, \buildrel\gt \over \sim \,\,}2$| from the cosmic web, finding that streams likely shatter upon entering the hot CGM through the virial shock. This could explain the large clumping factors and covering fractions of cold gas around such galaxies, and may be related to galaxy quenching by preventing cold streams from reaching the central galaxy.

1 INTRODUCTION

Only a small fraction of the Universe’s baryons are found in galaxies, including both stars and interstellar gas (e.g. Peeples et al. 2014; Tumlinson, Peeples & Werk 2017; Wechsler & Tinker 2018). The majority of baryons, and also the majority of metals, reside in the circumgalactic medium (CGM, gas outside galaxies but within dark matter, DM haloes), and the intergalactic medium (IGM, gas outside DM haloes). Besides their importance for the cosmic baryon budget, the physical properties and chemical composition of the C/IGM offer valuable insight into galaxy evolution, since they supply galaxies with fresh gas and also act as a reservoir for their ejected, enriched gas (the cosmic baryon cycle, e.g. Putman, Peek & Joung 2012; McQuinn 2016; Tumlinson et al. 2017). Moreover, the distribution of neutral hydrogen (H i) in the high-z IGM can be used to constrain cosmic reionization, structure formation, and the nature of DM through the Lyman-|$\alpha$| (Ly|$\alpha$|⁠) forest (e.g. Rauch 1998; Viel et al. 2013; Lidz & Malloy 2014; McQuinn 2016; Eilers, Davies & Hennawi 2018).

Gas in the C/IGM is highly diffuse and difficult to directly observe. It has traditionally been traced using absorption line spectroscopy along lines of sight to distant quasi-stellar objects (QSOs) or galaxies (e.g. Bergeron 1986; Hennawi et al. 2006; Steidel et al. 2010; Lehner et al. 2022). In recent years, the advent of new integral field unit spectographs such as Keck Cosmic Web Imager (KCWI) on Keck and Multi Unit Spectroscopic Explorer (MUSE) on the Very Large Telescope (VLT) have enabled emission-line studies of the CGM and IGM around galaxies at |$z\lower.5ex\rm{\,\, \buildrel\gt \over \sim \,\,}3$| (Steidel et al. 2000; Cantalupo et al. 2014; Martin et al. 2014a, b; Umehata et al. 2019). Both emission- and absorption-line studies reveal that the gas in and around galaxy haloes has a complex multiphase structure (Tumlinson et al. 2017). Surprisingly, large amounts of cold (⁠|$\sim 10^4 {\rm K}$|⁠) gas have been observed in the outskirts of galaxy haloes, which cannot be in hydrostatic equilibrium with the halo gravitational potential. During cosmic noon, at |$z\sim (2-6)$| near the peak of cosmic galaxy formation, this cold gas is inferred to be denser than the ambient hot gas within which it is embedded by a factor |$\chi \equiv \rho _{\rm c}/\rho _{\rm h}\sim 10^2-10^3$|⁠, and to be composed of tiny clouds with sizes of |$l\sim N_{\rm H}/n_{\rm H}\lower.5ex\rm{\,\, \buildrel\lt \over \sim \,\,}50\, {\rm pc}$| (Cantalupo et al. 2014; Hennawi et al. 2015; Borisova et al. 2016). The cold gas has order unity area covering fractions, |$f_{\rm C}\sim O(1)$|⁠, and mass fractions with respect to the hot CGM mass, |$f_{\rm M}\sim M_{\rm cold}/M_{\rm hot}\sim O(1)$| (Pezzulli & Cantalupo 2019). However, its volume filling factor is tiny, |$f_{\rm V}\sim f_{\rm M}/\chi \sim 10^{-3}$|⁠, making it extremely clumpy (Cantalupo et al. 2019) and its apparent abundance difficult to explain (Faucher-Giguère & Oh 2023).

Recent theoretical advances have shed new light on these issues. A major new insight (McCourt et al. 2018) is that when the cooling time of a gas cloud is much less than its sound crossing time such that it cannot cool isobarically, it does not cool isochorically as had been presumed (Field 1965; Burkert & Lin 2000). Rather, the cloud ‘shatters’ into many small fragments that lose sonic contact, causing them to contract independently and subsequently disperse, similar to a terrestrial fog (McCourt et al. 2018; Gronke & Oh 2020b). The typical size of the resulting cloudlets is expected to be of order the minimal cooling length, |$\ell _{\rm shatter} \sim {\rm min}(c_{\rm s} t_{\rm cool})$|⁠, with |$c_{\rm s}$| and |$t_{\rm cool}$| the sound speed and cooling time, and the minimal value obtained at |$T \sim 10^4\,{\rm K}$|⁠. For typical CGM conditions at |$z\sim (2-3)$|⁠, this is |$\lower.5ex\rm{\,\, \buildrel\gt \over \sim \,\,}10\, {\rm pc}$|⁠, consistent with inferred cloud sizes. This would explain the vastly different area covering and volume filling factors, since for N droplets of size l dispersed throughout a system of size R, |$f_{\rm C}/f_{\rm V}\sim R/l\gg 1$| (Faucher-Giguère & Oh 2023). A ‘fog’ can also explain a host of additional observations in the CGM, high velocity clouds, quasar broad-line regions, and the interstellar medium (ISM, Gronke et al. 2017; McCourt et al. 2018; Stanimirović & Zweibel 2018; Faucher-Giguère & Oh 2023; Sameer et al. 2024).

Despite the many appealing features of the shattering model, numerous puzzles remain. It is unclear under what conditions a large cooling cloud will shatter, with some suggesting this depends on the final overdensity between the cold and hot gas (Gronke & Oh 2020b) and others that it depends on the thermal stability conditions in the initial cloud (Waters & Proga 2019a; Das, Choudhury & Sharma 2021). Even when clouds do shatter in 3D simulations, they do not appear to do so hierarchically as was initially proposed by McCourt et al. (2018). Rather, if the initial cloud is large (⁠|$r_{\rm cl}\gg \ell _{\rm shater}$|⁠) and highly non-linear (⁠|$\delta \rho /\rho \gg 1$|⁠) when it loses sonic contact it initially cools isochorically, then becomes strongly compressed by its surroundings until its central pressure overshoots, and finally it explodes into many small fragments (Gronke & Oh 2020b). This process is sometimes referred to as ‘splattering’ (Waters & Proga 2019a, 2023), and seems to be due to vorticity generated by Richtmyer–Meshkov instabilities (RMI; Richtmyer 1960; Meshkov 1969; Zhou 2017a, b), which explains why it is not seen in 1D simulations (Waters & Proga 2019a; Das et al. 2021). Additional fragmentation mechanisms have been proposed, such as shredding by collisions with larger fragments (Jennings & Li 2021), or rapid rotation of clumps (‘splintering’; Farber & Gronke 2023). We will hereafter use the term ‘thermal fragmentation’ to refer to the general process of cold-gas fragmentation into numerous small clouds, as the process may be very different than that originally proposed by McCourt et al. (2018).

Even if a cloud initially fragments, the resulting cloudlets may recoagulate to form larger clouds. Gronke & Oh (2020b) found that a cloud remained fragmented if its final overdensity

where |$\rho _\mathrm{s,f}$|¹ is the cloud density at the temperature floor and |$\rho _\mathrm{bg}$| is the background density, was above a critical value of |$\chi _\mathrm{f,crit}\sim 300$| with a weak dependence on cloud size of |$(r_{\rm cl}/\ell _{\rm shatter})^{1/6}$|⁠. The origin of this threshold remains unclear. Several coagulation mechanisms have been discussed in the literature (see summary in Faucher-Giguère & Oh 2023). These include direct collisions, similar to dust grain growth in protoplanetary discs, and coagulation due to the advective flow generated by hot gas condensing onto a cold cloud. In the latter, the inflow velocity can be set by thermal conduction (or numerical diffusion; Elphick, Regev & Spiegel 1991; Elphick, Regev & Shaviv 1992; Koyama & Inutsuka 2004; Waters & Proga 2019b) or, more relevant for our purposes, by the so-called mixing velocities through turbulent mixing layers (Gronke & Oh 2023). Of particular interest is that in the case of turbulent mixing layers, the coagulation can be modelled as an effective force between two clouds, which scales as |$r^{-n}$|⁠, with r the distance between the clouds and |$n=2$|⁠, 1, or 0 for clouds in 3, 2, or 1 dimensions, similar to the gravitational force (Gronke & Oh 2023).

Studying this process in a cosmological context is extremely challenging, since even in the most advanced cosmological simulations employing novel methods to enhance the resolution in the CGM (Hummels et al. 2019; Peeples et al. 2019; Suresh et al. 2019; van de Voort et al. 2019) or the IGM (Mandelker et al. 2019b, 2021), cell sizes are still significantly larger than |$\ell _{\rm shatter}$|⁠. As a result, the amount and extent of cold, dense, low-ionization gas in the CGM increases with resolution and is not converged. Furthermore, different simulations disagree on the magnitude of the effect of enhanced CGM refinement, at least in part due to the different subgrid models employed for galaxy formation physics, such as stellar and active galactic nucleus (AGN) feedback, galactic winds, and gas photoheating and photoionization. This has obscured the details of why higher resolution leads to more cold gas in the CGM, where this cold gas comes from, and what a meaningful convergence criterion for the formation of multiphase gas might be. Numerically, it seems that the formation of multiphase gas requires resolving the cooling length at |$T\sim 10^5{\rm K}$| where isochoric cooling modes are stable (Das et al. 2021; Mandelker et al. 2021), though it is unclear how generic this is and other convergence criteria have been proposed (Hummels et al. 2019; Gronke et al. 2022).

For these reasons, thermal fragmentation and coagulation are most commonly studied using idealized simulations and numerical models. In the vast majority of cases, such models assume a spherical or quasi-spherical cloud or distribution of clouds. However, many systems in the C/IGM where these processes are important are filamentary (cylindrical) or planar in nature. Modern cosmological simulations reveal strong accretion shocks around intergalactic filaments (Ramsøy et al. 2021; Lu et al. 2024) and sheets (Mandelker et al. 2019b, 2021) that make up the ‘cosmic web’ of matter on Mpc-scales, similar to virial accretion shocks around massive DM haloes (Rees & Ostriker 1977; White & Rees 1978; Birnboim & Dekel 2003; Stern et al. 2021). The post-shock gas in the high-z cosmic web can fragment (Mandelker et al. 2019b, 2021; Lu et al. 2024), with the resulting cold cloudlets in intergalactic sheets potentially explaining observations of extremely metal-poor Lyman-limit systems (e.g. Robert et al. 2019; Lehner et al. 2022). The small-scale structure of cosmic filaments and sheets have additional important consequences for a wide variety of issues, including how gas is accreted onto DM haloes, interpretations of Ly|$\alpha$| forest statistics, measured dispersions of Fast radio bursts (FRBs), radiative transfer and the self-shielding of photoionized gas, and the overall cosmic census of baryons (see discussion in Mandelker et al. 2021). On smaller scales, gas accretion onto massive galaxies at high-z is thought to be dominated by cold streams flowing along cosmic-web filaments, which penetrate the halo virial shock and flow freely towards the central galaxy (e.g. Dekel & Birnboim 2006; Dekel et al. 2009). The interaction of these cold streams with the hot CGM can lead to fragmentation and the formation of small-scale cold clouds (Mandelker et al. 2020a; Lu et al. 2024). Finally, filamentary structures are expected around both inflowing and outflowing gas clouds in the CGM due to cloud–wind interactions (e.g. Banda-Barragán et al. 2016, 2019; Gronke & Oh 2018; Li et al. 2020; Sparre, Pfrommer & Ehlert 2020; Gronke & Oh 2020a; Tan, Oh & Gronke 2023; Tan & Fielding 2024).

An additional complication arises due to the metallicity of the gas, which affects the cooling rates and therefore |$\ell _{\rm shatter}$| and the resulting cloud sizes, as well as the strength of coagulation forces. While most studies in the literature assume solar metallicity gas (though see Das et al. 2021), gas in the high-z cosmic web is expected to have much lower metallicity, |$Z\sim (10^{-4}-0.1)Z_{\odot }$| (Mandelker et al. 2021). Similarly, the presence of a ultraviolet (UV) background (UVB) is typically not included in studies of thermal fragmentation although it too can influence cooling rates and cloud sizes.

In this paper, we explore the effects of geometry and metallicity on thermal fragmentation and coagulation using idealized 3D simulations of spherical clouds, cylindrical filaments, and planar sheets with metallicity values in the range |$(0.03-1.0)Z_{\odot }$|⁠. In Section 2, we introduce our numerical tools and simulation methods. In Section 3, we compare the evolution of fragmentation and coagulation processes among planar, cylindrical, and spherical geometries at solar metallicity. In Section 4, we extend the cylindrical geometry to lower metallicity and include a UVB. In Section 5, we discuss the size distribution of cloudlets formed through thermal fragmentation. We present a model for the coagulation criteria in both streams and spheres in Section 6, and tentatively apply this to the case of cold streams penetrating the halo virial shock in Section 7. In Section 8, we address caveats caused by additional physical processes not included in our analysis, before concluding in Section 9.

2 NUMERICAL METHODS

In this section, we describe the details of our simulation setup.

2.1 Simulation code

We use Eulerian adaptive mesh refinement (AMR) code ramses² (Teyssier 2002) to perform 3D idealized numerical simulations. We adopt the multidimensional MonCen limiter (van Leer 1977) for the piecewise linear reconstruction, the Harten–Lax–van Leer–Contact (HLLC) approximate Riemann hydro solver (Toro, Spruce & Speares 1994) for calculating fluxes at cell interfaces, and the Monotone Upstream-centered Schemes for Conservation Laws (MUSCL)–Hanchock scheme (van Leer 1984) for the Godunov integrator. The adiabatic index is |$\gamma =5/3$|⁠, and the Courant factor is 0.6.

2.2 Radiative cooling

We utilize the standard ramses cooling module, which accounts for atomic and fine-structure cooling for our assumed metallicity values. When comparing the three geometries in Section 3, we assume solar metallicity for both cold and hot phases and set a temperature floor at |$T_\mathrm{floor}=1.68\times 10^{4}\,$| K. When focusing on cold streams in Section 4, we follow Mandelker et al. (2020a) and assume metallicity values of |$Z_\mathrm{bg}=0.1\, Z_\odot$| for the background CGM and |$Z_\mathrm{s}=0.03\, Z_\odot$| for the streams, and include photoheating and photoionization from a |$z=2$| (Haardt & Madau 1996) UVB. At our assumed densities, the equilibrium temperature between the radiative cooling and the UV heating is roughly at |$T_\mathrm{floor}$| which we adopt in Section 3. In all cases, we shut-off cooling above |$0.8\, T_\mathrm{bg}$| to prevent the cooling of background (e.g. Gronke & Oh 2018; Mandelker et al. 2020a).

2.3 Clump finder

To quantify the degree of cold-gas fragmentation, we utilize the built-in ramses clump finder module phew (Bleuler et al. 2015), which is a parallel segmentation algorithm for 3D AMR data sets. It detects connected regions above a certain density threshold based on a ‘watershed’ segmentation of the computational volume into dense regions, followed by a merging of the segmented patches based on the saddle point density. Basically, each clump is centred on a local density peak (whose density is above the density threshold) and includes all surrounding gas with densities above both the density threshold and any local saddle points or density minima. Two neighbouring peaks separated by a saddle point whose density is above the secondary saddle density threshold are then ‘merged’ into a single clump. If the saddle density is below this secondary threshold, the two peaks represent two distinct clumps. We choose the clump density threshold to be the initial density of warm gas (see Section 2.5 below), while the saddle density threshold is the geometric mean of the initial warm gas density and the final cold-gas density.

2.4 Grid structure and boundary conditions

The coordinates can be generalized by |$(x_1,x_2,x_3)$|⁠, which represents |$(x,y,z)$|⁠, |$(r,\phi ,z)$|⁠, and |$(r,\theta ,\phi)$| in sheets, streams, and spheres, respectively. The sheets are aligned with the |$yz$| plane, while the stream axis is aligned with the z-axis. Sheets are initially confined to the region |$|x|\le r_\mathrm{s,i}$|⁠, while streams and spheres are initially confined to |$r\le r_\mathrm{s,i}$|⁠. Here, |$r_\mathrm{s,i}$| represents the initial sheet half-thickness, cylindrical radius, and spherical radius for the three different geometries, with cold gas always occupying the region |$|x_1|\le r_\mathrm{s,i}$|⁠.

The simulation region is a cubic box with size |$L=32\, r_\mathrm{s,i}$|⁠. We use a statically refined mesh with higher resolution around the cold gas. In our fiducial setup, the maximal refinement level is 10 corresponding to a minimal cell size |$\Delta =L/1024=r_{\rm s,i}/32$| valid in the region |$|x_1|\lt 3\, r_\mathrm{s,i}$|⁠. The cell size increases by a factor of 2 at |$|x_1|=[3,6,9,12]r_{\rm s,i}$|⁠, reaching a maximal value of |$r_{\rm s,i}/2$| at |$|x_1|\gt 12\, r_\mathrm{s,i}$|⁠. Our fiducial grid structure for stream geometry is illustrated in Fig. 1.

We adopt outflow conditions for all boundaries, such that the gradients of all hydrodynamic variables are set to zero. We note that while periodic boundary conditions along the stream axis and within the plane of the sheet would have been preferable to model the idealized cases of infinitely long streams and sheets, for technical reasons to do with our clump finder we were forced to adopt the same boundary conditions on all six boundaries. We opted to implement outflow boundary conditions everywhere as these are necessary to allow correct entrainment flows to develop perpendicular to the sheet plane and the stream axis, which are necessary for properly modelling coagulation. While this has no impact on our analysis of spheres, we find that streams contract along their axis and sheets within the plane due to coagulation along these axes induced by entrainment flows that develop after the initial fragmentation. To avoid any potential effects of these boundary conditions on our analysis in streams and sheets, we restrict our analysis of these geometries to a narrower box, excluding gas within |$10\, r_{\rm s,i}$| of the boundaries along the stream axis |$(x_3)$| and within the sheet plane |$(x2, x3)$|⁠. We find this narrower box to be unaffected by this contraction over the run time of our simulations, |$\sim (10-20)t_{\rm sc}$|⁠, where

is the sound crossing time of the initial cloud radius at the sound speed of cold gas, |$c_{\rm s,c}$|⁠, with a temperature of |$T_{\rm c}\sim T_{\rm floor}\lower.5ex\rm{\,\, \buildrel\gt \over \sim \,\,}10^4\,{\rm K}$|⁠.

2.5 Initial conditions

We initialize the simulations with a static, warm component (sheet/stream/sphere) of density |$\rho _\mathrm{s,i}$| in pressure equilibrium with a static, hot background of density |$\rho _\mathrm{bg}$|⁠. The initial overdensity is thus |$\chi _\mathrm{i}\equiv \rho _\mathrm{s,i}/\rho _\mathrm{bg}$|⁠. The temperature of the warm gas, |$T_\mathrm{s,i}$|⁠, is lower than that of the hot background, |$T_\mathrm{bg}$|⁠, but higher than the temperature floor |$T_\mathrm{floor}$|⁠. As the warm component cools towards thermal equilibrium at |$T_\mathrm{floor}$|⁠, a pressure contrast is generated between the cold gas and the hot background. We define |$\eta \equiv \rho _\mathrm{s,f}/\rho _\mathrm{s,i}=(T_\mathrm{s,i}/\mu _\mathrm{s,i})/(T_\mathrm{floor}/\mu _\mathrm{s,f})$|⁠, where |$\rho _\mathrm{s,f}$| denotes the final density of cold gas once pressure equilibrium has been re-established at |$T_\mathrm{floor}$|⁠, and |$\mu _{\rm s,i}$| and |$\mu _{\rm s,f}$| are the mean molecular weight in the initial warm and final cold gas, respectively. Consequently the final overdensity between the cold and hot gas is |$\chi _\mathrm{f}\equiv \rho _\mathrm{s,f}/\rho _\mathrm{bg}=\eta \chi _\mathrm{i}$|⁠, which is expected to be the key parameter determining whether cold gas remains fragmented or coagulates (Gronke & Oh 2020b). We fix |$\rho _\mathrm{s,f}=0.01\, \mathrm{m_p}$| for all simulations and vary |$\chi _\mathrm{f}$| by changing |$\rho _\mathrm{bg}$|⁠.

Table 1 summarizes the parameters of the simulations analysed throughout this work, and the sections where they are discussed.

2.6 Perturbations

We introduce density perturbations in the initial warm gas component. In units of the initial mean density of warm gas, the density follows a Gaussian distribution with |$(\mu ,\sigma)=(1,0.01)$|⁠, truncated at |$3\, \sigma$|⁠, similar to Gronke & Oh (2020b). We further introduce shape perturbations at the interface between the warm and hot components, as described below. Such interface perturbations have been shown to suppress the carbuncle instability (Quirk 1994), which is a numerical instability affecting strong shock fronts on the grid scale in multidimensional simulations, by misaligning the interface and the grid, and by generating additional vorticity and turbulence when the shock overtakes the interface.

For the stream geometry, we adopt the same interface perturbations as implemented in Mandelker et al. (2019a, 2020a),

Here, |$\delta r=0.1\, r_\mathrm{s,i}$| is the rms amplitude of perturbations, and |$k_j=2\pi n_j$| where |$n_j$| is an integer corresponding to a wavelength |$\lambda _j=1/n_j$|⁠. We include all wavenumbers in the range |$n_j=16-64$|⁠, corresponding to wavelengths in the range from |$r_\mathrm{s,i}/2$| to |$2\, r_\mathrm{s,i}$|⁠. |$m_j$| is the azimuthal mode number of the perturbation, with |$m=0$| corresponding to axisymmetric modes, |$m=1$| to helical modes, and |$m\ge 2$| to high-order fluting modes (Mandelker et al. 2016, 2019a). For each longitudinal wavenumber, |$n_j$|⁠, we include two azimuthal modes, |$m_j=0,1$|⁠. This results in a total of |$N_\mathrm{pert,str}=2\times 49=98$| modes. Each mode has a random phase, |$\phi _j \in [0,2\pi)$|⁠.

For sheets, we use analogous interface perturbations

where |$k_{x,j}=2\pi n_{x,j}$|⁠, |$k_{y,j}=2\pi n_{y,j}$|⁠. |$n_{x,j}$|⁠, and |$n_{y,j}$| are the wavenumbers along the x- and y-directions, chosen from 16 to 64 with an interval of 6, i.e. |$[16,\, 22,\, ...,\, 58,\, 64]$|⁠. This corresponds a total of |$N_\mathrm{pert,sht}=9\times 9=81$| modes, each with a random phase |$\phi _j \in [0,2\pi)$|⁠.

For spheres, we use spherical harmonics to perturb the interface:

where |$Y_l^{m}$| is the spherical harmonic given by

with |$P_l^{m}$| the associated Legendre polynomial. We choose |$l_{\rm max}=13$|⁠, yielding |$N_\mathrm{pert,sph}=1/2\times 14\times 15=105$| modes.

3 GEOMETRICAL EFFECTS ON THERMAL FRAGMENTATION AND COAGULATION

In this section, we present our results comparing thermal fragmentation and coagulation in the three geometries at solar metallicity. We begin in Section 3.1 by addressing the initial implosion and explosion phases of the fragmentation process. Then, in Section 3.2, we discuss the subsequent coagulation of the resulting cloudlets.

3.1 Cold-gas fragmentation

3.1.1 The implosion and explosion of cold gas

We begin with a detailed physical description of the implosion of the initial cloud, triggered by a lack of pressure support due to cooling, and the subsequent explosion triggered by a reverse shock reflected off the cloud centre. While this process is interesting in its own right, the main outcome relevant for our discussion of coagulation which follows is the ‘explosion velocity’, |$v_{\rm ex}$|⁠, the characteristic velocity of the cold gas once the reverse shock reaches scales of order the final cloud radius. In Appendix  A, we present a detailed mathematical discussion of these processes and a derivation of |$v_{\rm ex}$| for sheets. In this section, we offer more general considerations valid for all three geometries, and demonstrate these using results from one simulation for each geometry (Fig. 2).

Consider a warm gas cloud with cooling time,

with n the gas number density, k the Boltzmann constant, |$\gamma =5/3$| the adiabatic index of the gas, and |$\Lambda$| the cooling function. The sound crossing time is

with P the gas pressure. Mass conservation during the collapse tells us that |$r_{\rm s}\propto \rho ^{-1/(n+1)}$| with |$n=2$|⁠, 1, and 0 for spheres, streams, and sheets, respectively. Thus, |$t_{\rm cool}/t_{\rm sc}\propto \rho ^{-(3n+1)/(2n+2)}$|⁠, which decreases as the density rises. A cloud for which |$t_{\rm cool}\gt t_{\rm sc}$| initially will cool isobarically, contracting, and growing denser as it cools, until it reaches a radius |$r^{*}$| where |$t_{\rm cool}$| becomes shorter than |$t_{\rm sc}$|⁠. At this stage, the cloud loses sonic contact and proceeds to cool isochorically (Burkert & Lin 2000; Gronke & Oh 2020b; Waters & Proga 2023). From this point, the pressure within the cloud drops as it cools, driving it further away from equilibrium with the surrounding background. The resulting pressure gradient causes the cold gas to implode, with an isothermal shock (due to strong cooling) propagating inwards and a rarefaction wave outwards. This implosion decelerates and eventually reverses when the shock nears the centre and the central pressure becomes large compared to the background pressure, causing the cold gas to explode outwards. Eventually, the cold gas regains pressure equilibrium with its surroundings such that |$\rho _\mathrm{s,f}T_\mathrm{floor}/\mu _\mathrm{s,f}=\rho _\mathrm{s,i}T_\mathrm{s,i}/\mu _\mathrm{s,i}$|⁠, or |$\rho _\mathrm{s,f}=\eta \rho _\mathrm{s,i}$| (Section 2.5). Assuming that most of the cold-gas mass in this final state is in a single cloud, mass conservation yields |$\rho _\mathrm{s,f}r_\mathrm{s,f}^{n+1}=\rho _\mathrm{s,i}r_\mathrm{s,i}^{n+1}$|⁠, with |$r_\mathrm{s,f}$| the final radius of the cold cloud. We thus obtain |$r_\mathrm{s,f}=r_\mathrm{s,i}\, \eta ^{-1/(n+1)}$|⁠, showing that at a given |$\eta$|⁠, the final contraction radius decreases from spheres to streams to sheets.

If the peak central pressure is sufficiently large compared to the background pressure (see Fig. 2, panel c3, discussed below), a strong reflected shock propagates outwards from the cold to the hot gas, and a rarefaction wave propagates inwards causing the cold gas to expand. The post-shock cold gas is accelerated outward by the shock, and the peak cold-gas expansion velocity occurs when the central pressure subsides and pressure equilibrium between the cold and hot gas has been regained. Below, we provide three explanations for why this peak velocity, |$v_\mathrm{ex}$|⁠, is of order the cold-gas sound speed, |$c_{\rm s,c}$|⁠.

We first consider energy conservation as the thermal energy after the implosion gets converted into kinetic energy of the expanding gas. The total internal energy after the implosion is

where |$m_\mathrm{s,i}$| is the initial mass in the cooling cloud, and |$T_{\rm c}$| is the final equilibrium cloud temperature.³ When the shock propagation time is shorter than the cooling time of post-shock gas, we can assume that most of this internal energy becomes the kinetic energy of cold gas, |$E\sim 1/2m_\mathrm{s,i}v_\mathrm{ex}^{2}$|⁠. We thus obtain

Note that once thermal and pressure equilibrium have been re-established at the very end of the fragmentation process, the internal energy of the cold gas is the same as it is after the implosion since the same mass of gas is at the same temperature, |$T_{\rm c}$|⁠. It would thus seem that one cannot convert the internal energy in equation (9) to kinetic energy. However, the peak pressure at the end of the implosion phase initially causes the cloud to adiabatically expand, with the central temperature dropping to as low as |$\lower.5ex\rm{\,\, \buildrel\lt \over \sim \,\,}0.3\, T_{\rm c}$| before rising back up to |$T_{\rm c}$| due to mixing and compression of the hot gas. This drop in temperature seen in our simulations cannot be due to cooling, since we do not allow cooling below |$T_{\rm c}$|⁠. While energy is not conserved throughout the implosion–explosion process due to strong cooling, the radiative losses primarily come at the expense of the kinetic energy of imploding gas, which declines without a corresponding increase in thermal energy due to strong cooling keeping the cold gas roughly isothermal during the implosion (Fig. 2, panels b2 and b3, and Fig. 3, dashed and dot–dashed blue lines). On the other hand, after the end of the implosion phase, the internal energy of post-implosion gas is converted to kinetic energy of exploding gas (Fig. 3). It is during this adiabatic expansion phase that the cold gas is accelerated to |$v_{\rm ex}$|⁠. The energy and velocity evolution of the cloud during the first |$\, {t_{\rm sc}}$| is illustrated in Fig. 3.

An alternative way to see that the explosion velocity must be of order the cold-gas sound speed is to consider the shock-jump conditions. For isothermal shocks, we know that |$v_1 v_2 \sim c_{\rm s,c}^2$|⁠, where |$v_1$| and |$v_2$| are the pre- and post-shock velocities in the shock frame. |$v_1$| must be several times larger than |$c_{\rm s,c}$| due to the shock speed exceeding |$c_{\rm s,c}$| and the negative velocity of the still imploding pre-shock cold gas. Thus, |$v_2$| must remain small, indicating that the post-shock gas should have a velocity close to |$c_{\rm s,c}$| in the lab frame.

Yet a third way to understand why |$v_{\rm ex}\sim c_{\rm s,c}$| is as follows. One can think of the contracting cloud as a spring which is compressed and then released. Thus, from energy conservation, the explosion velocity cannot exceed the implosion velocity (in general, it will be smaller, because of radiative losses). The implosion velocity, or velocity of the cloud-crushing shock |$v_{\rm s}$|⁠, is given by |$\rho _{\rm c} v_{\rm s}^2 \sim \delta P \sim P \sim \rho _{\rm h} c_{\rm s,h}^2$|⁠, or |$v_{\rm s} \sim c_{\rm s,h}/\sqrt{\chi } \sim c_{\rm s,c}$| (Klein, McKee & Colella 1994). Hence, the cold-gas velocity is a characteristic expansion velocity. In detail, the implosion velocity is somewhat larger than we have estimated (since the overdensity during contraction is less than |$\rho _{\rm c}/\rho _{\rm h}$|⁠), and the expansion velocity is somewhat smaller (due to radiative losses). However, the estimate |$v_{\rm expand} \sim c_{\rm s,c}$| is robust.

All these estimates suggests that |$v_\mathrm{ex}\sim c_\mathrm{s,c}$| for all geometries, regardless of the initial conditions (see Appendix  A for a detailed derivation of |$v_{\rm ex}$| for sheet geometries).

In Fig. 2, we demonstrate the key features of the implosion and explosion processes in each of our three geometries. From top to bottom, we show radial profiles of density, temperature, pressure, and radial velocity, taken from simulations with |$r_\mathrm{s,i}=3\, \mathrm{kpc}$|⁠, |$\eta =10$|⁠, and |$\chi _\mathrm{f}=100$| (first row of Table 1). The density and pressure profiles are volume-weighted averages within each radial bin, while the temperature and velocity profiles are mass-weighted. We show the profiles at four times, from left to right these are at the initial condition, during the implosion, at the central shock collision, and near the end of the explosion phase when pressure equilibrium has been re-established and the explosion velocity has reached its peak value. Different colour lines mark the different geometries, while black dashed and dotted lines show results for cold and hot gas respectively (separated at |$T=10^5\, {\rm K}$|⁠) in stream geometry. Separating cold and hot gas for spheres yields similar results as for streams and is not shown, while the results for sheets are presented in Appendix  A.

Initially, the cold gas is in pressure equilibrium with the hot background, with a density contrast of |$\chi _{\rm i}=10$|⁠, and the fluid velocity is zero everywhere. Note that the density and temperature profiles exhibit a smooth transition between the cloud and the background. This is due to the shape perturbations we include (Section 2.6) and not to any explicit smoothing or ramp function in the initial conditions.

The initial cloud properties yield |$t_\mathrm{cool}\lt t_\mathrm{sc}$|⁠, so the cold gas loses sonic contact as soon as it starts cooling, meaning |$r^{*}=r_\mathrm{s,i}$|⁠. Consequently, the central pressure rapidly drops, within a cooling time, forming a pressure gradient between the cold cloud and the hot background. This results in an isothermal shock propagating inwards and a rarefaction wave outwards. The shock is visible in panel (c2) of Fig. 2 as a jump in pressure. Note that the contact discontinuity (CD), where the gas temperature begins rising in panel (b2), is outside the implosion shock, while cold gas both inside and outside the shock has |$T\sim T_{\rm c}$|⁠. In sheets, the implosion shock propagates inwards with a roughly constant Mach number of |$\sim 2$|⁠, reaching the centre in roughly |$0.5\, {t_{\rm sc}}$| (see Appendix  A). However, in streams and spheres, the Mach number increases as the shock radius decreases due to geometrical effects (Guderley 1942; Modelevsky & Sari 2021). The implosion shocks thus reach the centre faster in these geometries, as can be seen by the shock positions in panels (a2) and (c2).

This also causes the density and velocity in the post-shock region, between the shock radius and the CD, to increase from sheets to streams to spheres due to geometrical effects. At the CD, hot gas mixes with cold gas through a turbulent mixing layer, causing an entrainment flow to develop in the hot gas outside the cloud. This is why the hot gas inflow velocity is roughly twice as large as that of the cold gas. The mass flux of the entrainment flow, |${\dot{M}} \propto \rho _{\rm h} v_{\rm h} r^{n}$|⁠, is constant, with |$n=0$|⁠, 1, and 2 for sheets, streams, and spheres, respectively. Since the density in the hot gas is roughly constant with radius (panel a2), this implies that the hot gas velocity scales as |$v_{\rm h}\propto r^{-n}$|⁠, which is broadly consistent with the hot gas velocities seen in panel (d2), which increase in magnitude from spheres to streams to sheets. Finally, we note that the density, temperature, and pressure in the hot gas at |$r\gt r_{\rm s,i}$| all decrease during the implosion phase, due to the outward propagating rarefaction wave.

At |$t\lower.5ex\rm{\,\, \buildrel\lt \over \sim \,\,}0.5\, t_\mathrm{sc}$|⁠, the implosion shock reaches the cloud centre. As this occurs, the central density and pressure reach very large values, though the central temperature remains at |$T_{\rm c}$|⁠. We note that the boost in central pressure and density is much larger in streams compared to sheets, and much larger still in spheres (panels a3 and c3). Formally, one can show that for self-similar collapse the central pressure in spheres and streams diverges (Guderley 1942), while it reaches a finite value in sheets (Toro 2013, see also Appendix  A). In practice, however, the maximal pressure has a finite value due to the finite thickness of the shock, which in our simulations is further limited by grid resolution. The peak pressure occurs once the implosion shock reaches the centre, while the explosion phase begins after this, once the reflected shock reaches the cold–hot interface. In between these two times, the inner regions are expanding while the outer regions are still contracting. We are thus never in a situation where all the initial cold gas is condensed into a single cell, explaining why the peak pressure is smaller than |$\sim (r_{\rm s,i}/\Delta)^2P_0=32^2P_0$| for streams and |$32^3P_0$| for spheres.

The rightmost column shows the situation once pressure equilibrium between the hot and cold phases has been re-established, and the outward velocity of the post-shock cold-gas, |$v_{\rm ex}$|⁠, has reach its peak value. In both spheres and streams, we find |$v_{\rm ex}\sim c_{\rm s,c}$| as expected from equation (10) (panel d4). However, in sheets we find |$v_{\rm ex}\sim 0.4\, c_{\rm s,c}$|⁠, due to the relatively small peak central pressure in sheets (⁠|$\sim 3$| times larger than the equilibrium pressure) compared to streams and spheres (⁠|$\sim 30-200$| times larger than the equilibrium pressure). This implies that only a small fraction of the peak internal energy in sheets goes into kinetic energy before the system regains pressure balance, leading to an explosion velocity smaller than |$c_{\rm s,c}$|⁠. At the same time, the central temperature in spheres and streams is |$\sim 0.3\, T_{\rm c}$| due to the adiabatic expansion phase described above, while it is |$\sim T_{\rm c}$| for sheets implying no adiabatic expansion of the cold phase (panel b4). As discussed in Appendix A2.3, this seems to be due to some fragmentation occurring in the sheet already during the implosion, thus decreasing the central overpressure and strength of the collision shock. For cylinders and spheres, even if such fragmentation occurs during the implosion, geometrical focusing will still enhance the collision shock and the corresponding central overpressure. We note that this result does not appear to be an artefact of limited resolution, as we obtain the same result for sheet simulations with cell sizes 2 and 4 times smaller (see Appendix B1).

3.1.2 The number of clumps

As the reflected shock sweeps over density inhomogeneities at the interface of the two phases created during the contraction, the local density and pressure gradients become misaligned leading to RMI which, in conjunction with radiative cooling, drives the fragmentation of cold gas (see additional discussion in Section 6.1 and Appendix A2.3). In the weak shock limit, RMI can be modelled as a form of Rayleigh–Taylor instability in which the gravitational force is impulsive, i.e. |$g\sim \Delta v\delta (t-t_{\rm shock})$|⁠, where |$\delta (t)$| is the Kronicker delta function, |$t_{\rm shock}$| is the time when the shock overtakes the interface, and |$\Delta v$| denotes the interface velocity jump (Zhou 2017a). We assume |$\Delta v \sim v_\mathrm{ex} \sim c_\mathrm{s,c}$|⁠, and further assume that due to density inhomogeneities and shape perturbations throughout the cold gas and across the interface, the effective gravitational acceleration is better modelled as |$g\sim \Delta v/t_{\rm cross}$|⁠, where |$t_{\rm cross} \sim \ell /\Delta v$| with |$\ell \sim r_{\rm s,f}$|⁠, the characteristic cloud size after cooling and contraction, and likely the dominant perturbation wavelength. The growth time-scale of the RMI is thus |$t_\mathrm{RM}\sim (\ell /g)^{1/2}\sim r_{\rm s,f}/c_\mathrm{s,c}$|⁠, comparable to the sound crossing time of the collapsed cloud. However, the constant of proportionality can deviate from unity and depends on the initial cloud geometry, as discussed in Section 6.1. The fragmentation time-scale, over which the number of clumps/cloudlets rapidly increases, is proportional to |$t_{\rm RM}$|⁠. However, we stress that fragmentation into discrete cloudlets occurs is not expected during the linear growth of RMI, or even in non-linear RMI for adiabatic gas. Rather, this seems to be a consequence of non-linear RMI together with strong cooling and the generation of vorticity due to the initial shape perturbations (see Appendix A2.3). A deeper analysis of the detailed microphysics behind the initial formation of the clumps is beyond the scope of the current work, where we mostly focus on the evolution of these cloudlets after the initial fragmentation.

The left-hand panel of Fig. 4 shows the number of clumps as a function of time, normalized by the initial cold cloud sound crossing time, |$t_{\rm sc}\equiv r_{\rm s,i}/c_{\rm s,c}$|⁠. Different colour lines represent different geometries, while the line thickness grows with increasing final overdensity, |${\chi _{\rm f}}$|⁠. In streams and spheres, the number of clumps peaks at |$N_\mathrm{c}\sim 100$| immediately after the simulation starts, due to a combination of RMI and thermal instabilities during the implosion. The number of clumps then decreases due to coagulation enhanced by further contraction, before rapidly rising again to values of a few |$10^3-10^4$| during the explosion. In sheets, the coagulation during the implosion is much weaker due to the lack of geometrical focusing, so the number of clumps monotonically increases until the end of the explosion phase. In all cases, |$N_\mathrm{c}$| stops growing rapidly by |$t\sim 4\, {t_{\rm sc}}$|⁠, when fragmentation stops and/or coagulation begins.

While the peak number of clumps increases from spheres to streams to sheets, this is proportional to the total amount of cold material. To factor this out, we present in the right-hand panel of Fig. 4 the evolution of |$N_{\rm c}/N_{\rm 0}$|⁠, where |$N_{\rm 0}=m_{\rm s,i}/m_{\rm cell}$| with |$m_{\rm s,i}$| the initial cold-gas mass in the analysis region (more than |$10\, {r_{\rm s,i}}$| from any boundary, Section 2.4) and |$m_{\rm cell}=\rho _{\rm s,f}\Delta ^{3}$| the mass of a cell at the equilibrium density of cold gas at |$T_{\rm floor}$|⁠. |$N_0$| thus represents the maximal number of cold-gas clumps possible if the cold-gas mass does not increase due to entrainment. With this normalization, we see that the efficiency of fragmentation increases from sheets to streams to spheres. The rate of fragmentation and clump formation is similar in all three geometries, while the time-scale for |$N_{\rm c}$| to reach its peak and saturate increases from sheets to streams to spheres. This will be further discussed Section 3.2 in the context of coagulation.

The number of clumps increases with |$\chi _{\rm f}$| for all geometries, though it tends to converge at |$\chi _{\rm f}\lower.5ex\rm{\,\, \buildrel\gt \over \sim \,\,}600$| for streams and spheres.⁴ In spherical geometry, the case of |${\chi _{\rm f}}=100$| recoagulates into a single large cloud after roughly |$10\, {t_{\rm sc}}$|⁠, while cases with |${\chi _{\rm f}}\ge 200$| remain fragmented with |$N_{\rm c}$| either continuing to grow or saturating at |$N_{\rm c}\sim 3000$|⁠. This is qualitatively similar to the results found in Gronke & Oh (2020b), only with a lower threshold in |${\chi _{\rm f}}$| for sustained fragmentation. This will be discussed further in Section 6.1. Streams exhibit qualitatively similar behaviour, with |$N_{\rm c}$| decreasing to order unity for |${\chi _{\rm f}}=100$| and remaining at large values for |${\chi _{\rm f}}\ge 200$|⁠. However, unlike the spherical case, streams with |${\chi _{\rm f}}\ge 200$| do exhibit some coagulation, with |$N_{\rm c}$| decreasing after an initial peak. This is particularly noticeable for |${\chi _{\rm f}}=200$|⁠, which we consider a borderline case (Table 1). Furthermore, unlike the spherical case which coagulated into a single cloud for |${\chi _{\rm f}}=100$|⁠, streams with |${\chi _{\rm f}}=100$| maintain several distinct clumps along the stream axis. The coagulation along the stream axis is suppressed compared to the radial direction due to opposing forces pulling clumps in either direction.⁵ At late times, the number of clumps fluctuates between |$N_{\rm c}\sim (1-4)$| due to centres of large clumps along the stream axis moving in and out of our analysis region, |$|z| \le 6\, {r_{\rm s,i}}$| (see Section 2.4).

Sheets exhibit even stronger coagulation than streams for |${\chi _{\rm f}}\ge 200$|⁠, with |$N_{\rm c}$| decreasing by more than |$\sim 50~{{\ \rm per\ cent}}$| between its peak and |$10\, {t_{\rm sc}}$| even for |${\chi _{\rm f}}=1000$|⁠. However, similar to streams, the coagulation is primarily in the radial direction and is suppressed in the plane of the sheet. As a result, several tens of clumps remain at the end of the |${\chi _{\rm f}}=100$| simulation.

We show density maps at |$t=10\, t_{\rm sc}$| in two orthogonal projections and for different values of |${\chi _{\rm f}}$|⁠, for sheets, streams, and spheres in Figs 5, 6, and 7, respectively. For sheets and streams, the projections correspond to face-on and edge-on, while for the spheres we simply show two orthogonal orientations. These maps show the maximal density along the line of sight, which highlights the small clumps resulting from fragmentation. Complementary to these, we show in Figs C1–C3 the average density along the line of sight for the same projections, which better highlights coagulation and the geometry of large clouds. Qualitatively, one sees radial coagulation grow stronger from spheres to streams to sheets, with the edge-on projection revealing strong coagulation in sheets even when |${\chi _{\rm f}}=1000$| (bottom-right panel of Figs 5 and C1). While some radial coagulation is always apparent in each geometry, we find more small clumps at larger radial distances as |${\chi _{\rm f}}$| increases. On the other hand, coagulation is suppressed along the stream axis and within the plane of the sheet even for |${\chi _{\rm f}}=100$|⁠, where |$N_{\rm c}=1$| for spheres, a few for streams, and a few tens for sheets. While this is hard to see in Figs 5 and 6, it becomes clear when examining Figs C1and C2.

3.1.3 The distance of clump propagation

One of the key features of thermal fragmentation is that cold-gas clumps get spread over an area which grows larger with time as the cloudlets spread out and fragmentation continues, whereas if coagulation is important the region occupied by cold gas reaches a maximum and then begins to shrink. In Fig. 8, we show the time evolution of |$d_{\rm max}$|⁠, the maximal radial distance of any clump whose mass is at least |$8\, {\rm m_{cell}}\equiv 8\, \rho _{\rm c}\Delta ^3$|⁠, where |$\Delta$| is the minimal cell size (valid in the region |$|x_1|\lt 3\, {r_{\rm s,i}}$|⁠, Section 2.4) and |$\rho _{\rm c}$| is the equilibrium density in the cold phase. We implement this threshold to reduce our sensitivity to resolution effects by removing small clumps at the grid scale. Using a threshold of 16 cells gives very similar results, while taking all clumps gives qualitatively similar results. |$d_{\rm max}$| is measured as the radial distance from the centre of simulation domain for spheres, the radial distance from the z- axis for streams, and the vertical distance from the |$yz$|-plane for sheets, and is normalized by the initial cloud radius, |$r_{\rm s,i}$|⁠.

Initially, |$d_{\rm max}=0$| because there is only one ‘clump’ whose centre is at the centre of the simulation volume. However, immediately following that we get |$d_{\rm max}\sim r_{\rm s,i}$| for all cases, as fragmentation begins near the initially perturbed cloud interface. |$d_{\rm max}$| then proceeds to shrink during the implosion before rapidly rising during the explosion phase. The growth rates during the explosion are similar for all cases, peaking at |$({\rm d}/{\rm d}t)d_{\rm max}\sim (0.5-2.0)c_{\rm s,c}\sim v_{\rm ex}$| (Section 3.1.1).

For sheets with |$\chi _{\rm f}= 100$|⁠, |$d_{\rm max}$| never exceeds |$r_{\rm s,i}$| due to very strong coagulation. In all other cases, |$d_{\rm max}$| peaks at values several times the initial cloud radius. For sheets with |$\chi _{\rm f}\ge 200$|⁠, |$d_{\rm max}$| peaks at |$\sim (2-3)\, {r_{\rm s,i}}$| at |$t\sim 4\, {t_{\rm sc}}$| and then noticeably decreases, indicative of strong coagulation and consistent with the decline in |$N_{\rm c}$| seen in Fig. 4. The same is true for streams and spheres with |$\chi _{\rm f}=100$|⁠. In all coagulating cases, |$d_{\rm max}$| appears to oscillate at late times, consistent with the pulsations observed at relatively low values of |$\chi _{\rm f}$| in previous work (Gronke & Oh 2020b, 2023). On the other hand, streams and spheres with |${\chi _{\rm f}}\ge 200$| exhibit |$d_{\rm max}$| values that either continue to rise or saturate until the end of the simulation. This suggests a critical |$\chi _{\rm f}\sim 200$| for sustained fragmentation in streams and spheres, consistent with Fig. 4. In the stream simulation with |${\chi _{\rm f}}=200$|⁠, |$d_{\rm max}$| slightly decreases at |$t\gt 10\, {t_{\rm sc}}$| consistent with this being more of a borderline case as noted above. We note that the strong saturation observed at |$d_{\rm max}\sim 10\, {r_{\rm s,i}}$| is likely numerical, because the cell size at |$d\gt 9\, {r_{\rm s,i}}$| grows to |$8\Delta$|⁠, causing even large clumps to artificially disrupt. However, this does not change the qualitative distinction between cases where |$d_{\rm max}$| decreases due to strong coagulation and cases where it does not.

3.2 Cold-gas coagulation

In the previous section, we saw that sheets were prone to strong coagulation at all overdensities, while spheres and streams were prone to coagulation for |$\chi _{\rm f}\lower.5ex\rm{\,\, \buildrel\lt \over \sim \,\,}200$|⁠. This was found to be the case based on the number of clumps (Fig. 4), their morphology (Figs 5–7) and their radial extent (Fig. 8). In this section, we study the coagulation process in detail.

3.2.1 Theoretical overview

We first review the theoretical model described in Gronke & Oh (2023), based on their analysis of quasi-spherical systems of clouds. Consider a cold clump, hereafter clump 1, embedded in a hot gaseous medium. Any perturbations at the interface of the cold and hot phases will induce a turbulent mixing layer, where efficient radiative cooling will drive an entrainment flow of hot background gas onto the clump with velocity

for radii |$r\ge r_\mathrm{cl,1}$|⁠. Here,

is the entrainment velocity at the surface of the clump (Gronke & Oh 2018, 2020a, 2023; Fielding et al. 2020; Sparre et al. 2020; Mandelker et al. 2020a; Tan, Oh & Gronke 2021; see also Appendix A2.2 and Fig. A4).

Now consider a static clump, hereafter clump 2, at a distance d from the central clump 1. The force acting on clump 2 is

where |$\Delta v_2$| is the velocity of clump 2 relative to the hot background. The first term on the right-hand side represents an effective force due to condensation, which causes the mass of the cold cloud to increase and thus its velocity to decrease due to momentum conservation,

The second term on the right-hand-side of equation (13) is the drag force at fixed cloud mass which is given by

where |$C_\mathrm{d}\approx 0.47$| is the drag coefficient for a sphere. The ratio of these two forces is thus

Initially, clump 2 is at rest, while the hot gas moves at velocity |$v_\mathrm{hot}(r=d)$|⁠. Over time, clump 2 becomes entrained in the flow and the relative velocity decreases. We thus have in general |$\Delta v_2\le v_\mathrm{hot}(r=d)$|⁠, yielding

where in the final equation we have used equation (12) and the fact that the two clouds have similar temperatures and densities, though not necessarily similar sizes. We conclude that |$F_\mathrm{cond}\gg F_\mathrm{drag}$| unless |$r_\mathrm{cl,1}/r_\mathrm{cl,2}\gt 10^{4}$| and d is not much larger than |$r_{\rm cl,1}$|⁠. This is true regardless of the details of the cooling function, and in particular regardless of the gas metallicity, so long as the two clumps are initially static with respect to the background.

So long as the velocity of clump 2 with respect to clump 1 is negligible compared to the entrainment velocity of hot gas onto clump 1, we have |$\Delta v_2\sim v_\mathrm{hot}$|⁠. The condensation force then becomes

where |$A_{\rm cl}=4\pi r_{\rm cl}^2$| is the cloud surface area. Neglecting the weak dependence on the ratio |$r_{\rm cl,2}/r_{\rm cl,1}$|⁠, this force is similar to gravity if we make the analogy that |$G\rightarrow \rho _\mathrm{h}v_\mathrm{mix,1}^{2}/4\pi$| and |$m\rightarrow A$|⁠, with both forces proportional to |$d^{-2}$|⁠.

If the central region is occupied by a collection of similar sized clumps each with |$r_{\rm cl}\sim r_{\rm cl,1}$|⁠, the entrainment velocity of the background flow towards the centre becomes

where |$A_{\rm cl,tot}(\lt r)$| is the total combined surface area of all cold clouds interior to r, and |$f_\mathrm{A}(r)\equiv A_\mathrm{cl,tot}(\lt r)/4\pi r^{2}$| is the area modulation factor. This factor can be greater or less than unity, as discussed in Section 3.2.2. Consequently, the total force acting on a test clump 2 at a distance d from the centre of clumps is

where we have approximated |$v_\mathrm{mix}\sim v_\mathrm{mix,1}\sim v_\mathrm{mix,2}$|⁠. The acceleration of clump 2 is thus

We can use this to derive a characteristic time-scale for coagulation, analogous to the gravitational free-fall time. Assuming that the acceleration is roughly constant, namely that |$f_{\rm A}(d)\sim \rm const.$|⁠, we have⁶

We can similarly derive the coagulation force and time-scale for stream and sheet geometries. It turns out that the only difference⁷ is the factor |$f_\mathrm{A}$|⁠, which can be defined as

|$\Lambda _\mathrm{tot}$| is the total area of all cold clouds per unit length, proportional to the cloud diameter, and |$\Sigma _\mathrm{tot}$| is the total area of cold clouds per unit area, proportional to the number of clouds. Note the different scaling with d in different geometries, similar to the gravitational force from a spherical, cylindrical, and planar distribution of mass.

The above considerations, based on Gronke & Oh (2023), are valid for a ‘test clump’ (clump 2) initially at rest with respect to the central collection of cold clouds (clump 1). Now let us consider the case where clump 2 is escaping the central region, while at the same time, there is an entrainment flow of hot gas towards the centre. This is the situation immediately after the explosion phase described in the previous section, where fragmented clumps were escaping the central region with a velocity |$v_{\rm ex}\sim c_{\rm s,c}$|⁠. The total time-scale for the escaping clumps to recoagulate is the sum of the deceleration time-scale, when the interaction of the clumps and the entrainment flow causes the clumps to decelerate and turn around, and the coagulation time-scale from equation (22) which measures the coagulation time of static clumps in the hot entrainment flow. The turnaround/deceleration process can be due either to the condensation force or the drag force, depending the ratio in equation (16). By combining equations (12) and (16), we have

Considering |$t_\mathrm{sc,2}/t_\mathrm{cool}\approx r_\mathrm{cl,2}/\ell _\mathrm{shatter}\gtrsim 1$|⁠, and |$\Delta v_2\sim c_\mathrm{s,c}$| (see Fig. 2), we conclude that regardless of geometry and metallicity |$F_\mathrm{cond}$| dominates over |$F_\mathrm{drag}$| in the deceleration process, as well as in the coagulation process discussed above.

The deceleration time-scale is thus determined by the condensation time-scale, which is the time-scale for cloud growth by condensation

The ratio of the deceleration and coagulation time-scales is thus

Since |$t_\mathrm{decel}$| has no explicit dependence on d while |$t_\mathrm{coag}$| increases with d, we obtain that |$t_\mathrm{decel}\gt t_\mathrm{coag}$| at small d where |$f_{\rm A}$| is large (equation 23). However, if a clump finds itself at large d with small |$f_{\rm A}$|⁠, then |$t_\mathrm{coag}$| will dominate the remaining recoagulation time-scale even if the clump is still in the process of decelerating.

Before significantly decelerating, clumps can reach a maximal distance of |$d_\mathrm{max}\sim \xi c_\mathrm{s,c}t_\mathrm{decel}$|⁠, where |$v_{\rm ex}\sim \xi c_{\rm s,c}$| is the explosion velocity and |$\xi \sim 0.4$| or 1.0 for sheets or streams and spheres, respectively (Section 3.1.1). The maximum value of |$t_\mathrm{coag}$| along the clump’s orbit is reached at |$d=d_\mathrm{max}$|⁠,

Reaching |$d_{\rm max}$| and turning around is a necessary condition for the clump to coagulate. The question of whether the final coagulation is efficient is determined by the ratio in equation (26),

Neglecting the second term with the |$1/8$|-power, we see that larger values of |$f_\mathrm{A}(d_\mathrm{max})$| result in more efficient coagulation. We now turn to quantify |$f_\mathrm{A}$| as a function of d in different geometries, to gain a better understanding of the geometrical effects on coagulation.

3.2.2 Cold-gas area

We can crudely estimate the efficiency of coagulation by measuring the coagulation forces on the outermost clump. To this end, we estimate |$f_{\rm A,dmax}\equiv f_{\rm A}(d_{\rm max})$|⁠, with |$d_{\rm max}$| as in Fig. 8. At early times, when |$d_{\rm max}$| is still small and the fragmentation process is still ongoing, we expect |$A_{\rm tot}$| to increase with d and therefore |$f_{\rm A}$| to either increase or decrease depending on the details. However, at later times, once |$d_{\rm max}$| grows beyond the central concentration of cold gas, both the number of clumps (Fig. 4) and their radial extent (Fig. 8) are roughly constant, consistent with the fragmentation process ending and/or coagulation affecting the inner region. At this stage, we expect |$A_{\rm tot}$| to be roughly constant and |$f_{\rm A}\propto d_{\rm max}^{-n}$|⁠, with |$n=0$|⁠, 1, or 2 for sheets, streams, and spheres (equation 23).

In Fig. 9, we show the total surface area of cold gas, |$A_{\rm tot}(\lt d_{\rm max})$|⁠, on the left, and the area modulation factor, |$f_{\rm A}$|⁠, on the right. Here, |$A_{\rm tot}$| refers to the total surface area of cold gas, regardless of the size of the clumps, since this is the relevant quantity governing coagulation. To measure this in the simulations, we first interpolate the gas density onto a uniform grid at the highest resolution of our refined grid using yT (Turk et al. 2011), and then extract a 2D surface mesh from a 3D volume using the python package scikit-image (van der Walt et al. 2014). We extract the density isosurface corresponding to |$\rho _{\rm iso}=\rho _{\rm s,i}/1.05$|⁠, where the factor 1.05 is to ensure that this captures all of the initial cloud including any density fluctuations present in the initial conditions (Section 2.5). We then normalize |$A_{\rm tot}$| by its initial value, |$A_0$|⁠, which lacks an analytical form due to the shape perturbations on the initial cloud interface. As noted in Fig. 8, |$d_{\rm max}=0$| at the initial condition, and is only self-consistently defined once fragmentation begins near the cloud interface. Therefore, |$A_0=A_{\rm tot}(t=0)$| simply represents the total surface area of the initial cloud, and |$A_{\rm tot}$| is only limited by |$d_{\rm max}$| from the second snapshot.

During the implosion phase, |$A_{\rm tot}$| decreases rapidly for spheres and streams due to contraction, which reduces the surface area of the cloud. This effect is stronger for spheres than for streams because of the different scaling of cloud area with radius. However, for sheets |$A_{\rm tot}$| actually increases during the implosion phase, because the overall area of the central sheet is independent of its ‘radius’ (thickness), while fragmentation increases the total cold-gas surface area.

Following the explosion phase, |$A_{\rm tot}$| reaches values of |$\sim (2-4)A_0$| at |$t\sim (2-4)\, {t_{\rm sc}}$|⁠. In cases with strong coagulation (based on Figs 4 and 8), namely spheres with |${\chi _{\rm f}}=100$|⁠, streams with |${\chi _{\rm f}}\lower.5ex\rm{\,\, \buildrel\lt \over \sim \,\,}200$|⁠, and all sheet simulations, |$A_{\rm tot}$| proceeds to decrease towards values of |$\lower.5ex\rm{\,\, \buildrel\lt \over \sim \,\,}A_0$| at late times. On the other hand, in cases which show strong fragmentation, namely spheres with |${\chi _{\rm f}}\ge 200$| and streams with |${\chi _{\rm f}}\ge 400$|⁠, |$A_{\rm tot}$| saturates following the explosion phase at values of |$A_{\rm tot} \sim (2-3)A_0$|⁠. We note that the borderline nature of the stream simulation with |${\chi _{\rm f}}=200$| is particularly evident here. While a boost in the surface area by a factor of |$\sim 3$| may not seem like much, we recall that without fragmentation the final equilibrium configuration has a radius |$\, {r_{\rm s,f}}$|⁠, which is |$\sim 10$|⁠, 3, and 2 times smaller than |$\, {r_{\rm s,i}}$| for sheets, streams, and spheres, respectively (Table 1), yielding a final equilibrium surface area which is |$\sim 1$|⁠, 3, and 5 times smaller than |$A_0$|⁠. The actual increase in cold-gas surface area with respect to the equilibrium configuration due to thermal fragmentation is thus a factor of |$\sim 6$| for streams and 10 for spheres.

The precise late-time value of |$A_{\rm tot}$| is somewhat sensitive to the late-time value of |$d_{\rm max}$| (Fig. 8) and therefore to our threshold of |$m_{\rm cl}\ge 8m_{\rm cell}$|⁠. Likewise, the decrease in the late-time value of |$A_{\rm tot}$| with increasing |${\chi _{\rm f}}$| for spheres seems to also be an artefact of our refinement strategy which causes even large clumps to artificially disrupt at distances of |$d\gt 9\, {r_{\rm s,i}}$| thus decreasing the cold-gas surface area. Nevertheless, there is a real qualitative distinction between |$A_{\rm tot}$| saturating at |$\sim (2-3)A_0$| for fragmenting cases, versus decreasing to |$\lower.5ex\rm{\,\, \buildrel\lt \over \sim \,\,}A_0$| for coagulating cases.

On the right, we show |$f_{\rm A,dmax}$| by combining |$A_{\rm tot}(\lt d_{\rm max})$| shown on the left with equation (23). We show this as a function of |$d_{\rm max}/r_{\rm s,i}$| itself, so each point on this plot represents the situation of the farthest clump at the corresponding time, and can be inserted into equation (22) or (26). Note that initially, when |$d_{\rm max}=\, {r_{\rm s,i}}$| prior to the implosion, |$f_{\rm A,dmax}=1$| for spheres and streams without shape perturbations, though |$f_{\rm A,dmax}=2$| for sheets with no perturbations, because of the definition of |$\Sigma$| in equation (23) and the fact that a sheet has two surfaces. The initial perturbations increase this value, such that spheres and streams begin at |$f_{\rm A,dmax}\sim 1.6$| while sheets begin at |$f_{\rm A,dmax}\sim 3.2$|⁠, which is consistent with each surface adding |$\sim 0.6$| to |$f_{\rm A,dmax}$|⁠. During the implosion phase |$d_{\rm max}$| decreases, while |$A_{\rm tot}$| increases for sheets and decreases for streams and spheres. This leads to fairly chaotic behaviour in the the plane of |$d_{\rm max}/r_{\rm s,i}$| and |$f_{\rm A,dmax}$| (grey shaded region). However, the value of |$f_{\rm A}$| during the implosion phase is uninteresting, since the distinction between coagulation and fragmentation is only meaningful after the explosion phase, once |$d_{\rm max}\gt \, {r_{\rm s,i}}$|⁠.

For sheets, as well as for streams with |${\chi _{\rm f}}\le 200$| and spheres with |${\chi _{\rm f}}=100$|⁠, coagulation is important and |$d_{\rm max}$| decreases back towards |$\, {r_{\rm s,i}}$| after peaking at somewhat larger values (Fig. 8). This is also seen in the right-hand panel of Fig. 9, where the curves turn around towards lower |$d_{\rm max}$| after first reaching values |$d_{\rm max}/\, {r_{\rm s,i}}\gt 1$|⁠. For these cases, while the value of |$f_{\rm A}$| is smaller during the coagulation then during the initial expansion, there is no clear trend of |$f_{\rm A}$| with d during either phase. On the other hand, for cases which fragment, namely spheres with |${\chi _{\rm f}}\ge 200$| and streams with |${\chi _{\rm f}}\ge 400$|⁠, |$d_{\rm max}$| never decreases (Fig. 8), while |$A_{\rm tot}$| is roughly constant at late times. This results in |$f_{\rm A}\propto d^{-1}$| for streams and |$f_{\rm A}\propto d^{-2}$| for spheres (equation 23), as highlighted in the figure. Overall, |$f_\mathrm{A,dmax}$| at late times is largest in sheets and smallest in spheres, as is the coagulation efficiency (equation 28), as demonstrated in previous sections. We will discuss this further in Section 6.1.

4 METALLICITY AND UVB EFFECTS ON THERMAL FRAGMENTATION AND COAGULATION

Our analysis in the previous section focused on solar metallicity gas in collisional ionization equilibrium, similar to previous studies of thermal fragmentation versus coagulation in quasi-spherical clouds (Gronke & Oh 2020b, 2023). While solar metallicity may be reasonable for gas in the inner CGM at |$z\sim 0$|⁠, the metallicity in the high-z cosmic web which is our primary focus is much lower. This lowers the cooling rates for intermediate temperature gas in the turbulent mixing layers between the cold clouds and hot background, thus lowering the efficiency of phase mixing and entrainment (Fig. 10). Moreover, intergalactic gas is affected by the ionizing UVB, and photoionization is more important than collisional ionization over a wide temperature range (e.g. Strawn et al. 2021; Strawn, Roca-Fàbrega & Primack 2023). This lowers the cooling rates near the cold phase, |$T\lower.5ex\rm{\,\, \buildrel\gt \over \sim \,\,}10^4\, {\rm K}$| (Fig. 10). These changes to the cooling curve affect both the shattering length-scale, |$\ell _{\rm shatter}$|⁠, and the strength of coagulation forces. Therefore, in this section, we revisit some of our previous analysis focusing on low-metallicity gas exposed to a UVB. We focus here specifically on stream geometry and in Section 7 below we apply these results to cosmic-web filaments at high-z.

As stated in Section 2, we assume metallicity values of |$Z_{\rm bg}=0.1Z_{\odot }$| for the background and |$Z_{\rm s}=0.03Z_{\odot }$| for the initial stream, and apply the ionizing UVB of Haardt & Madau (1996) at |$z=2$|⁠. We do not include self-shielding of dense gas. At our fiducial densities and metallicities, the equilibrium temperature of the cold stream is |$1.68\times 10^4 \, {\rm K}$|⁠, which was our chosen |$T_{\rm floor}$| in the previous section. In this section, we do not use an artificial temperature floor since the UVB sets an effective cooling floor. Hereafter, when we refer to simulations of low-metallicity (low-Z) streams, it is understood that these are also exposed to a UVB, while high-metallicity (high-Z) streams are not.

While our initial conditions are the same as in Section 3 in terms of stream size, density, and temperature, the initial cooling rates are much lower (Fig. 10). Therefore, unlike our simulations in Section 3 where the streams immediately cool isochorically, in our current simulations, the streams initially contract isobarically and only loose sonic contact when they reach a radius of |$r^{*}\sim 2.5\, \mathrm{kpc}$| (see Table 1). After this point, the implosion and explosion processes occur similarly at low- and high-Z. We note that despite the weaker radiative cooling, the explosion velocities are still roughly |$v_{\rm ex}\sim c_\mathrm{s,c}$|⁠, and the dominant deceleration force as clumps escape outward is still |$F_\mathrm{cond}$| (equation 24). However, the corresponding |$t_{\rm decel}$| at a fixed |${\chi _{\rm f}}$| is longer for low-Z streams, due to |$v_\mathrm{mix}$| being smaller (equation 25).

Fig. 11 shows maps of the maximal density along the line of sight in two orthogonal projections for low-Z streams with different values of |${\chi _{\rm f}}$|⁠, similar to Fig. 6. Complementary to this, we show in Fig. C4, the average density along the line of sight in the same projections. These simulations exhibit sustained fragmentation at |${\chi _{\rm f}}\lower.5ex\rm{\,\, \buildrel\gt \over \sim \,\,}60$|⁠, in contrast to |${\chi _{\rm f}}\lower.5ex\rm{\,\, \buildrel\gt \over \sim \,\,}200$| at solar metallicity with no UVB, with |${\chi _{\rm f}}=60$| at low-Z being similarly borderline to |${\chi _{\rm f}}=200$| at high-Z. This is further demonstrated in Fig. 12, which shows the number of clumps, |$N_{\rm c}$|⁠, on the left and the maximal clump distance, |$d_{\rm max}$|⁠, on the right for different |${\chi _{\rm f}}$| values. Solid lines show results for low-Z streams, while dashed lines show the same high-Z simulations as Figs 4 and 8. Note that the low-Z runs extend to lower |${\chi _{\rm f}}$| values. In general, the explosion phase and subsequent rapid rise of both |$N_{\rm c}$| and |$d_{\rm max}$| takes |$\sim 2\, t_\mathrm{sc}$| longer for the low-Z streams because of the longer cooling times and the initial phase of isobaric cooling. We find a lower |$\chi _\mathrm{f,crit}$| for fragmentation in these simulations, with the |${\chi _{\rm f}}=60$| run demonstrating weak to no coagulation, and the |${\chi _{\rm f}}=100$| run displaying strong fragmentation. We will discuss this further in Section 6.1.

The shattering length-scale in the low-Z runs is |$\sim 80$| times larger than in the high-Z runs, as is the cooling time near the cold phase. From equation (25), we thus expect low-Z clumps to propagate to distances |$\sim 3$| times greater than high-Z clumps. This is indeed the case for |${\chi _{\rm f}}=100$|⁠, as can be seen by comparing the solid and dashed red lines in the right-hand panel of Fig. 12. However, this is not evident in streams with higher |${\chi _{\rm f}}$| where |$d_{\rm max}\gt 10\, {r_{\rm s,i}}$|⁠, due to artificial clump disruption in low-resolution regions, as discussed in the context of Fig. 8. We also note the similarity between the evolution of |${\chi _{\rm f}}=30$| (60) at low-Z and |${\chi _{\rm f}}=100$| (200) at high-Z. These simulations have a similar ratio of |$\chi _\mathrm{f}/v_\mathrm{mix}$| (equation 12) and a similar distribution of clump sizes (Section 5) yielding a similar deceleration time-scale (equation 25).

Fig. 13 shows the total surface area of cold gas, |$A_{\rm tot}(\lt d_{\rm max})$| on the left, and the area modulation factor |$f_{\rm A,dmax}$| (equation 23) on the right, as in Fig. 9. We compare simulations with different |${\chi _{\rm f}}$| at both low (solid) and high (dashed) metallicity. Despite some differences during the initial isobaric contraction phase, the overall behaviour of the cold-gas surface area is similar at low- and high-Z. Following the explosion, the area increases rapidly, saturating at |$A_{\rm tot}\sim (2-3)A_0$| for fragmented streams, and at |$A_{\rm tot}\lower.5ex\rm{\,\, \buildrel\lt \over \sim \,\,}A_0$| for coagulated streams. As discussed following Fig. 9, this corresponds to a cold-gas surface area |$\sim 6$| times larger than the expected surface area in the final equilibrium state without fragmentation, namely a single stream of radius |$\, {r_{\rm s,f}}$|⁠.

Note that |$A_{\rm tot}$| does not grow monotonically with |$\chi _\mathrm{f}$| at low-Z, but rather reaches a maximum around |${\chi _{\rm f}}\sim (100-400)$|⁠. This is due to the disruption of cold-gas clumps with large |${\chi _{\rm f}}$| by hydrodynamic instabilities caused by the interaction with the surrounding hot wind (i.e. cloud crushing). Radiatively cooling clouds can survive these instabilities and grow in mass by entrainment if their overdensity obeys (Gronke & Oh 2018)

where |$P_{2.3}=P/(200{\rm cm^{-3}\, {\rm K}})$| is the thermal pressure of the hot gas, |$T_{\rm cl,4.3}=T_{\rm cl}/(2\times 10^4\, {\rm K})$| is the temperature of the cold cloud, |$\Lambda _{\rm mix,-22.4}=\Lambda _{\rm mix}/(10^{-22.4}{\rm erg\, s^{-1}\, cm^{3}})$| is the cooling rate of gas in the mixing layer with |$T_{\rm mix}\sim (T_{\rm c}T_{\rm h})^{1/2}$| and |$n_{\rm mix}\sim (n_{\rm c}n_{\rm h})^{1/2}$|⁠, |$M_{0.1}=v_{\rm cl}/(0.1c_{\rm s,h})$| is the cloud Mach number with respect to the hot gas, |$r_{\rm cl}$| is the cloud radius normalized to |$560\, {\rm pc}$| which is roughly six high-resolution cells, and |${\tilde{\alpha }}\sim 1$| in the ‘wind-tunnel’ setup. To demonstrate this, we show the evolution of the total cold-gas mass in Fig. 14. The cold mass increases with time for low-Z streams with |${\chi _{\rm f}}\lower.5ex\rm{\,\, \buildrel\lt \over \sim \,\,}1000$|⁠, while it decreases for larger |${\chi _{\rm f}}$| as clouds move into the disruption regime.⁸ For high-Z streams, the cold mass increases for all values of |${\chi _{\rm f}}$| due to the order of magnitude larger cooling rates.⁹ Both low- and high-Z streams exhibit |$M_{0.1}\sim (1-3)$| for all cases, with lower |${\chi _{\rm f}}$| corresponding to slightly lower values of |$M_{0.1}$| but no systematic difference between low- and high-Z. The increased cold-gas disruption in low-Z streams is thus only due to the lower cooling rates in these simulations.

The area modulation factor, |$f_\mathrm{A}$| (right-hand panel of Fig. 13), behaves similarly in the low- and high-Z runs, and is as expected based on Fig. 9. At early times when |$d_{\rm max}\lower.5ex\rm{\,\, \buildrel\lt \over \sim \,\,}\, {r_{\rm s,i}}$| (grey region), |$f_{\rm A}$| displays chaotic behaviour, stemming from a competition between a decrease in cold-gas surface area due to the contraction and an increase due to fragmentation. Sustained fragmentation and coagulation are easily distinguishable by the behaviour of |$f_{\rm A}$| at later times, once |$d_{\rm max}\gt \, {r_{\rm s,i}}$|⁠. Coagulation manifests as a turnaround or saturation in the distance at |$d_{\rm max}\sim (1-2)\, {r_{\rm s,i}}$| with no clear trend of |$f_{\rm A}$| with d. Sustained fragmentation, on the other hand, has |$f_{\rm A} \propto d^{-1}$| in the range |$d_{\rm max}\sim (3-10)\, {r_{\rm s,i}}$|⁠. Similar to |$A_{\rm tot}$|⁠, |$f_{\rm A}$| in low-Z streams decreases at overdensities |${\chi _{\rm f}}\gt 100$| due to cold-gas disruption. However, |$f_{\rm A}$| for |${\chi _{\rm f}}=100$| at low-Z is as large as it is for fragmented high-Z streams.

One of the key observational indicators of thermal fragmentation, and an important parameter in estimating the cold-gas content of the CGM/cosmic web, is the area covering fraction of cold gas, |$f_{\rm C}$|⁠. In fragmented gas, this can be very large even if the volume filling fraction of the cold gas is small (McCourt et al. 2018; Faucher-Giguère & Oh 2023). While |$f_{\rm C}$| overall behaves similarly to |$A_{\rm tot}$|⁠, two clumps aligned along the same line of sight would both contribute to |$A_{\rm tot}$| but not to |$f_{\rm C}$|⁠, which can cause large differences if there are many small clumps in the background/foreground of large clumps. Nevertheless, we crudely estimated |$f_{\rm C}$| perpendicular to the stream axis in a square region of side |$10\, {r_{\rm s,i}}$|⁠, and found similar enhancements as seen in Fig. 13 for |$A_{\rm tot}$|⁠. Namely, typical values at late times roughly |$\sim 2$| times larger than the covering fraction of the initial stream, and |$\sim 6$| times larger than the covering fraction of a stream with radius |$\, {r_{\rm s,f}}$| in thermal and pressure equilibrium. This was found to be true for both low- and high-Z fragmented streams. We defer a more detailed study of the cold-gas covering fraction and clumping factor, particularly in the context of cold streams feeding massive haloes at high-z (see Section 7), to future work with more realistic simulation setups and additional physical processes, as discussed in Section 8.

5 CLUMP SIZES AND MASSES

Our estimate of the coagulation time (equation 22) and the deceleration time (equation 25), as well as their ratio (equation 28, which determines the efficiency of coagulation in our model) all depend on the radii of clumps resulting from the initial fragmentation process. It therefore behooves us to quantify the distribution of clump sizes. In the initial formulation of the ‘shattering’ model, the final scale of cold-gas clouds was predicted to be |$r_{\rm c}\sim \ell _{\rm shatter}\equiv {\rm min}(c_{\rm s}t_{\rm cool})$| (McCourt et al. 2018). The emergence of this as a characteristic length-scale for cold gas was also discussed in previous analytic models of non-linear thermal instability (Burkert & Lin 2000; Waters & Proga 2019a). Motivated by these insights, some cosmological simulations have implemented cooling-length-based refinement models in order to resolve cold gas in the CGM (Peeples et al. 2019), while others have implemented subgrid models for the existence of unresolved cold-gas clouds of size |$\ell _{\rm shatter}$| (Butsky et al. 2024). Numerical simulations of non-linear thermal instability in 1D find that while the minimum cloud size was of order |$\ell _{\rm shatter}$|⁠, larger clouds were common due to merging of smaller clouds (Das et al. 2021). However, the size distribution of post-fragmentation cloudlets has not been constrained in 3D simulations of thermal instability which resolve |$\ell _{\rm shatter}$|⁠.

Several works have studied the size distribution of cold clumps in multiphase gas in other contexts. Gronke et al. (2022) studied the mass distribution of cloudlets forming not through an implosion–explosion process as in this work, but rather by placing large clouds that are only slightly out of thermal equilibrium (by a factor of |$\sim 2$|⁠) in a (driven) turbulent box. They found that when the clouds were in the growth regime (see equation 29), the clump mass function in their simulations converged to |$\mathrm{ d}N/\mathrm{ d}m\propto m^{-2}$| over a wide range of simulation parameters. This implies roughly constant mass per logarithmic mass bin, and was found to extend down to the resolution limit of the simulation. Similar results were found in larger scale simulations of AGN jets in cool-core clusters (Li & Bryan 2014) and in magnetohydrodynamic (MHD) simulations of a multiphase ISM (Fielding et al. 2023). However, these works did not resolve |$\ell _\mathrm{shatter}$| so they could not comment on whether this would serve as a minimal or characteristic size of cold-gas cloudlets (see also Jennings et al. 2023). Cloud-crushing simulations of large, thermally stable cold clouds interacting with a hot wind that marginally resolve |$\ell _\mathrm{shatter}$| find that the size distribution of cold clumps forming in turbulent mixing layers downstream does not appear converged, with clouds smaller than |$\ell _\mathrm{shatter}$| common (Sparre, Pfrommer & Vogelsberger 2019; Gronke & Oh 2020a). On the other hand, Liang & Remming (2020) performed 2D cloud crushing simulations where the initial cloud was thermally unstable with a fractal structure consisting of gas with temperatures ranging from |$(10^{4.3}-10^{6.5})\, {\rm K}$| and a median of |$\sim 10^5\, {\rm K}$|⁠. They found that the resulting cold-gas clouds had a characteristic column density of |$N_{\rm H}\sim 10^{17.5}\, {\rm cm}^{-2}$|⁠, consistent with predictions for a characteristic cloud size of |$\ell _{\rm shatter}$| (McCourt et al. 2018). However, |$\ell _{\rm shatter}$| was only marginally resolved in their simulations and the characteristic scale they discovered may have in fact been tied to the grid scale. Moreover, the cloud crushing setup is in general quite different from our own, and its relevance for the pure thermal instability picture discussed in McCourt et al. (2018) is unclear. To summarize, while it appears clear that in order for thermally unstable clouds to fragment they must be much larger than |$\ell _\mathrm{shatter}$| (McCourt et al. 2018; Sparre et al. 2019; Gronke & Oh 2020b; Farber & Gronke 2023), and that both the total mass and clumpiness of cold gas in multiphase media increase as the resolution approaches |$\ell _\mathrm{shatter}$| (Peeples et al. 2019; Mandelker et al. 2021), it remains unclear whether this is indeed a lower limit or a characteristic value of the cloud size distribution.

Our simulations offer a unique opportunity to study this issue, since unlike most previous works we focus on low-metallicity gas at relatively low pressures of |$P/k_{\rm B}\sim 100\, {\rm K}\, {\rm cm}^{-3}$| exposed to a strong |$z=2$| UVB (Section 4). These result in |$\ell _{\rm shatter}\sim 300\, {\rm pc}$| (Table 1), a factor of |$\sim (100-1000)$| larger than typical in most other works. For instance, in our simulations with solar metallicity presented in Section 3, |$\ell _{\rm shatter}\sim 4\, {\rm pc}$|⁠. In Fig. 15, we study the distribution of clump sizes in simulations of low-Z streams with |$\eta =10$| exposed to a |$z=2$| UVB, as in Section 4. All these simulations have |$\ell _{\rm shatter}\sim 300\, {\rm pc}$|⁠. We present results for streams with initial radii |$r_{\rm s,i}=1\, \mathrm{ and}\, 3\, {\rm kpc}$| (these are the simulations presented in Section 4), and |$30\, {\rm kpc}$|⁠. At our fiducial resolution of |$\Delta =r_{\rm s,i}/32$|⁠, these correspond to cell sizes of |$\sim 31,\,94,\,\mathrm{ and}\, 938 \, {\rm pc}$| in the high-resolution region, corresponding to |$\Delta /\ell _{\rm shatter}\sim 0.1$|⁠, 0.32, and 3.2, respectively. The simulations with |$r_{\rm s,i}=30\, {\rm kpc}$| thus do not quite resolve |$\ell _{\rm shatter}$|⁠, while our fiducial simulations with |$r_{\rm s,i}=3\, {\rm kpc}$| marginally resolve it. In the simulations with |$\, {r_{\rm s,i}}=1\, {\rm kpc}$|⁠, |$\ell _{\rm shatter}$| is well resolved, though these simulations have |$r^{*}\sim 2\ell _{\rm shatter}$| (Table 1) due to the initial phase of isobaric cooling discussed in Section 4, so fragmentation is expected to be weak (Gronke & Oh 2020b).

For each value of |$r_{\rm s,i}$|⁠, we perform additional higher resolution simulations with |$\Delta =r_{\rm s,i}/64$| and |$r_{\rm s,i}/128$|⁠, to test convergence. In these simulations, |$\ell _{\rm shatter}$| is well resolved for |$r_{\rm s,i}=3\, {\rm kpc}$| and marginally resolved even for |$r_{\rm s,i}=30\, {\rm kpc}$|⁠. We note that in these high-resolution simulations, we simulate a smaller portion of the stream with uniform resolution throughout the entire box. The stream radius here is only |$1/4$| of the box size as opposed to |$1/32$| in our fiducial simulations. The total number of clumps is thus not directly comparable between our fiducial resolution and the two higher resolution runs. However, we have verified that this change does not affect the physics of fragmentation versus coagulation nor the distribution of clump properties.

In the left-hand panel of Fig. 15, we show the cumulative clump mass distribution for these simulations, obtained directly from the clump finder. The main panel shows results from simulations with varying resolution, while the inset displays results from all simulations with |$\, {r_{\rm s,i}}/\Delta =32$| listed in the legend. We show the distributions at |$t=10\, {t_{\rm sc}}$| for fragmented streams, and at |$t\sim 5\, {t_{\rm sc}}$| for coagulating streams which is when |$N_{\rm c}$| is at its peak. In all cases, the number of clumps is clearly dominated by low-mass clumps, though a single massive clump contains anywhere from |$\sim (15\,\mathrm{ per\,cent}-95\,\mathrm{ per\,cent})$| of the total mass, with more strongly fragmented cases having a smaller maximal clump mass. Recall, however, that the total cold-gas mass can be a factor of |$\lower.5ex\rm{\,\, \buildrel\lt \over \sim \,\,}2$| larger or smaller than the original cloud, depending on its size and overdensity (Fig. 14). In fragmented streams, a power-law distribution of |$N(\gt m_{\rm cl})\propto m_{\rm cl}^{-1}$| develops over |$\sim (2-3)$| orders of magnitude in clump mass, consistent with previous studies (Gronke et al. 2022; Fielding et al. 2023). Such a scale-free distribution is emblematic of a much more general phenomenon known as Zipf’s law (Gronke et al. 2022). It requires a large dynamic range of clump sizes to develop, and has been argued to stem from the fact that for a highly fragmented cold-gas distribution, the cold-gas mass grows as |${\dot{m}}_{\rm cold,tot}\propto m_{\rm cold,tot}$| while the growth of the lowest mass clumps is dominated by mergers (Gronke et al. 2022).

At the high-mass end, we see deviations from the power-law behaviour as the distribution becomes dominated by a few large clumps along the stream axis onto which many of the intermediate-mass clumps have coagulated. The turnover at low masses is simply due to the resolution limit which sets a minimal mass for clumps. It may be affected in part by artificial disruption of small clumps in low-resolution regions far from the stream, though the same feature is present in simulations with uniform resolution as well, namely those with |$r_{\rm s,i}/\Delta =64$| and 128. Even in cases which coagulate, |$(r_{\rm s,i},{\chi _{\rm f}})=(3\, {\rm kpc},30),\, (30\, {\rm kpc},100)$|⁠, this power-law holds over a narrow range of intermediate masses where fragmentation is important. However, the global distribution in these cases is not well described by Zipf’s law since they are not highly fragmented and their distribution is predominanly shaped by coagulation.

Larger streams, with larger |$r^{*}/\ell _{\rm shatter}$|⁠, exhibit more violent fragmentation with many more clumps (compare the runs with |${\chi _{\rm f}}=1000$| and |$\, {r_{\rm s,i}}=3\, \mathrm{ and}\, 30\, {\rm kpc}$| in the insets). Comparing the curves with different resolutions in the main panels, we find that higher resolution runs produce more low-mass clumps indicative of stronger fragmentation. This is consistent with previous studies of the formation of multiphase media in other contexts (Sparre et al. 2019; Mandelker et al. 2021). However, the shape of the distribution is quite similar in all cases, developing a power law of |$N(\gt m)\propto m^{-1}$| over a wide range of clump masses, which turns over at the resolution scale, and not at some fixed scale despite |$\ell _{\rm shatter}$| being resolved.

However, whether the clump mass distribution is truly an example of Zipf’s law is not completely clear. We note that the full distribution down to the minimal clump mass is much better described by a lognormal distribution, rather than a power law with a turnover (see Appendix  D). The mode of the distribution is always at a clump mass slightly larger than the mass of a single cell at the cold density, even when |$\ell _{\rm shatter}$| is well resolved, supporting our conclusion that |$\ell _{\rm shatter}$| does not set a characteristic size for clumps.

The right-hand panel of Fig. 15 shows the cumulative distribution of clump radii. These are estimated as¹⁰|$r_{\rm cl}=[3m_{\rm cl}/(4\pi \rho _{\rm cl})]^{1/3}$|⁠, where |$\rho _{\rm cl}$| is the mean density in the clump, provided by the clump finder. The vertical dashed lines mark |$\ell _{\rm shatter}/r_{\rm s,i}$| for |$r_{\rm s,i}=1,\, 3,\, \mathrm{ and}\, 30\, {\rm kpc}$|⁠, from right to left. Focusing on simulations with |$r_{\rm s,i}=1\, \mathrm{ and}\, 3\, {\rm kpc}$| where |$\ell _{\rm shatter}$| is resolved, we see that there is no feature in the distribution of clump sizes at |$\ell _{\rm shatter}$| for fragmented streams. Rather, the distribution is a power law with |$N(\gt r_{\rm cl})\propto r_{\rm cl}^{-\beta }$| with |$\beta \sim (3-4)$|⁠, which breaks at the resolution limit. Considering that |$N\propto m_\mathrm{cl}^{-1}$|⁠, this indicates a mild size-dependence of clump density of |$\rho _\mathrm{cl}\propto r_\mathrm{cl}^{\beta -3}$|⁠, which is driven by smaller clumps containing more mixed gas at intermediate densities. As we increase the resolution, the break in the power law and the minimal clump size tend towards smaller sizes, rather than being tied to |$\ell _{\rm shatter}$|⁠. This suggests that clump sizes are not set by a hierarchical process where clumps continuously fragment to the local cooling-length, |$l_{\rm cool}\sim c_{\rm s}t_{\rm cool}$|⁠, until they reach |$\ell _{\rm shatter}={\rm min}(l_{\rm cool})$| as was originally proposed by McCourt et al. (2018, see also Das et al. 2021). Rather, while |$\, {r_{\rm s,i}}\gt \ell _{\rm shatter}$| may be a necessary condition for thermally unstable clouds to fragment, the actual formation of clumps is not itself driven by thermal instability. Rather, these are formed by a combination of RMI during the initial implosion and explosion processes (Gronke & Oh 2020b), and shredding (Jennings & Li 2021) and/or rotational fragmentation (‘splintering’) at late times (Farber & Gronke 2023). The clump sizes are thus determined by these processes and clumps can always exist at the smallest possible scales in a highly perturbed environment. Without additional processes such as thermal conduction (Koyama & Inutsuka 2004; Sharma, Parrish & Quataert 2010) or strong external turbulence (Tan et al. 2021; Gronke et al. 2022), we will always have clumps down to the grid scale when thermal fragmentation is present.

Consistent with the fact that |$\ell _{\rm shatter}$| is not a characteristic size for clumps, we also find that |$n_{\rm c}\ell _{\rm shatter}$| is not a characteristic column density for clumps, where |$n_{\rm c}$| is the number density in the cold phase. Rather, the distribution of column densities always peaks at the grid scale, |$N\sim n_{\rm c}\Delta$| (see Fig. D3). While this may seem to contradict the results of Liang & Remming (2020), we again note that |$\ell _{\rm shatter}$| was only marginally resolved in their simulations, so it is difficult to disentangle this from the grid scale, and moreover their simulation setup is sufficiently different from ours to render a detailed comparison difficult and beyond the scope of this paper.

Finally, we note that the power law describing the clump size distribution seems to break at |$r_{\rm cl}\sim r_{\rm s,i}/3\sim r_{\rm s,f}$|⁠. This seems to be the characteristic size of large clumps along the stream axis onto which many intermediate sized clumps have coagulated (Figs 6, 12, C2 and C4). A similar feature is seen in the cumulative distribution of clump sizes in sheets, where a power law is present from the resolution scale up to |$r_{\rm cl}\sim r_{\rm s,f}$| above which it flattens, decreasing again at |$r_{\rm cl}\lower.5ex\rm{\,\, \buildrel\gt \over \sim \,\,}r_{\rm s,i}$|⁠. The small clumps are visible in the plane of the sheet, within islands of hot gas in between large cold clumps (Figs 5 and C1). The large cold clumps surrounding these hot regions have typical sizes of order |$r_{\rm s,i}$| rather than |$r_{\rm s,f}$|⁠. In both streams and sheets, the size of the largest clumps seems rather constant in time, and is not merely a consequence of clump growth due to coagulation. The reasons for this are unclear, though further investigation of it is beyond the scope of this paper and is left for future work.

6 THE COAGULATION CRITERIA FOR STREAMS AND SPHERES

6.1 Fast versus slow coagulation

Our results in Sections 3 and 4 revealed a critical overdensity for sustained fragmentation in streams and spheres, |$\chi _{\rm f,crit}$|⁠, based on both the number of clumps, |$N_{\rm c}$|⁠, and their radial extent, |$d_\mathrm{max}$|⁠. Similar to previous work (Gronke & Oh 2020a; Farber & Gronke 2023), we find that |$\chi _{\rm f,crit}$| depends on the initial cloud size, or more specifically on |$r^{*}/\ell _{\rm shatter}$|⁠. In Fig. 16, we present a wide range of simulations with different initial sizes and metallicities and different values of |$\eta$| and |${\chi _{\rm f}}$|⁠, in both spherical and stream geometries (Table 1), in the plane of |$(r^{*}/\ell _{\rm shatter}, {\chi _{\rm f}})$|⁠. We distinguish those cases which remained fragmented from those that coagulated as in previous sections, and in Appendix B3, we show that this distinction is converged with resolution. We find that for both streams and spheres, |$\chi _{\rm f,crit}\propto (r^{*}/\ell _{\rm shatter})^{1/4}$|⁠, and that the normalization for streams is |$\sim 1.5$| times larger than for spheres. Below, we attempt to explain these trends and derive expressions for |$\chi _{\rm f,crit}$| in different geometries.

We recall that following the explosion (Section 3.1), the reflected shock reaches the interface between the cold and hot gas when the cloud radius is roughly |$r\sim r_{\rm s,f}$|⁠. At this point, linear RMI begins growing at the interface, seeded by perturbations in the collapsed cloud and at its interface, while on larger scales, the cloud continues expanding isotropically with its surface area growing as |$d_\mathrm{max}^n$|⁠, with |$n=1$| and 2 for streams and spheres, respectively. Thus, |$f_\mathrm{A}$| remains roughly constant during this phase, or may even increase due to perturbations and some fragmentation further increasing the surface area. Eventually, RMI grows to non-linear amplitudes characterized by familiar mushroom-shaped plumes and other non-isotropic features which, in the presence of cooling, quickly fragment into small clumps. The surface area then quickly becomes dominated by small clumps which do not further expand as they flow out, and the tight correlation between the total cloud surface area and |$d_{\rm max}$| breaks down. Once the gas has sufficiently fragmented, the total cold-gas surface area saturates and |$f_\mathrm{A}$| decreases as the fragmented clumps escape outwards, asymptotically approaching |$d_\mathrm{max}^{-n}$|⁠. Consequently, the coagulation time-scale increases significantly |$t_\mathrm{coag}\propto (d_\mathrm{max}/f_\mathrm{A})^{1/2} \propto d_\mathrm{max}^{(n+1)/2}$| (equation 22).

We can now formulate a condition for slow, inefficient, coagulation as follows. We assume that the fragmentation of the initial cloud into small clumps occurs at |$r\sim r_{\rm fA}$|⁠, the radius where |$f_{\rm A}$| begins monotonically declining following the explosion. If the cloud reaches this radius before significantly decelerating, then the clumps will continue expanding outwards with a velocity of |$\sim c_{\rm s,c}$| and will not coagulate, following the arguments presented in Section 3.2.1. If the cloud has significantly decelerated prior to fragmentation, then the resulting clumps will coagulate. In other words, for coagulation to become inefficient we demand

where |$t_\mathrm{decel,s}$| is the deceleration time-scale of the expanding cloud. If there is a turbulent mixing layer at the surface of the cloud, which is always the case in the presence of shape and/or density perturbations and in the absence of thermal conduction, then the cloud deceleration is dominated by condensation rather than ram pressure because the expansion velocity is |$v_{\rm ex}\sim c_{\rm s,c}$| (equation 16). We therefore have |$t_\mathrm{decel,s}=m_\mathrm{s,i}/\dot{m}$|⁠, where |$\dot{m}=4\pi r_\mathrm{avg}^2\rho _\mathrm{h}v_\mathrm{mix,s}$| for spheres and |$\dot{m}=2\pi r_\mathrm{avg}L\rho _\mathrm{h}v_\mathrm{mix,s}$| for streams, with |$v_\mathrm{mix,s}\sim 0.4c_\mathrm{s,c}(r_\mathrm{avg}/\ell _\mathrm{shatter})^{1/4}$| and |$r_\mathrm{avg}\sim r_\mathrm{fA}$|⁠.

All that remains is to estimate |$r_{\rm fA}$|⁠. In Fig. 17, we show |$f_{\rm A}(\lt d_{\rm max})$|⁠, similar to Figs 9 and 13, for a large number of simulations with a wide range of different properties in both spherical and stream geometries (see Table 1). Unlike previous figures, we here plot |$f_{\rm A}$| as a function of |$r/r_{\rm s,f}$| in order to distinguish cases with different |$\eta$|⁠. For this wide range of simulation parameters, we find that |$r_\mathrm{fA}\sim \alpha r_\mathrm{s,f}$|⁠, where |$\alpha \sim 6$| for streams and |$\alpha \sim 2$| for spheres. This can be more easily seen by examining the median profiles of |$f_{\rm A}$|⁠, which are shown with thick black lines. The larger |$\alpha$| for streams may be due in part to the Bell–Plesset effect (Bell 1951; Plesset 1954), which is a geometrical factor in the linear growth rate of RMI for converging or expanding shock waves, which can be |$\lower.5ex\rm{\,\, \buildrel\lt \over \sim \,\,}2$| times faster in spherical compared to cylindrical geometry. Additionally, if the total surface area saturates during the non-linear phase of RMI in both geometries, then |$f_{\rm A}\propto d_{\rm max}^{-n}$| declines faster in spherical geometry than in cylindrical geometry, yielding a larger |$r_{\rm fA}$| for the latter.

Inserting our expressions for |$t_\mathrm{decel,s}$| and |$r_\mathrm{fA}$| into the equation (30), we obtain a critical overdensity for sustained fragmentation,

for spheres, and

for streams, where |$\eta ^{*}\equiv \rho _\mathrm{s,f}/\rho (r=r^{*})$| and we have chosen a characteristic value of |$\alpha$| to set the normalization. These two relations are plotted in Fig. 16, ignoring the weak dependence on |$\eta ^{*}$|⁠, and are in very good agreement with the fragmentation thresholds for spheres and streams seen in the simulations. In fact, the characteristic values of |$\alpha$| in equations (31) and (32) were determined based on a fit to the threshold seen in Fig. 16, but are also nicely consistent with Fig. 17.

For our fiducial parameters presented in both Section 3 and the convergence tests in Appendix B3, namely |$r_{\rm s,i}=r^{*}=3\, {\rm kpc}$|⁠, |$\eta =\eta ^{*}=10$|⁠, |$Z=Z_{\odot }$|⁠, and no UVB, the critical overdensities for sustained fragmentation are |$\chi _{\rm f,crit}\sim 120$| for spheres and |$\sim 180$| for streams, consistent with our detailed results presented in those sections. We stress that the condition derived above for rapid coagulation does not explicitly depend on metallicity, or on the cooling rate more generally. This is because, as highlighted in Section 3.2.1, the condensation force dominates the deceleration process for all relevant parameters, including low-metallicity clouds/streams exposed to a UVB. However, the slower cooling for lower metallicity and/or in the presence of a UVB, and corresponding larger |$\ell _{\rm shatter}$| yield a lower critical overdensity for a given cloud/stream size.

To emphasize the transition between fast and slow coagulation, we plot in Fig. 18 the total time for outgoing clumps to turn around and coagulate, namely |$t_{\rm decel,s}+t_{\rm coag}$|⁠, as a function of |$d_{\rm max}$|⁠, for high-Z spheres and streams with |$\eta =10$| and |$r_{\rm s,i}=3\, {\rm kpc}$| (as in Section 3). These are all measured directly from the simulation, but |$d_{\rm max}$| here can be interpreted as the hypothetical maximal distance reached by a clump prior to turnaround, with the y-axis showing the total time to coagulation given this |$d_{\rm max}$|⁠. At small distances, this time-scale is dominated by |$t_{\rm decel}$| which is independent of d (equation 25). However, once the maximal distance exceeds |$r_{\rm max}\sim r_{\rm fA}$|⁠, the time-scale quickly becomes dominated by |$t_{\rm coag}\propto d_{\rm max}^{(n+1)/2}$| and coagulation becomes inefficient. Sheets only have a fast coagulation regime since |$f_\mathrm{A}$| never systematically decreases. Note that in this estimate of time-scales we neglect the continuous evolution of the clump velocity, which may slightly alter the actual time to turnaround.

As a final note, we stress that our distinction between ‘fast’ and ‘slow’ coagulation implies that given a long enough time and a large enough simulation volume, all cases should eventually coagulate. However, this is an artificial conclusion based on our highly idealized numerical setup. For example, external turbulence in the surrounding hot gas would drive the small clouds further away from the central cloud and/or outright disrupt them (Gronke et al. 2022). We thus expect that cases which are in our ‘slow coagulation’ regime will in practice not coagulate in realistic scenarios.

6.2 Comparison to previous work

Gronke & Oh (2020b) studied fragmentation versus coagulation of thermally unstable clouds using simulations similar to ours. Their initial conditions consisted of four identical spherical clouds, with radius |$r_{\rm s,i}$|⁠, overdensity |$\chi _{\rm i}$|⁠, thermal imbalance factor |$\eta$|⁠, and solar metallicity (without a UVB). They found a critical final overdensity for fragmentation of |$\chi _\mathrm{f}\gtrsim 300(r^{*}/\ell _\mathrm{shatter}/5000)^{1/6}$|⁠, while in our spherical simulations, we find |$\chi _\mathrm{f}\gtrsim 190(r^{*}/\ell _\mathrm{shatter}/5000)^{1/4}$| (equation 31, assuming a fixed |$\alpha$| and ignoring the weak dependence on |$\eta ^{*}$|⁠). These two criteria are compared in Fig. 16, where it is clear that our data are not well fitted by the Gronke & Oh (2020b) criterion. We predict a slightly lower |$\chi _\mathrm{f,crit}$| for |$r^{*}/\ell _\mathrm{shatter}\lesssim 10^{4}$|⁠, and a slightly steeper dependence on cloud sizes. Furthermore, their criterion depends only on the cloud radius when it loses sonic contact, |$r^{*}$|⁠, while ours has an additional dependence on the pressure contrast at this time, |$\eta ^{*}$|⁠. In practice, our criterion depends on the final cloud radius, |$r_{\rm s,f}=(\eta ^{*})^{-1/2}r^{*}$|⁠.

We suspect that the lower normalization is primarily due to our adopting a temperature floor of |$T_\mathrm{floor}=1.68\times 10^{4}\, \mathrm{K}$| compared to |$4\times 10^{4}\, \mathrm{K}$| in their simulations. Consequently, our |$\ell _\mathrm{shatter}$| is a factor of |$\sim 14$| smaller at the same density and metallicity, which results in a very similar normalization of |$\chi _\mathrm{f,crit}$| given the same initial cloud parameters. Moreover, since Gronke & Oh (2020b) initialized their simulations with four clouds compared to our one, they have four times the cold-gas mass given the same initial cloud parameters, which makes coagulation more likely. More difficult to explain is the different dependence of |$\chi _\mathrm{f,crit}$| on cloud size. We note that in fig. 5 of Farber & Gronke (2023), who performed similar experiments to Gronke & Oh (2020b), the slope of their atomic criterion is closer to |$1/4$| than |$1/6$|⁠. Regardless, these scalings are very weak, and other complications we have ignored (e.g. non-spherical geometry and background turbulence) could have a stronger impact on |${\chi _{\rm f}}$|⁠.

Gronke & Oh (2023) described an analogy between coagulation forces and gravity, summarized in Section 3.2.1 above (see equation 18). They used this to calculate the effective binding energy of a collection of clumps, and then evaluated a coagulation criterion by comparing this to the clumps’ kinetic energy, in analogy to the virial parameter of a self-gravitating system. While this analogy is compelling and has intriguing implications, it is intended to provide physical intuition rather than a precise quantitative criterion, as acknowledged in their paper. First, their model does not easily lend itself to calculating the dependence of |$\chi _{\rm f,crit}$| on |$r^{*}/\ell _{\rm shatter}$| without knowing the number of clumps, which is very difficult to model. Secondly, the analogy between coagulation forces and simple point-mass gravity may not strictly hold in our case. This analogy is expected to be valid when the velocities of clumps with respect to the central cloud are negligible compared to the hot background velocity, so that |$\Delta v_2\sim v_\mathrm{hot}\propto d^{-n}$| in equation (14). However, in the case of imploding/exploding thermally unstable clouds as studied in this work, the clumps escape at a velocity |$\sim c_\mathrm{s,c}$|⁠, which can be much larger than the hot gas velocity at large radii. Fig. 19 shows the radial velocities of cold clumps and hot gas near |$d_{\rm max}$| as a function of |$d_\mathrm{max}$|⁠, for high-Z spheres with |$\eta =10$| and |$r_{\rm s,i}=3\, {\rm kpc}$|⁠. For all |$\chi _\mathrm{f}$|⁠, hot gas velocities exhibit similar profiles of |$\propto d^{-2}$|⁠, while the velocities of cold clumps are roughly constant at |$c_\mathrm{s,c}$| during their escape. The larger relative velocity means that clumps feel a stronger coagulation force than if they were static, |$F=\dot{m}v$|⁠, facilitating ‘turnaround’ and eventual coagulation. Nevertheless, we agree with the qualitative conclusion of the Gronke & Oh (2023) energy argument that the coagulation efficiency increases from spheres to streams to sheets. A more detailed analysis of the effective virial parameter argument in fragmented systems would be interesting to pursue in future work.

7 COLD STREAMS PENETRATING VIRIAL SHOCKS

As described in Introduction, one of our main motivations for studying thermal fragmentation in low-Z streams is the application to cold streams feeding massive DM haloes at high-z from the cosmic web. In this section, we attempt to apply our results from previous sections to such streams, by constructing a toy model for how their interaction with virial shocks around DM haloes may cause them to fragment. The basic picture is that even though the cold streams themselves are not expected to shock as they enter the virial radius (Dekel & Birnboim 2006; Dekel et al. 2009), their confining pressure can increase by a large factor (Aung et al. 2024; Lu et al. 2024), and the resulting pressure contrast that develops between the hot CGM and the cold streams can cause streams to fragment. While this pressure contrast is not caused by radiative cooling of an intermediate-temperature stream as assumed in previous sections, we show in Appendix A2 below (Fig. A5) that, at least in sheets, the implosion and explosion processes are similar in cases with initially cold and underpressureized gas. In any event, this section serves as a proof of concept and a prelude to a more in depth study where we will explore this interaction using simulations (Yao et al., in preparation).

Cosmological simulations show that at high-z, cosmic-web filaments in the IGM have a three-zone structure (Lu et al. 2024). Their outer regions contain hot, diffuse gas in virial equilibrium within the potential well set by the DM filament. Interior to this is a zone of multiphase gas with high turbulence and vorticity. The innermost region is a dense, isothermal core. Lu et al. (2024) dubbed the outer two zones the circumfilamentary medium (CFM), while the innermost region is the cold stream that penetrates the hot CGM around massive galaxies. The relative size of each zone, and in particular how much of the CFM mass is hot at the virial temperature, depend on the profiles of the cooling and free-fall times within the filament (Birnboim, Padnos & Zinger 2016; Stern et al. 2021; Aung et al. 2024). Important for our purposes is that cosmological simulations suggest that these three zones are in approximate pressure equilibrium, with gas cooling isobarically from the hot CFM towards the cold stream (Ramsøy et al. 2021; Lu et al. 2024). Therefore, the pressure in cold streams prior to their entering the halo is roughly the virial pressure of cosmic-web filaments, while after entering the halo they become confined by the halo virial pressure. A schematic cartoon of this setup is presented in Fig. 20. We can thus estimate the pressure contrast between the cold stream and the hot CGM by comparing the virial pressure in cosmic-web filaments versus DM haloes.¹¹

The virial temperature and density of hot CGM gas are given by

(e.g. Dekel et al. 2013), where |$M_{12}=M_{\rm v}/10^{12}\, {\rm M}_\odot$| is the halo virial mass, |$(1+z)_3=(1+z)/3$| is the redshift, |$f_{\rm b}\sim 0.17$| is the Universal baryon fraction, |$\rho _{\rm u}(z)$| is the mean matter density in the Universe at redshift z, and |$\Delta _{\rm v,h}\sim 200$| is the halo virial overdensity in the spherical collapse model. In practice, the post-shock values of the temperature and density at the virial radius will differ from these values, which represent averages over the whole halo and neglect non-thermal pressure support of hot CGM gas (Komatsu & Seljak 2001; Lochhaas et al. 2021). Nevertheless, |$P_{\rm v,h}\propto \rho _{\rm v,h}\, T_{\rm v,h}$| using equations (33) and (34) is a good representation of the characteristic thermal pressure of hot CGM gas.

The corresponding values for filaments are (Lu et al. 2024),

Here, |$f_{\rm acc,3}=f_{\rm acc}/(1/3)$| is the fraction of total accretion onto the halo flowing along the given filament normalized to the typical value for a halo fed by three prominant streams (Dekel et al. 2009; Danovich et al. 2012), |$\mathcal {M}_{\rm v}=V_{\rm s}/V_{\rm v,h}$| is the inflow velocity along the stream normalized by the halo virial velocity, and |$\Delta _{\rm v,f}\sim (15-35)$| is the virial overdensity assuming cylindrical collapse (Mandelker et al. 2018; Lu et al. 2024).¹² We adopt |$\Delta _{\rm v,f}\sim 20$| as our fiducial value.

For completeness, we note that equation (35) stems from the fact that

where |$\Lambda _{\rm tot,f}$| is the total line mass of the DM filament,

where the specific accretion rate onto haloes is roughly (Neistein & Dekel 2008; Fakhouri, Ma & Boylan-Kolchin 2010; Dekel et al. 2013)

Combining equations (33)–(36), we obtain the pressure contrast between the hot CGM and the cold stream,

At |$z=4$| this results in |$\eta \sim 20$|⁠, which seems consistent with some cosmological simulations (Lu et al. 2024). In practice, neither the CGM nor the CFM are perfectly isobaric. The cold stream at the centre of the filament will have pressures slightly larger than |$P_{\rm v,f}$|⁠, while the hot CGM near the outskirts of the halo will have pressures slightly smaller than |$P_{\rm v,h}$|⁠. Lu et al. (2024) find the pressure to increase by a factor of |$\sim 2$| from the outer filament to the cold stream at the centre. Assuming a similar variation in the CGM pressure, |$\eta$| in equation (40) can be overestimated by a factor of |$\lower.5ex\rm{\,\, \buildrel\gt \over \sim \,\,}4$|⁠.

Gronke & Oh (2020b) found that clouds can thermally fragment if |$\eta \gt 2$|⁠, while for smaller values, they simply pulsate. Even if |$\eta$| is overestimated as described above, it is thus still expected to be large enough to cause cold streams to fragment upon entering the CGM of massive haloes, and this tendency is expected to increase towards lower redshifts. This may explain the large observed covering factions and clumping factors of cold gas in the CGM of |$M_{\rm v}\lower.5ex\rm{\,\, \buildrel\gt \over \sim \,\,}10^{12.5}\, {\rm M}_\odot$| haloes at |$z\lower.5ex\rm{\,\, \buildrel\gt \over \sim \,\,}3$| (Cantalupo et al. 2014, 2019; Borisova et al. 2016). It can also have important implications for gas accretion onto the central galaxy, and may lead to galaxy quenching by shutting off the cold-gas supply if the stream remains fragmented. Whether the stream remains fragmented or recoagulates depends on its size and overdensity compared to the critical overdensity for fragmentation.

Assuming pressure equilibrium between the cold stream and the hot CGM, their density ratio is (Mandelker et al. 2020b)

where |$\Theta _\mathrm{h}= T_\mathrm{h}/T_\mathrm{v}$| is the hot CGM temperature in units of the halo virial temperature, and |$\Theta _\mathrm{s}= T_\mathrm{s}/1.5\times 10^{4}\, {\rm K}$| is the stream temperature normalized to approximate thermal equilibrium with the UVB. For |$M_{\rm v}\gt 10^{12}\, {\rm M}_\odot$| haloes with a hot CGM, the radius of the stream in pressure equilibrium with the CGM is (Mandelker et al. 2020b)

with no explicit dependence on halo mass. Here, |$f_{\rm c,s}$| is the cold-gas mass fraction in the filament and |$f_{\rm h,0.3}$| is the hot gas mass fraction in the halo normalized to 0.3. Inserting equation (42) into equation (32), we obtain the critical overdensity for coagulation of cold streams,

where we have used |$(\eta ^{*})^{-1/2}r^{*}=r_{\rm s,f}$|⁠. Comparing this to equation (41), we find for cold streams

We thus find that cold streams can marginally sustain fragmentation for haloes with mass |$M_{\rm v}\lower.5ex\rm{\,\, \buildrel\gt \over \sim \,\,}10^{12}\, {\rm M}_\odot$| at redshifts |$z\lower.5ex\rm{\,\, \buildrel\gt \over \sim \,\,}2$|⁠, more so in more massive haloes and higher redshift. Coagulation becomes more likely towards lower redshift, which is further enhanced as the streams become more metal enriched thus decreasing |$\ell _{\rm shatter}$|⁠. On the other hand, |$\ell _{\rm shatter}\propto \rho ^{-1}\propto (1+z)^{-3}$| (McCourt et al. 2018; Mandelker et al. 2020b), which mitigates this effect somewhat.

We note that in more massive filaments, |$f_{\rm c,s}$| can be significantly lower than unity as larger virial temperatures and longer cooling times in the CFM increase the hot component of cosmic-web filaments at the expense of the cold streams (Aung et al. 2024; Lu et al. 2024). This becomes important in more massive haloes at a given redshift, and higher redshift for a given halo mass, and will increase |${\chi _{\rm f}}/\chi _{\rm f,crit}$| making fragmentation even more likely in these cases.

The final question we ask is whether sustained fragmentation will have time to manifest before streams reach the central galaxy. In our simulations, the full process consisting of the implosion, explosion, and fragmentation, takes |$\sim (1-2)t_{\rm sc}$| to manifest, where |$t_{\rm sc}$| is the sound crossing time of the initial stream radius (see Figs 4 and 12). The stream sound speed is (Mandelker et al. 2020b)

where |$\mu _\mathrm{s,0.6}=\mu _\mathrm{s}/0.6$| is the mean molecular weight. The uncertainties in |$\eta$| from equation (40), as well as our having ignored self-gravity and stream rotation, both of which are important for setting the size of streams in the IGM (Lu et al. 2024), make the initial stream radius prior to contraction difficult to constrain. We thus take as a strict lower limit the final stream radius, |$R_{\rm s,f}$| from equation (42), to obtain a lower limit for the stream sound crossing time near |$R_{\rm v}$| and thus the time-scale for fragmentation,

The streams reach the central galaxy in roughly a virial crossing time, |$t_{\rm v}\sim R_{\rm v}/V_{\rm v}$|⁠, which is given by (e.g. Mandelker et al. 2020b)

We thus have

This ratio decreases with an increasing z making stream fragmentation more likely to manifest in the CGM at lower redshift, though still at |$z\lower.5ex\rm{\,\, \buildrel\gt \over \sim \,\,}2$| due to equation (44). However, since the ratio is typically less than 1 it is unclear if streams will have time to fragment prior to reaching the central galaxy. Before we neglect this possibility, we list several issues that may help. First, the stream sound crossing time decreases rapidly as it flows towards the central galaxy due to its shrinking radius, while the stream velocity increases by a factor of |$\lower.5ex\rm{\,\, \buildrel\lt \over \sim \,\,}2$| from |$R_\mathrm{v}$| to |$0.1\, R_\mathrm{v}$| (Aung et al. 2024). As a result, the sound crossing time decreases faster than the inflow time, and at some point they should become comparable. Furthermore, the virial shock can extend up to |$\sim 2\, R_\mathrm{v}$| (Zinger et al. 2018; Aung, Nagai & Lau 2021). Finally, as shown in fig. 1 of Mandelker et al. (2020b), there is an order of magnitude uncertainty in |$R_\mathrm{s,f}$| due to the various additional parameters in the equations. Thinner streams, which also tend to be denser, will be more prone to fragmentation and disruption. Considering other physical mechanism neglected here, such as the halo potential, a stratified CGM, turbulence and shear, the dynamics of fragmentation and the fate of streams as they penetrate the CGM remain unclear. This will be explored in detail in future work.

8 CAVEATS AND ADDITIONAL PHYSICS

While our analysis has been thorough in terms of the impact of cloud geometry and metallicity on thermal fragmentation and coagulation, we have neglected a number of important physical processes which can impact these issues, as well as our application to cold streams penetrating the hot CGM. We briefly discuss these here, leaving more in depth analysis to future work.

• Turbulent environments. In all our simulations, the background was initially static. However, the CGM and the high-z cosmic web are both highly turbulent environments, driven by inflows, outflows, and galaxy interactions. Turbulence plays a critical role in the dynamics of multiphase gas, sustaining turbulent radiative mixing layers around clumps that facilitate condensation and coagulation, while also contributing to shattering of large clouds and disruption of small clumps. The latter one is particularly significant because the clump velocities in our simulations, |$\sim c_\mathrm{s,c}$|⁠, are much smaller than the expected turbulent velocities (e.g. Mandelker et al. 2021; Gronke et al. 2022; Fielding et al. 2023; Das & Gronke 2024). Momentum coupling between the cold clumps and the turbulent hot environment quickly entrains the clumps in the turbulent flows, if it does not first destroy them due to mixing. This supports our claim from Section 6.1 that our ‘slow coagulation’ regime is likely an idealization. Such clouds will likely either remain fragmented with the resulting clumps highly dispersed and with a minimal size set by the turbulence (Gronke et al. 2022), or else the resulting small clumps will all be destroyed making it a matter of definition whether we consider this cloud fragmented or not. Turbulence may also impact the critical overdensity for fragmentation, reducing the efficiency of coagulation.

• Magnetic fields. Most astrophysical plasma is magnetized. While the outer CGM and cosmic web at high-z are not expected to be highly magnetized, the magnetic fields can amplify during cloud contraction and fragmentation, and small-scale cloudlets resulting from thermal fragmentation can therefore be magnetically dominated and out of thermal pressure equilibrium with their surroundings (Nelson et al. 2020). This non-thermal pressure support leads to streaming motions along field lines and filamentary clumps, as the total pressure gradient becomes unbalanced along the field direction (Wang, Oh & Jiang, in preparation). Magnetic draping affects the kinematics and morphology of clumps (Hidalgo-Pineda, Farber & Gronke 2024; Ramesh et al. 2024). Magnetic fields also suppress mixing in the turbulent mixing layer, reducing the entrainment rate (Ji, Oh & Masterson 2019; Sparre et al. 2020; Gronke & Oh 2020a; Grønnow et al. 2022). However, if turbulence is externally driven, then for a given rms turbulent velocity hydrodynamic and MHD mixing rates are similar (Das & Gronke 2024).

• Cosmic rays. Similar to magnetic field, cosmic rays are ubiquitous and play a crucial role in the CGM. Their relatively long cooling times and strong dynamical coupling to the gas make them a significant factor in galaxy formation processes (Ruszkowski & Pfrommer 2023). Cosmic rays provide substantial non-thermal pressure in the CGM, helping to counteract gravity and prolong the retention of cold gas. This, in turn, supports the survival of fragmented clumps and increases the cold mass fraction in the CGM (Butsky et al. 2020; Ji et al. 2020).

• Thermal conduction. Thermal conduction is a universal dissipation process that smooths temperature gradients along magnetic field lines, leading to the evaporation of clouds smaller than the field length when competing with radiative cooling (Brüggen & Scannapieco 2016; Armillotta et al. 2017; Li et al. 2020). As a result, the field length can act as a lower limit for clump size, similar to turbulence, and may change the extent of small-scale clumps resulting from thermal fragmentation.

• Self-gravity. Self-gravity affects the dynamics of cold clouds in hot gas only when the cloud mass exceed the Jeans mass (Li et al. 2020), which can occur in large-scale cosmic-web structures. It increases the binding energy of explosive clouds, thereby potentially inhibiting sustained fragmentation and promoting coagulation in these larger clouds. More over, strong implosions like those we see prior to fragmentation, could directly trigger star formation, offering insights into star formation processes far from galaxies. Self-gravity is also important for cold streams (Mandelker et al. 2018; Lu et al. 2024), and contributes to their radial structure as well as to their interactions with the CGM (Aung et al. 2019).

• Shear. Shear plays a critical role in cold streams penetrating virial shocks. The shear between the cold stream and the hot background gas can drive Kelvin–Helmholtz instabilities and enhance mixing. This will either promote stream growth through entrainment of hot gas in the mixing layer, or lead to stream disruption, depending on the ambient conditions (Mandelker et al. 2020a; Aung et al. 2024; Ledos, Takasao & Nagamine 2024). Additionally, shear can induce streaming motions in explosive clumps, as well as induce additional fragmentation within the turbulent mixing layer, leading to the formation of filamentary clump structures. However, the effect of shear on the fragmentation and coagulation processes, particularly in a stratified medium as found in the hot CGM in galaxy haloes, is unknown.

9 SUMMARY AND CONCLUSIONS

We studied the effects of cloud geometry and metallicity on the thermal fragmentation and coagulation of cold-gas clouds in the CGM or the high-z cosmic web by utilizing 3D idealized hydro simulations. We initialized a thermally unstable warm cloud embedded in a hot medium in pressure equilibrium. We explored three different cloud geometries (planar sheets, cylindrical streams, and spherical clouds), varying the final overdensity of the cloud with respect to the background after thermal and pressure equilibrium have been reestablished, |${\chi _{\rm f}}\equiv \rho _{\rm s,f}/\rho _{\rm bg}$|⁠, with |$\rho _{\rm s,f}$| the final density of the cold cloud (sheet, stream, or sphere) and |$\rho _{\rm bg}$| the density of hot background. We also varied the size of the initial cloud, |$r_{\rm s,i}$|⁠, and the ratio of the initial (warm) cloud density to the final (cold) cloud density after thermal and pressure equilibrium have been re-established, |$\eta \equiv \rho _{\rm s,f}/\rho _{\rm s,i}$|⁠, which also represents the ratio of the background pressure to the pressure in the cloud after thermal equilibrium. Initially, we explored a setup where both the cloud and background have solar metallicity and there is no ionizing UV background (our ‘high-Z’ case, Section 3). We then extended our study of streams specifically to the low-metallicity regime where the streams and background have metallicities |$Z_{\rm s}=0.03Z_{\odot }$| and |$Z_{\rm bg}=0.1Z_{\odot }$|⁠, respectively, and both are exposed to a |$z=2$| (Haardt & Madau 1996) UVB (our ‘low-Z’ case, Section 4). We examine the distribution of clump sizes in a regime where the shattering length-scale, |$\ell _{\rm shatter}\equiv {\rm min}(c_{\rm s}t_{\rm cool})$|⁠, is well resolved to explore its importance as a characteristic size for cold gas (Section 5) and derived analytical criteria for the fragmented cloudlets to recoagulate (Section 6). Finally, we applied our results for low-Z streams to the case of cold streams feeding massive haloes (⁠|$M_{\rm v}\lower.5ex\rm{\,\, \buildrel\gt \over \sim \,\,}10^{12}\, {\rm M}_\odot$|⁠) at high redshift (⁠|$z\lower.5ex\rm{\,\, \buildrel\gt \over \sim \,\,}2$|⁠) from the cosmic web through virial accretion shocks (Section 7).

The general evolution of such thermally unstable clouds consists of two stages: (1) the initial implosion and explosion processes, triggered by a rapid loss of thermal pressure in the cloud followed by a converging and subsequent reflected shock that induces fragmentation of cold gas through a combination of RMI, strong cooling, and vorticity induced by surface perturbations; and (2) the deceleration and recoagulation processes as escaping clumps begin comoving with the entrainment flow of hot gas converging onto the central cloud. We schematically summarize our main results and key insights into these two stages in the cartoon presented in Fig. 21. In detail, our results can be summarized as follows:

• The explosion velocity, which determines the growth time-scale of RMI, the fragmentation time-scale of the initial cloud, and (eventually) the spatial distribution of fragmented clumps, is close to the sound speed of cold gas for spheres and streams, regardless of metallicity, while for sheets it is |$\sim 0.4\, c_\mathrm{s,c}$| (Figs 2 and A5).

• Due to this relatively low explosion velocity (⁠|$v_{\rm ex}\sim c_\mathrm{s,c}$|⁠), the dominant deceleration mechanism for escaping clumps is always momentum exchange with the hot gas through condensation rather than drag or ram pressure. This is true regardless of initial cloud geometry and metallicity (equation 16).

• We distinguish cases where clouds remain fragmented from cases where they recoagulate based on both the number of clumps and their maximal radial distance from the centre, |$d_{\rm max}$|⁠. Based on these metrics, both streams and spheres display sustained fragmentation only above a critical final overdensity (⁠|${\chi _{\rm f}}\lower.5ex\rm{\,\, \buildrel\gt \over \sim \,\,}200$| at high-Z). Sheets, on the other hand, show coagulation at all |${\chi _{\rm f}}$| (Figs 4 and 8).

• At the end of their evolution, spherical clouds either remain fragmented into many small clumps or else coagulate into a single cloud at the centre. This is not the case for streams and sheets. Rather, even if such clouds coagulate radially, coagulation is suppressed along the stream axis and within the plane of the sheet. The final number of clumps thus increases with |$\chi _\mathrm{f}$|⁠, even if the farthest clump falls back to the central plane/line (Figs 5–7 and C1–C3).

• The coagulation time-scale is longest for spheres and shortest for sheets. This is determined primarily by the area modulation factor, |$f_{\rm A}$|⁠, which is the ratio of the total surface area of cold gas to the area of the cloud-centric shell at the distance of the furthest clump (equation 23). This decreases with the distance of escaping clumps as |$d^{-n}$|⁠, where |$n=0,1,2$| for sheets, streams, and spheres, respectively (Fig. 9).

• Low-Z streams behave similarly to high-Z streams with a similar ratio of |${\chi _{\rm f}}/v_\mathrm{mix}$|⁠, which determines the deceleration time-scale (equation 25). For our chosen metallicities and UVB, this translates into a factor |$\sim 3$| increase in |${\chi _{\rm f}}$| from low-Z to high-Z. This similar behaviour is seen in the morphology, number of clumps, maximum clump distance, and total cold-gas area (Figs 11–13).

• As |${\chi _{\rm f}}$| increases, the fragmented clumps can enter the ‘disruption’ regime and proceed to mix with the hot background, thus decreasing the total cold-gas mass. For our low-Z runs, this occurs at |${\chi _{\rm f}}\sim 1000$|⁠, though for our high-Z runs clumps survive at all simulated values of |${\chi _{\rm f}}$| (Fig. 14 and equation 29).

• For fragmented clouds, the cumulative clump mass distribution, |$N(\gt m)$|⁠, closely follows a power law with an index of |$\sim (-1)$|⁠, in accordance with Zipf’s law, over |$\sim (2-3)$| dex in clump mass (Fig. 15, left). However, the full distribution of clump masses down to the smallest and largest scales is better fit by a lognormal distribution (Appendix  D).

• This power law extends from roughly the resolution scale, even if this is much smaller than |$\ell _{\rm shatter}$|⁠, until |$r_{\rm s,f}$|⁠, the size the initial cloud would have had assuming monolithic collapse to thermal and pressure equilibrium. The distribution of clump radii is a power law over the same range, showing that |$\ell _{\rm shatter}$| is not a characteristic size for cold-gas clumps resulting from thermal-instability-induced fragmentation (Fig. 15, right).

• Exploring a wide range of simulations in both spherical and stream geometries, we find a critical overdensity for sustained fragmentation in streams and spheres that scales as |$\chi _{\rm f,crit}\propto (r_{\rm s,f}/\ell _{\rm shatter})^{1/4}$| (Fig. 16). For |$\chi _{\rm f}\lt \chi _{\rm f,crit}$| the fragmented cloudlets recoagulate. The normalization of |$\chi _{\rm f,crit}$| for streams is larger than for spheres by a factor of |$\sim 1.5$| (Fig. 16), while sheets always recoagulate.

• Our model for sustained fragmentation differs from previous models (e.g. Gronke & Oh 2020b), where |$\chi _{\rm f,crit}$| depends on the radius where the cloud loses sonic contact, |$r^{*}$|⁠, rather than |$r_{\rm s,f}$|⁠.

• We propose a model for ‘slow’ versus ‘fast’ coagulation, based on the competition between fragmentation of the initial cloud and deceleration and entrainment of the clumps, which agrees well with our empirical coagulation criterion. Sheets are always in the fast coagulation regime (Section 6.1).

• Many of these conclusions may change when additional physics are included, such as externally driven turbulence, magnetic fields, thermal conduction, self-gravity, or shear between the initial cloud and the background. These effects should be explored in future work.

The application of our model to cold streams penetrating the hot CGM of massive haloes at high-z from the cosmic web (Section 7) rests on the realization that the confining pressure of these streams in intergalactic filaments is lower than their confining pressure in the CGM (Lu et al. 2024). This could cause streams to fragment as they enter the virial shock, in much the same way as clouds which are underpressurized due to cooling fragment in our simulations. We evaluate this pressure contrast (⁠|$\eta$|⁠), the final density contrast between the cold streams and hot CGM (⁠|${\chi _{\rm f}}$|⁠), and the critical overdensity for fragmentation (⁠|$\chi _{\rm f,crit}$|⁠) as a function of halo mass and redshift, based on previous models for the properties of cold streams in the CGM (Mandelker et al. 2020b). We find that fragmentation can be important for haloes with |$M_{\rm v}\gt 10^{12}\, {\rm M}_\odot$| at |$z\gt 2$|⁠. By comparing the fragmentation time-scale to the inflow time of the streams, we find that fragmentation is more likely to manifest towards the lower end of this redshift range. This offers a possible explanation for the large clumping factors and covering fractions of cold gas in the CGM around such galaxies, and may be related to galaxy quenching by providing a mechanism to prevent cold streams from reaching the central galaxy. These issues will be explored using direct simulations in upcoming work (Yao et al., in preparation).

ACKNOWLEDGEMENTS

We thank the anonymous referee for helpful comments that have improved the quality of this manuscript. We thank Romain Teyssier for technical support and help with the ramses clump finder. We thank Andrei Antipov, Yuval Birnboim, Frank Van den Bosch, Orly Gnat, Max Gronke, Elisha Modelevsky, and Daisuke Nagai for helpful comments and interesting discussions. ZY and NM acknowledge support from BSF grants 2020302 and 2022281 and NSF-BSF grant 2022736. SPO thanks NASA grant 19-ATP19-0205 and NSF grant 240752 for support. AD has been partly supported by the grants ISF861/20, NSF-BSF 2023730, and NSF-BSF 2023723. The simulations were performed on the Moriah cluster at the Hebrew University.

Software:matplotlib (Hunter 2007), numpy (Harris et al. 2020), astropy (Astropy Collaboration 2022), yT (Turk et al. 2011), scikit-image (Walt et al. 2014), and powerlaw (Alstott, Bullmore & Plenz 2014).

DATA AVAILABILITY

The data supporting the plots within this article are available on reasonable request to the corresponding author.

Footnotes

REFERENCES

APPENDIX A: THE IMPLOSION AND EXPLOSION PROCESSES IN SHEETS

Complementary to the qualitative discussion regarding the implosion and explosion processes in all three geometries presented in Section 3.1.1, we here offer a rigorous mathematical description of these processes in sheets.

A1 Adiabatic underpressurized sheets

Consider a cold, underpressurized sheet of density |$\rho _{\rm c}$| and pressure |$P_{\rm c}$|⁠, surrounded by hot gas with density |$\rho _{\rm h}=\rho _{\rm c}/\chi _{\rm i}$| and pressure |$P_{\rm h}=\eta P_{\rm c}$|⁠. This setup, shown by the black lines in Fig. A1 for the case |$\chi _{\rm i}=\eta =10$|⁠, is unstable, as the negative pressure gradient drives the contraction of the sheet towards its centre. The evolution during the contraction, shown in Fig. A1 by the blue lines, is similar to the Sod shock tube.

Inside the sheet (the left-hand side of the profiles in Fig. A1), a shock wave develops according to the Rankine–Hugoniot conditions, describing the conservation of mass, momentum, and energy across the shock:

Here |$v_\mathrm{L}=0$|⁠, |$\rho _\mathrm{L}=\rho _\mathrm{c}$|⁠, and |$P_\mathrm{L}=P_\mathrm{c}$| are the physical quantities in the sheet (on the left-hand side). The velocities in equations (A1)–(A3) are in the lab frame, so we subtract the shock velocity, |$v_{\rm sh}$|⁠, to describe the equations in the comoving shock frame. The star symbol, |$*$|⁠, represents the region that develops between the initial cold and hot media (between the left- and right-hand sides of the initial condition). The density and temperature in the starred region have sudden transitions across the CD, while the velocity and pressure are constant. Therefore, |$\rho _\mathrm{*,L}$| denotes the density on the left-hand side of the starred region, while |$v_*$| and |$P_*$| represent the velocity and pressure throughout the whole starred region. These three equations can be combined to eliminate |$v_{\rm sh}$| and obtain expressions for |$v_{\rm *}$| and |$\rho _{\rm *,L}$|⁠,

We find that if |$P_*/P_\mathrm{L}\gg 1$| and |$v_{\rm L} \ll c_{\rm s,L}$|⁠, then |$v_*\propto c_\mathrm{s,L}\sqrt{P_*/P_\mathrm{L}}$| and |$\rho _\mathrm{*,L}\sim \rho _\mathrm{L}(\gamma +1)/(\gamma -1)$|⁠.

In the background (the right-hand side of the profiles in Fig. A1), adiabatic expansion produces a rarefaction wave that conserves entropy and the generalized Riemann invariant:

Here |$v_\mathrm{R}=0$|⁠, |$\rho _\mathrm{R}=\rho _\mathrm{h}$|⁠, and |$P_\mathrm{R}=P_\mathrm{h}$| are the physical quantities in the background (on the right-hand side). |$\rho _\mathrm{*,R}$| and |$c_\mathrm{s,*}$| are the density and the sound speed on the right-hand side of the starred region, respectively. These equations yield

Equations (A4), (A5), (A8), and (A9) offer solutions for the implosion quantities, |$v_*=v_\mathrm{im}\lt 0$|⁠, |$P_*=P_\mathrm{im}$|⁠, |$\rho _\mathrm{*,L}=\rho _\mathrm{im,c}$|⁠, and |$\rho _\mathrm{*,R}=\rho _\mathrm{im,h}$|⁠. Furthermore, from equation (A7), we get

and subsequently

Finally, the shock velocity |$v_\mathrm{sh}$| can be derived from equation (A1)

These solutions for the implosion quantities are shown by the blue lines in Fig. A2. We show these as functions of |$\eta$|⁠, while fixing |$\chi _\mathrm{f}\equiv \eta \chi _\mathrm{i}=100$|⁠. The blue dots show results from 3D sheet simulations for comparison. We find that they agree with each other very well.

The implosion shock propagates towards the centre of the sheet where it collides with the shock from the other side at |$t_{\mathrm{coll}}=r_{\mathrm{s,i}}/v_{\mathrm{sh}}$|⁠. This head-on shock collision, shown in Fig. A1 by the orange lines, is similar to the planer Noh problem. Applying the shock solutions from equations (A4) and (A5) to both sides, and inserting |$v_*=0$|⁠, |$v_\mathrm{R}=-v_\mathrm{L}=v_\mathrm{im}\lt 0$|⁠, |$\rho _\mathrm{L}=\rho _\mathrm{R}=\rho _\mathrm{im,c}$|⁠, |$P_\mathrm{L}=P_\mathrm{R}=P_\mathrm{im}$|⁠, and |$c_\mathrm{s,L}=c_\mathrm{s,R}=c_\mathrm{s,im,c}$| from symmetry arguments, we can derive expressions for |$P_*=P_\mathrm{coll}$|⁠, |$\rho _*=\rho _\mathrm{coll}$|⁠, and |$c_\mathrm{s,*}=c_\mathrm{s,coll}$|⁠. We obtain

where |$\mathcal {M}_{\rm im,c}\equiv v_\mathrm{im}/c_\mathrm{s,im,c}$|⁠. For |$\mathcal {M_{\rm im,c}}\gg 1$| and |$\gamma =5/3$|⁠, we have |$P_\mathrm{coll}/P_\mathrm{im} \sim \mathcal {M}_{\rm im,c}^2$| and |$\rho _\mathrm{coll}/\rho _\mathrm{im,c} \sim 4$| as expected for strong shocks. Together these imply that |$c_\mathrm{s,coll}\sim v_\mathrm{im}/2$|⁠. The orange lines and points in Fig. A2 compare these solutions with simulation results and find excellent agreement.

Once all the cold gas has been shock-heated, the inflowing hot gas is unable to prevent the expansion of shocked cold gas. At this stage, a shock wave develops on the right-hand side, propagating into the background, while a rarefaction waveforms on the left-hand side, propagating into the cloud. This stage is shown by the red lines in Fig. A1. Following similar steps to those applied above when deriving the implosion quantities, but switching the left- and right-hand sides and using the conditions of |$v_\mathrm{L}=0$|⁠, |$v_\mathrm{R}=v_\mathrm{im}$|⁠, |$c_\mathrm{s,L}=c_\mathrm{s,coll}$|⁠, |$c_\mathrm{s,R}=c_\mathrm{s,im,h}$|⁠, |$P_\mathrm{L}=P_\mathrm{coll}$|⁠, and |$P_\mathrm{R}=P_\mathrm{im}$|⁠, we can obtain expressions for |$P_*=P_\mathrm{ex}$|⁠, |$v_*=v_\mathrm{ex}$|⁠, |$\rho _\mathrm{*,L}=\rho _\mathrm{ex,c}$|⁠, and |$\rho _\mathrm{*,R}=\rho _\mathrm{ex,h}$|⁠. During this phase, the reversal of the direction of the rarefaction wave modifies equation (A8) so that we have

This equation provides a maximum value for the explosion velocity of |$v_{\rm ex,max}=2c_\mathrm{s,coll}/(\gamma -1)$|⁠. In practice, this is always an overestimate because |$P_*$| is bounded from below by |$P_{\rm im}$|⁠, and furthermore the exponent |$(\gamma -1)/(2\gamma)= 0.2$| is weak. The red lines and points in Fig. A2 compare our model predictions for the explosion quantities with simulation results, again finding excellent agreement.

Overall, we find that our model for the implosion and explosion properties predicts the simulation results quite well. Both the implosion and explosion velocities increase with |$\eta$|⁠, while the latter velocity is smaller than the former because part of the kinetic energy is converted to internal energy when the shock heats the cold gas during the central collision, increasing its temperature.

A2 Underpressurized sheets with radiative cooling

We expand upon the simple adiabatic problem explored above by including radiative cooling, which clearly plays an important role in the fragmentation process we study in the main text. We consider identical initial conditions to those explored in Appendix A1, and assume that both the cold and the hot gas have solar metallicity. We set a temperature floor at the temperature of initial cold gas, |$T_\mathrm{c}\sim 10^4~\mathrm{K}$|⁠, and prohibit radiative cooling above |$0.8~T_\mathrm{h}$| as in the main text. Similar to the adiabatic case, the initial implosion has a shock solution on the left and a rarefaction wave on the right. However, since radiative cooling is efficient in the cold and dense in the sheet, the post-shock gas cools down immediately after being heated up, effectively forming an isothermal shock, as shown by the blue lines in Fig. A3.

A2.1 Isothermal shocks

In isothermal shocks, energy conservation (equation A3) is replaced by

where |$c_\mathrm{s,L}$| is the adiabatic¹³ sound speed in the left-hand fluid. Consequently, the velocity, density, and sound speed in the left-hand starred region are given by

Since radiative cooling is prohibited in the hot gas, the rarefaction wave on the right-hand side has a similar solution to the adiabatic case in equations (A8) and (A9). The four equations (A8), (A9), (A17), and (A18), thus provide the solutions for the implosion quantities in the case of an isothermal shock |$v_*=v_\mathrm{im,iso}$|⁠, |$P_*=P_\mathrm{im,iso}$|⁠, |$\rho _\mathrm{*,L}=\rho _\mathrm{im,iso,c}$|⁠, |$\rho _\mathrm{*,R}=\rho _\mathrm{im,iso,h}$|⁠, and |$c_\mathrm{s*,L}=c_\mathrm{s,floor}$|⁠. Comparing |$v_\mathrm{im,iso}$| in equation (A17) to |$v_\mathrm{im}$| in equation (A4), we find that for |$P_*/P_\mathrm{L}\gg 1$|⁠, the isothermal implosion velocity is roughly a factor of |$\sqrt{(\gamma +1)/2}\sim 1.15$| larger than the adiabatic one.

The head-on shock collision is also isothermal, shown by the orange lines in Fig. A3. By applying equations (A17) and (A18) on both sides, we have |$P_\mathrm{coll,iso}/P_\mathrm{im,iso}\sim \rho _\mathrm{coll,iso}/\rho _\mathrm{im,iso,c}\sim \gamma (v_\mathrm{im,iso}/c_\mathrm{s,floor})^2$| and |$c_\mathrm{s,coll}=c_\mathrm{s,floor}$|⁠.

However, the expansion phase is adiabatic for both the shock propagating into the background on the right and the rarefaction wave propagating into the cloud on the left, shown by the red lines in Fig. A3. This is because the shock heats up the hot gas in the background (on the right) to the regime where radiative cooling is prohibited, while the cold gas in the cloud (on the left) is already at the temperature floor and cannot cool, though its temperature can decrease through adiabatic expansion. Therefore, we expect similar solutions to the adiabatic case, and in particular the maximum explosion velocity is predicted to be |$v_{\rm ex,max}\sim 2\, c_\mathrm{s,floor}/(\gamma -1)\sim 3\, c_{\rm s,floor}$| from equation (A15), larger by a factor of |$\sim 2$| than our estimate of |$v_{\rm ex}\sim 1.34\, c_{\rm s,c}$| from equation (10). However, as noted above, the maximum limit provided by equation (A15) is necessarily an overestimate, while the limit provided by equation (10) is not.

A2.2 Radiative cooling in mixing layers

In addition to isothermal shocks, radiative cooling also takes place at the interfaces between cold and hot phases, where gas is at intermediate temperatures and consequently has efficient cooling. In the context of our shock model, we assume this occurs at the CD. Radiative cooling causes enthalpy loss and pulls both cold and hot gas into this transition region, inducing turbulent mixing that smooths the CD. Therefore, we hereafter refer to this region as the ‘mixing layer’ instead of the CD. The effects of radiative cooling in the mixing layer can be seen in the blues lines in the pressure and velocity panels of Fig. A3. There is a small dip of pressure at the mixing layer due to cooling. This pressure difference pulls both hot and cold gas into the mixing layer, causing the hot gas to accelerate and the cold gas to decelerate with respect to their implosion velocities in the lab frame. During the explosion phase, the cold gas is accelerated and the hot gas is decelerated near the mixing layer, as can be seen by the red line in the velocity panel.

To study this process quantitatively, we begin with the momentum and energy equations 

where |$\epsilon$| is the specific internal energy. The terms on the right-hand side of equation (A21) are (from left to right) the advection of thermal energy, adiabatic expansion/compression, radiative cooling, and viscous dissipation (caused by turbulence, viscosity, etc.). In steady state, the thermal pressure is roughly the same on either side of the mixing layer so the ram pressure exerted by cold and hot gas inflowing into the mixing layer should roughly balance,

This can also be deduced from equation (A20). Here, |$v_\mathrm{ent,c}$| and |$v_\mathrm{ent,h}$| are the entrainment velocities of cold and hot gas into the mixing layer, respectively. Assuming that the advection of thermal energy into the mixing layer from the hot phase is balanced by the advection of thermal energy out of the mixing layer into the cold phase and that viscous dissipation is negligible, then radiative cooling is only balanced by adiabatic compression (equation A21),

where H is the width of the mixing layer. Combining equations (A22) and (A23), we can derive the entrainment velocities of cold and hot gas with respect to the mixing layer

where |$\chi =\rho _\mathrm{c}/\rho _\mathrm{h}$| is the density ratio between the cold and hot phases.

We can now evaluate the velocity of cold and hot gas in the lab frame. Both phases share the same implosion velocity |$v_\mathrm{im,iso}$|⁠, and we assume the mixing layer has a velocity |$v_{\rm ml}$| relative to |$v_\mathrm{im,iso}$|⁠. If radiative cooling is efficient, the mixing layer moves outwards into the hot phase by accumulating cold gas due to the cooling of intermediate-temperature gas in the mixing layer, |$v_{\rm ml}\gt 0$|⁠. On the other hand, if radiative cooling is inefficient, the mixing layer moves inwards into the cold phase due to the heating by mixing, |$v_{\rm ml}\lt 0$|⁠. Thus, in the lab frame, the cold-gas velocity during implosion is |$v_\mathrm{im,c}=v_\mathrm{im,iso}+v_\mathrm{ml}+v_\mathrm{ent,c}$|⁠, while the hot gas velocity is |$v_\mathrm{im,h}=v_\mathrm{im,iso}+v_\mathrm{ml}+v_\mathrm{ent,h}$|⁠. Therefore, the velocity difference is

which only depends on the properties of the mixing layer.

In the simulation shown in Fig. A3, there is an initial pressure jump at the boundary of cold and hot phases, limited by the grid scale, |$\Delta$|⁠. The width of the mixing layer and the velocity difference thus initially satisfy |$|v_{\rm im,c}-v_{\rm im,h}|\propto H\sim \Delta$|⁠. However, we find that the velocity difference in this case is proportional to |$\Delta ^{1/2}$| rather than |$\Delta$| (blue dots in Fig. A4). This is consistent with Tan et al. (2021), who found that the surface brightness, which is proportional to the entrainment velocity, scales as |$\Delta ^{1/2}$| in the absence of thermal conduction.

The result |$v_\mathrm{ent}\propto H$| was first proposed by Begelman & Fabian (1990), but was later found to not hold by recent 3D turbulent mixing layer simulations (Ji et al. 2019; Tan et al. 2021). The fundamental reason is that viscous dissipation cannot be neglected. In the absence of strong turbulence, the dissipation on the grid scale (the scale of the initial phase transition) is mainly numerical. According to the Navier–Stokes equations, the viscous dissipation rate can be described by

where |$S_{ij}=\frac{1}{2}(\frac{\partial v_i}{\partial x_j}+\frac{\partial v_j}{\partial x_i})$| is the strain-rate tensor, which has a value of approximately |$|v_\mathrm{im,h}-v_\mathrm{im,c}|/\Delta$| across the mixing layer. The numerical kinetic viscosity |$\nu$| is estimated by |$\nu =\alpha c_\mathrm{s,c}\Delta$|⁠. Equating the viscous dissipation term to the cooling term, we have

Using |$n^2\Lambda =P/[(\gamma -1)t_{\rm cool}]$| and considering the cold phase yields

where |$\mathcal {M}_\mathrm{im,h}$| and |$\mathcal {M}_\mathrm{im,c}$| are the Mach numbers of hot and cold gas with respect to |$c_\mathrm{s,c}$|⁠, respectively. As shown by the blue line in Fig. A4, equation (A29) with |$\alpha =1.5$| is an excellent fit to simulation results, which use ramses with a MUSCL scheme and HLLC Riemann solver. This is true whether we vary |$\Delta$| by keeping the sheet thickness fixed and varying the number of cells across the sheet (points with a black boundary), or whether we keep the number of cells across the sheet fixed and vary the sheet thickness (points without a black boundary).

A2.3 The impact of turbulence

In more realistic cases, turbulence is ubiquitous and shapes the properties of mixing layers along with radiative cooling. In this case, the width of the mixing layer is comparable to the largest eddy, which has a comparable size of the width of the sheet. According to Gronke & Oh (2020a), the entrainment velocity of the hot phase is

when cooling dominates over turbulent mixing, where |$\beta =0.2-0.5$| (Gronke & Oh 2020a, 2023). Combining this with equation (A22), the velocity difference is

In order to develop turbulence during the implosion, we add interface perturbations to the surfaces of sheets and density perturbations inside the sheets, as described in Section 2. As shown by the orange circular dots and orange line in Fig. A4, we find that the velocity differences in simulations are well fitted by equation (A31) with |$\beta =0.4$|⁠. This result is independent of resolution, measured by the number cells across the sheet thickness, as shown by the orange points with a black boundary, where we vary the number of cells by factors of 2, 4, and 8 while keeping the sheet size fixed.

Until now, we have considered a case with an initial pressure discontinuity between the cold and hot gas. We now consider a case as in the main text where the sheet is initially in pressure equilibrium with the hot background, but its temperature is a factor of |$\eta$| higher than the temperature floor. Assuming |$t_\mathrm{cool}\ll t_\mathrm{sc}$|⁠, which is true in our case, the sheet rapidly cools isochorically and induces the pressure contrast. As shown by the orange triangles in Fig. A4, the properties of the mixing layer in this case are similar to the case with an initial pressure contrast.

If the movement of mixing layer with respect to the implosion shock can be neglected, then the total implosion velocities for the cold and hot gas are given by |$v_\mathrm{im,c}=v_\mathrm{im,iso}+v_\mathrm{ent,c}$| and |$v_\mathrm{im,h}=v_\mathrm{im,iso}+v_\mathrm{ent,h}$|⁠, respectively. Furthermore, the thermal pressure of imploding gas decreases near the mixing layer, which drives the inflow into the mixing layer. In other words, the pressure and density are modified due to the entrainment flow

where we have assumed the density varies adiabatically because the cold gas is at the cooling floor and the hot gas cannot cool (recall that it is primarily intermediate gas that cools in the mixing layer). We have thus obtained the properties of the implosion, accounting for the effects of the radiatively cooling mixing layer. Similar modifications are applicable to the explosion phase as well, where the entrainment accelerates the cold gas and decelerates the hot gas.

We compared our model predictions with simulations in Fig. A5. Here, we vary |$\chi _\mathrm{f}$| while fixing |$\eta =10$| and |$\, {r_{\rm s,i}}=3\,$| kpc. Circles represent simulations with an initial pressure contrast, while triangles represent simulations where the pressure contrast was induced by rapid cooling. Note that the triangles at |${\chi _{\rm f}}=100$| correspond to the sheet case shown in Fig. 2. For all values of |${\chi _{\rm f}}$|⁠, the model is a good fit to the density and temperature values at all stages – implosion, peak collision, and explosion – and a good fit to the pressure and radial veloity during the implosion. However, we overpredicted the peak collision pressure and the explosion velocity, especially at large |${\chi _{\rm f}}$|⁠. We suspected that this discrepancy arises from the fragmentation of the sheet during the implosion, as illustrated in Fig. A6. This fragmentation seems to be due to the surface perturbations and a combination of RMI and vorticity induced at the curved surface of the perturbations together with strong cooling. Clumps begin to form and detach from the sheet very early on in the evolution, before the outer cold-gas surface has gained much velocity. At high overdensity, these detached clumps are unable to catch up to the implosion front since the time-scale for clumps to become entrained in the flow increases linearly with |$\chi$| (Section 3.2.1). These fragmented clumps increase the cold-gas surface area, thereby driving stronger inflows of hot background gas into the mixing layers. Furthermore, because of this fragmentation the imploding gas that reaches the centre is not fully collisional, with small-scale clumps ‘missing’ each other. This leads to a lower collisional density and pressure than predicted, ultimately resulting in a reduced explosion velocity.

Following the explosion, the cold gas escapes outwards at a velocity of |$\sim c_\mathrm{s,c}$|⁠, while a mild shock propagates outwards into the hot gas. In the adiabatic case (Fig. A1, red lines), the velocity of the post-shock hot gas is also |$\sim c_\mathrm{s,c}$|⁠. However, radiative cooling at the surface of the central cold-gas sheet induces an entrainment flow through a turbulent mixing layer, as described in Section 3.2.1 of the main text. This causes the hot gas to flow towards the centre in an entrainment flow, while the cold clumps escape outwards.

APPENDIX B: CONVERGENCE TESTS

B1 The peak pressure and explosion velocity of sheets

In Fig. 2, we found that the peak pressure and the explosion velocity in sheets are both much smaller than in streams and spheres. Since the density and pressure in collapsing isothermal clouds increases fastest for spheres and slowest for sheets, we wished to verify that we would not see a more dramatic increase in the central pressure (and therefore in explosion velocity) in sheets if the resolution were increased allowing further collapse. To this end, we performed sheet simulations with the same parameters as in Fig. 2, with three different resolutions within |$r\lt 3\, r_\mathrm{s,i}$|⁠. The results are shown in Fig. B1. We found that the peak pressure and explosion velocity do not change with resolution.

B2 |$f_\mathrm{A}$| and |$r_\mathrm{fA}$|

In Section 6.1, we derived a coagulation criterion that depends on |$r_{\rm fA}$|⁠, the radius where |$f_{\rm A}$| began to decrease. In order to verify that this is not sensitive to the artificial disruption of cold gas induced by our statically refined mesh (Section 2, Fig. 1), we performed two tests with a larger high-resolution region, such that the resolution first decreases at |$r=6\, r_\mathrm{s,i}$| rather than |$3\, r_{\rm s,i}$|⁠. In other words, we increased the resolution by one refinement level (corresponding to cell sizes half as large) for all cells outside the initial high-resolution region, |$r\gt 3\, \, {r_{\rm s,i}}$|⁠, while keeping the same high-resolution interior to this radius. As a result, the first drop in resolution now occurs at |$6\, \, {r_{\rm s,i}}$| instead of |$3\, \, {r_{\rm s,i}}$|⁠, with subsequent drops occurring at |$[9,12]\, {r_{\rm s,i}}$|⁠. We did this for both streams and spheres with |$\chi _\mathrm{f}=1000$|⁠, as shown in Fig. B2. We found that the size of the high-resolution region does not affect |$f_\mathrm{A}$| significantly.

B3 The coagulation criterion

Our definition of fragmentation or coagulation rests primarily on the evolution of the number of clumps. To verify that this is not sensitive to resolution, we chose four simulations that are close to the borderline and lowered their resolutions everywhere by factors of 2 and 4, as shown in Fig. B3. While the number of clumps obviously depends on resolution, the distinction between fragmentation and coagulation, namely whether |$N_{\rm c}$| decreases to order unity or remains large until the end of the simulation, does not.

 

APPENDIX C: VOLUME-WEIGHTED AVERAGE OF DENSITY PROJECTION MAPS

In Figs C1–C4, we presented maps of the volume-weighted average density in sheets, streams, and spheres, with different |${\chi _{\rm f}}$| and metallicity values. These are meant to complement, respectively, Figs 5–7 and 11 from the main text, which show the maximal density along the line of sight. As described in the text, the maximal density better highlights thermal fragmentation and small clumps, while the average density better highlights coagulation and large clumps. This better shows how coagulation is suppressed along the stream axis and in the plane of the sheet, and also how in general radial coagulation is stronger in sheets compared to streams compared to spheres.

APPENDIX D: LOGNORMAL FITS TO THE CLUMP MASS DISTRIBUTIONS

In Section 5 and Fig. 15, we discussed the distribution of clump sizes, suggesting these are well described by a power law, |$N(\gt m)\propto m^{-1}$|⁠. However, many power laws are actually lognormal distributions in disguise (Mitzenmacher 2004), and pure fragmentation, in particular, is expected to produce a lognormal due to the central limit theorem (Kolmogorov 1941). We compared these distributions using the python package powerlaw (Alstott et al. 2014), which fits both power laws and lognormals while estimating the p-values of the significance via log-likelihood ratios. When we used all clumps down to the minimal clump mass, a lognormal distribution is significantly favoured over a power, law with a p-value less than |$10^{-20}$|⁠. However, when allowing the minimal clump mass to vary as part of the fitting procedure, the p-value can rise to as high as 0.3 when fitting over the range |$m_{\rm cl}/m_0\gt 10^{-4}$|⁠. This suggests that while the lognormal distribution is still favoured, a power law is also a good description of the data over this mass range. In either case, the mode of the lognormal distribution is sensitive to resolution even when |$\ell _\mathrm{shatter}$| is resolved, consistent with our main conclusion that |$\ell _\mathrm{shatter}$| does not set a characteristic value for clump size in our simulations.

Fig. D1 shows the lognormal fits to the case of a low-Z stream with |$\, {r_{\rm s,i}}=3\,$| kpc, |$\eta =10$|⁠, and |${\chi _{\rm f}}=100$| at different resolutions. Fig. D2 shows the PDF of all low-Z cases, which can be fit by lognormal distributions similar to Fig. D1. Note that the peak of the distribution tracks the numerical resolution. Fig. D3 shows the probability density function (PDF) of column densities of all low-metallicity streams. The peaks of the column densities are determined by the resolution rather than by the column density at |$\ell _\mathrm{shatter}$|⁠.