Shattering of Cosmic Sheets due to Thermal Instabilities: a Formation Channel for Metal-Free Lyman Limit Systems

TL;DR

This paper uses high-resolution cosmological simulations to reveal how thermal instabilities in large-scale structure lead to the formation of metal-free Lyman limit systems outside galaxies, providing a new formation channel.

Contribution

It introduces a novel simulation approach showing the formation of pristine, cold gas clouds via sheet shattering due to thermal instabilities in the intergalactic medium.

Findings

Cold clouds have N_{HI}>10^{17.2}cm^{-2} and are detectable as Lyman limit systems.

Clouds are metal-free and located outside galaxy halos.

Sheet shattering occurs from merger-induced shocks leading to thermal instabilities.

Abstract

We present a new cosmological zoom-in simulation, where the zoom region consists of two halos with virial mass M_v~5x10^{12}M_{sun} and a ~Mpc long cosmic filament connecting them at z~2. Using this simulation, we study the evolution of the intergalactic medium in between these two halos at unprecedented resolution. At 5>z>3, the two halos are found to lie in a large intergalactic sheet, or "pancake", consisting of multiple co-planar dense filaments along which nearly all halos with M_v>10^9M_{sun} are located. This sheet collapses at z~5 from the merger of two smaller sheets. The strong shock generated by this merger leads to thermal instabilities in the post-shock region, and to a shattering of the sheet resulting in <~kpc scale clouds with temperatures of T>~2x10^4K and densities of n>~10^{-3}cm^{-3}, which are pressure confined in a hot medium with T~10^6K and n>~10^{-5}cm^{-3}. When the sheet is viewed face on, these cold clouds have neutral hydrogen column densities of N_{HI}>10^{17.2}cm^{-2}, making them detectable as Lyman limit systems, though they lie well outside the virial radius of any halo and even well outside the dense filaments. Their chemical composition is pristine, having zero metalicity, similar to several recently observed systems. Since these systems form far from any galaxies, these results are robust to galaxy formation physics, resulting purely from the collapse of large scale structure and radiative cooling, provided sufficient spatial resolution is available.

Full text

††slugcomment: Submitted to ApJL

1. Introduction

Only a small fraction of the Universe’s baryons and metals belong to galaxies (e.g. Tumlinson, Peeples & Werk, 2017; Wechsler & Tinker, 2018). The rest reside in the circumgalactic medium (CGM), the space outside galaxies but within their host dark matter halo, and the intergalactic medium (IGM), the space in between dark matter halos. Both of these baryonic reservoirs are intimately linked to galaxy evolution through cycles of gas accretion, star-formation, galactic outflows, and eventual re-accretion (e.g. Putman, Peek & Joung, 2012; McQuinn, 2016; Tumlinson, Peeples & Werk, 2017). Thus, the physical properties and chemical composition of the IGM and CGM offer valuable insight into processes related to galaxy formation and evolution.

In recent decades, the low density gas in the IGM and CGM has been probed using absorption line spectroscopy along lines of sight to distant QSOs or galaxies (e.g. Lynds, 1971; Hennawi et al., 2006; Steidel et al., 2010). Intervening gas clouds with low neutral hydrogen column densities, N_{\rm HI}\lower 2.15277pt\hbox{;\buildrel<\over{\sim};}10^{15}\,{\rm cm}^{-2}, are understood to reside in the IGM and comprise the Lyman-$\alpha$ Forest, hereafter LyAF. This gas is thought to trace fluctuations in the underlying dark matter distribution which are still in the linear regime, making diagnostics of the LyAF a powerful tool to constrain cosmology (see Rauch, 1998 and McQuinn, 2016 for reviews). Clouds with high column densities, $N_{\rm HI}>10^{17.2}\,{\rm cm}^{-2}$, are optically thick bluewards of the Lyman limit, $\lambda<$ 912Å, and are referred to as Lyman Limit Systems (LLSs). At redshifts 2\lower 2.15277pt\hbox{;\buildrel<\over{\sim};}z\lower 2.15277pt\hbox{;\buildrel<\over{\sim};}5, LLSs exhibit a broad distribution of metalicities. The bulk of the population has $Z\sim 10^{-2}\,Z_{\odot}$, while a handful of systems have $Z<10^{-3}\,Z_{\odot}$ (Fumagalli, O’Meara & Prochaska, 2016; Lehner et al., 2016; Robert et al., 2019).

LLSs, particularly those with $Z>10^{-3}\,Z_{\odot}$, are commonly thought to reside in the CGM rather than the IGM (Sargent, Steidel & Boksenberg, 1989; Fumagalli, O’Meara & Prochaska, 2016; Lehner et al., 2016). However, the recent discovery of several LLSs with $Z<10^{-3.4}\,Z_{\odot}$ as well as potentially pristine LLSs at $3<z<5$ has led some to question whether these may represent a separate population originating in the IGM (Fumagalli, O’Meara & Prochaska, 2011; Crighton, O’Meara & Murphy, 2016; Robert et al., 2019). A single PopIII supernova would pollute gas to higher metalicity values (Wise et al., 2012; Crighton, O’Meara & Murphy, 2016), and simulations of structure formation which include PopIII star formation suggest that such low metalicities exist only in isolated low-density patches of the IGM (e.g. Tornatore, Ferrara & Schneider, 2007; Wise et al., 2012). Alternatively, pristine LLSs may originate in cold accretion streams feeding massive galaxies from the IGM (Dekel et al., 2009; Fumagalli et al., 2011; Faucher-Giguère & Kereš, 2011; van de Voort et al., 2012). While cosmological simulations suggest that the typical metalicity in such streams is \lower 2.15277pt\hbox{;\buildrel>\over{\sim};}10^{-3}\,Z_{\odot} at $z<5$ (van de Voort & Schaye, 2012; Ceverino et al., 2015; Mandelker et al., 2018), lower metalicity clouds may still be present. However, it has also been suggested that the evolution of the number of LLSs per unit redshift at $z>3.5$ is inconsistent with a contribution from the CGM alone, indicating a growing contribution of LLSs in the IGM above this redshift (Fumagalli et al., 2013). All in all, the origin of extremely metal-poor Lyman limit systems in the IGM is not yet understood.

It is notoriously difficult to study the detailed properties of gas in the IGM and CGM with cosmological simulations. The resolution in most state-of-the-art simulations is adaptive in a quasi-Lagrangian sense, such that the effective mass resolution is fixed. Consequently, the spatial resolution becomes very poor in the low density CGM and even worse in the IGM (Nelson et al., 2016), orders of magnitude larger than the cooling length of $T\sim 10^{4}\,{\rm K}$ gas, $l_{\rm cool}=c_{\rm s}t_{\rm cool}\sim 100\,{\rm pc}~{}(n/10^{-3}\,{\rm cm}^{-3})^{-1}$, where $c_{\rm s}$ is the sound speed and $t_{\rm cool}$ is the cooling time (McCourt et al., 2018, hereafter M18; Sparre, Pfrommer & Vogelsberger, 2019). While several groups have recently introduced different methods to better resolve the CGM (Hummels et al., 2018; Corlies et al., 2018; Peeples et al., 2019; Suresh et al., 2019; van de Voort et al., 2019), we are unaware of similar attempts to better resolve the IGM.

In this letter, we introduce a new cosmological simulation where we zoom-in on a large region of the IGM in between two massive galaxies at $z\sim 2.3$, with a comoving separation of $\sim 3\,{\rm Mpc}/h$. This is the highest resolution simulation of such a large patch of the IGM to date. Using this simulation, we investigate the evolution of the IGM and show how thermal instabilities triggered by shocks during the collapse of large-scale structure can lead to the formation of pristine LLSs, far from any galaxies. The simulation is described in §2. In §3 we present our results, and we conclude in §4. Throughout, we assume a flat $\Lambda$CDM cosmology with $\Omega_{\rm m}=1-\Omega_{\Lambda}=0.3089$, $\Omega_{\rm b}=0.0486$, $h=0.6774$, $\sigma_{8}=0.8159$, and $n_{\rm s}=0.9667$ (Planck Collaboration et al., 2016).

2. Simulation Method

We perform simulations using the quasi-Lagrangian moving-mesh code AREPO (Springel, 2010). To select our target halos, we first consider the 200 most massive halos at $z\sim 2.3$ in the Illustris TNG100111http://www.tng-project.org magnetohydrodynamic cosmological simulation (Pillepich et al., 2018a; Nelson et al., 2018; Springel et al., 2018). These span a mass range of $M_{\rm v}\sim(0.7-27)\times 10^{12}{\rm M}_{\odot}/h$, where $M_{\rm v}$ is the virial mass defined using the Bryan & Norman (1998) spherical overdensity. We then select all pairwise combinations of them with a comoving distance in the range $(2.5-4.0)\,{\rm Mpc}/h$, finding 48 such halo pairs. Visual inspection revealed each such halo pair to be directly connected by a dark matter cosmic web filament, with comparable radius to the halo virial radii. One such pair was randomly chosen for resimulation, consisting of two halos with $M_{\rm v}\sim 3.4\times 10^{12}{\rm M}_{\odot}/h$ each, separated by a proper distance of $D\sim 0.8\,{\rm Mpc}/h$. At $z=0$, the two halos have masses $M_{\rm v}\sim(1.1-1.3)\times 10^{13}{\rm M}_{\odot}/h$ and are $\sim 1.8\,{\rm Mpc}/h$ apart, so their comoving distance has decreased by \lower 2.15277pt\hbox{;\buildrel<\over{\sim};}30\%.

We define $R_{\rm ref}=1.5\times R_{\rm v,max}\sim 163\,{\rm kpc}/h$, with $R_{\rm v,max}$ the larger of the two virial radii at $z=2.3$. The zoom-in region is the union of a cylinder with radius $R_{\rm ref}$ and length $D$ extending between the two halo centers, and two spheres of radius $R_{\rm ref}$ centred on either halo. We trace all dark matter particles within this volume back to the initial conditions of the simulation, at $z=127$, refine the corresponding Lagrangian region to higher resolution, and rerun the simulation to $z=2$, when the region of interest by construction becomes contaminated by low resolution material from outside the refinement region. The simulations were performed with the same physics model used in the TNG100 simulation, described in detail in Weinberger et al. (2017) and Pillepich et al. (2018b). We briefly summarize below the implementation of the ionizing radiation field and of cooling, which are most relevant to our current work.

We follow the production and evolution of nine elements (H, He, C, N, O, Ne, Mg, Si, and Fe). These are produced in supernovae Type Ia and Type II and in AGB stars according to tabulated mass and metal yields. Metal line cooling is included using pre-calculated values as a function of density, temperature, metalicity and redshift, with corrections for self-shielding (Wiersma, Schaye & Smith, 2009). The metal enriched gas radiatively cools in the presence of a spatially uniform but redshift dependent ionizing UV background (UVB; Faucher-Giguère et al., 2009), which is instantaneously switched on at $z=6$. To minimize any potential influence of this instantaneous switching on of the UVB, we limit our current analysis to $z\leq 5$. Cooling is further modulated by the radiation field of nearby active galactic nuclei (AGN) by superimposing the UVB with the AGN radiation field (Vogelsberger et al., 2013).

We performed five simulations with different resolution within the refinement region. A detailed convergence study will be presented in an upcoming paper (Mandelker et al., in prep.). In the current letter we focus on our highest resolution simulation, which has a dark matter particle mass of $m_{\rm dm}=8.2\times 10^{4}\,{\rm M}_{\odot}$ and a Plummer-equivalent gravitational softening of $\epsilon_{\rm dm}=250\,{\rm pc}$ comoving. Gas cells are refined such that their mass is within a factor of 2 of $m_{\rm gas}=1.5\times 10^{4}{\rm M}_{\odot}$, and have a minimal gravitational softening $\epsilon_{\rm gas}=0.5\epsilon_{\rm dm}$. We compare this to a simulation with comparable resolution to TNG100, having $m_{\rm dm}=5.3\times 10^{6}{\rm M}_{\odot}$, $m_{\rm gas}=1.0\times 10^{6}{\rm M}_{\odot}$, and $\epsilon_{\rm dm}=2\epsilon_{\rm gas}=1000\,{\rm pc}$ comoving.

3. Results

In LABEL:NH we show the evolution of the large scale structure surrounding our system, at $z\sim 5$, $4$, and $3$. The left and center columns show the total hydrogen column density, $N_{\rm H}$, in two orthogonal projections, with the intergalactic sheet containing the two halos shown edge-on and face-on, respectively. At $z>5$ the system actually consists of two sheets initially inclined to one another, marked by dashed lines in panel A. These merge at $z\sim 5$, with only a single sheet visible in panels D and G. The sheet contains several prominent coplanar filaments, with end-points at either of the two main halos and along which lie nearly all halos with $M_{\rm v}>10^{9}M_{\odot}$. Most of these filaments merge at $z<3$, leaving behind the single giant filament selected at $z=2.3$. The beginning of this merger is visible in panel H. We note that the configuration of our system at $z\sim 3$ is remarkably similar to a system recently observed at $z\sim 3.2$ with MUSE (Lusso et al., 2019).

Following the sheet collision at $z\sim 5$, several regions between the filaments in the post-merger sheet develop a granular morphology. As the merging sheets were initially inclined, the collision and resulting granular structure propagate from left-to-right in LABEL:NH, as can be seen by comparing panels B and E. In panels B, E, and H we highlight such granular regions, denoted C, F, and I, selected to contain no halos with $M_{\rm v}>10^{9}{\rm M}_{\odot}$ and to not intersect any filaments. In the right-hand column we show the neutral hydrogen column density, $N_{\rm HI}$, within these regions. The granular structure is even more prominent in $N_{\rm HI}$, with many regions being classified as LLSs, $N_{\rm HI}>10^{17.2}\,{\rm cm}^{-2}$. Importantly, these regions do not coincide with the locations of dark matter halos, or with fluctuations in the dark matter distribution which is smooth in these regions. We examined several similar regions within the sheet in each snapshot, also marked in panels B, E, and H, and found the gas properties to be very similar in all regions within the post-merger sheet.

As we argue below, this granular morphology seems to be triggered by nonlinear thermal instabilities within the post-shock sheet. Based on the model presented in M18, nonlinear thermal instabilities in a rapidly cooling medium cause the medium to “shatter”, forming dense cloudlets with T\lower 2.15277pt\hbox{;\buildrel>\over{\sim};}10^{4}\,{\rm K} in pressure equilibrium with a more tenuous, hot background. The size of these cloudlets is set by the local cooling length, $l_{\rm cool}=c_{\rm s}t_{\rm cool}$. This procedure is hierarchical, in the sense that as the gas cools, $l_{\rm cool}$ decreases, causing existing cloudlets to shatter into even smaller cloudlets. We note that there are several differences between our system and the idealized study of M18. Firstly, the gas in our system is photo-heated by the UV background, while M18 considered a purely cooling system in collisional-ionization equilibrium. Second, in M18 the external pressure is set by the thermal pressure in the hot background, where the cooling time is assumed to be much longer than the shattering timescale in the cooling medium, while in our case it is set by the ram pressure of the infalling material, as discussed below. Finally, our system is in 3D while those studied in M18 were 2D. Nevertheless, as we argue below our results appear consistent with the M18 shattering model.

In LABEL:Temp we show the projected, density-weighted gas temperature in the same frame as panels D and E from LABEL:NH. A planar accretion shock around the sheet, triggered by the earlier collision, is clearly visible in the edge-on view, as are spherical accretion shocks around the two main halos. In the face-on view, the filaments appear cold, with $T\sim 2\times 10^{4}\,{\rm K}$, while the regions between filaments exhibit a multiphase structure, with hot and cold regions coexisting in a granular structure similar to that seen in the column density (LABEL:NH). In the right-hand panel of LABEL:Temp, we show the projected metalicity in the face-on view of the sheet at $z\sim 4$. While the filaments are enriched to Z\lower 2.15277pt\hbox{;\buildrel>\over{\sim};}10^{-2}{\rm Z}_{\rm\odot}, consistent with previous studies (van de Voort & Schaye, 2012; Ceverino et al., 2015; Mandelker et al., 2018), the regions in between the filaments retain near-pristine compositions with $Z<10^{-3}{\rm Z}_{\rm\odot}$, due to their large distance from any star-forming galaxies.

LABEL:phase_diagrams shows the distribution of gas within region F from LABEL:NH in density-temperature space, weighted by mass (left) and by volume (middle). The post-shock gas has roughly constant pressure, $P_{\rm th}\sim 50~{}\,{\rm K}~{}\,{\rm cm}^{-3}$, roughly the ram pressure of the infalling gas which has a density of $\sim 4\times 10^{-29}{\rm g~{}\,{\rm cm}^{-3}}$ and a velocity of $\sim 200\,{\rm km}\,{\rm s}^{-1}$. Two phases in approximate pressure equilibrium are apparent, with most of the mass at $(n_{\rm H},T)\sim(10^{-3}\,{\rm cm}^{-3},10^{4.5}\,{\rm K})$, and most of the volume at $\sim(10^{-4.5}\,{\rm cm}^{-3},10^{6}\,{\rm K})$. In the right panel we show the cooling length, $l_{\rm cool}=c_{\rm s}t_{\rm cool}$, for gas undergoing net cooling. The post-shock gas, which has a temperature of $T\sim 2\times 10^{5}\,{\rm K}$, has a cooling length of $l_{\rm cool}\sim 10\,{\rm kpc}$, while the cold phase along the same isobar has $l_{\rm cool}\sim 1\,{\rm kpc}$.

In LABEL:slice we show the hydrogen density in a slice through the sheet-midplane in region F from LABEL:NH. On the left we show our highest resolution simulation, while on the right we show a simulation with similar resolution to TNG100. In the former, the typical (minimal) cell size within this region is $\Delta\sim 0.8~{}(0.3)\,{\rm kpc}$, significantly smaller than the post-shock cooling length and comparable to the cooling length in the dense phase. We are thus able to resolve the onset of shattering into dense $\sim\,{\rm kpc}$ scale cloudlets (M18), resulting in the large neutral column densities seen in LABEL:NH. However, since the minimal cooling length is only marginally resolved, the end result is not converged and the actual cloudlets are expected to be smaller. On the other hand, in the simulation with TNG100 resolution, the typical (minimal) cell size within this region is $\Delta\sim 4.0~{}(2.5)\,{\rm kpc}$. The initial phases of the shattering are thus unresolved, and no dense cloudlets are formed. We note that the thermal Jeans length in the cold phase is $L_{\rm J}=[9c_{\rm s}^{2}/(4\pi G\rho)]^{1/2}\sim 30\,{\rm kpc}$, significantly larger than the cooling length, the cloud sizes, and the typical cell size. This implies that the clouds are not the result of gravitational instability in the sheet, and supports our hypothesis that they result from thermal instabilities and shattering.

In LABEL:convergence we show the covering fraction of neutral hydrogen as a function of $N_{\rm HI}$, in the two simulations shown in LABEL:slice. We show results at $z\sim 5$, $4$, and $3$, corresponding respectively to regions C, F, and I in LABEL:NH. In order to focus on gas with primordial composition that condenses due to thermal instabilities rather than fluctuations in the underlying dark matter distribution as is often assumed in studies of the LyAF, when evaluating $N_{\rm HI}$ we ignore all cells with metalicity $Z>10^{-4}\,Z_{\odot}$ and with a dark matter density greater than 5 times the Universal mean at the relevant redshift. At all redshifts, the covering fraction of LLSs is significantly larger in our high resolution simulation than in the simulation with TNG100 resolution. As discussed above, this is because the latter does not resolve the initial shattering of the post-shock medium within the sheet. At z\lower 2.15277pt\hbox{;\buildrel<\over{\sim};}4 we find covering fractions of order $\sim 1\%$ for pristine LLSs. At $z\sim 5$ the covering fraction is $\sim 35\%$. We find comparable covering fractions in other sheet regions that do not intersect any filaments or massive halos. Our results are qualitatively similar using metalicity thresholds as large as $0.1\,Z_{\odot}$ and dark matter overdensity thresholds in the range $2.5-10$.

A detailed convergence study of the covering fraction of neutral hydrogen in the IGM in our five simulations with varying resolutions, accounting for different viewing angles through the sheet, will be presented in an upcoming study (Mandelker et al. in prep.). Here, we wish to highlight in LABEL:convergence the fact that metal-free LLSs in the IGM occur naturally in our simulations with sufficient resolution, with non-negligible covering fractions. Furthermore, it is interesting to note the decline in the covering fraction of dense clouds with redshift, by a factor of \lower 2.15277pt\hbox{;\buildrel>\over{\sim};}70 from $z\sim 5$ when the clouds are formed following the sheet collision, to $z\sim 3$. This decline may be caused by cosmic expansion, which causes the pressure in the sheet to decline by a factor of \lower 2.15277pt\hbox{;\buildrel>\over{\sim};}5 from $z=(5-3)$, in rough agreement with the naive scaling of $P\propto(1+z)^{5}$ in the IGM. As the cold clouds all have approximately the same temperature, the typical cloud density declines by a similar factor, causing the neutral fraction to decline by an even larger factor. Alternatively, cold clouds moving rapidly through a hot medium are expected to be disrupted on a cloud-crushing timescale, $t_{\rm cc}=2R_{\rm cl}v_{\rm cl}^{-1}\delta^{1/2}$, where $R_{\rm cl}$ is the cloud radius, $v_{\rm cl}$ its velocity, and $\delta$ is the density ratio between the cloud and the background (e.g. Agertz et al., 2007). In our case, $R_{\rm cl}\sim 2\,{\rm kpc}$, $v_{\rm cl}\sim 50\,{\rm km}\,{\rm s}^{-1}$ comparable to the sound speed in the post-shock medium, and $\delta\sim 30$ (LABEL:phase_diagrams). This yields $t_{\rm cc}\sim 0.2\,{\rm Gyr}$, while the time between $z=(5-3)$ is $\sim 1\,{\rm Gyr}$. This would imply that clouds are continuously created during this period, presumably due to turbulence in the sheet.

4. Discussion and Conclusions

We have presented a new cosmological simulation which zooms in on two massive, $\sim 5\times 10^{12}\,{\rm M}_{\odot}$ halos connected by a Mpc-scale cosmic filament at $z\sim 2$. This large zoom-in region enables us to resolve the IGM, far from any massive halo, at unprecedented resolution. The simulation reveals the growth of the large scale structure and the cosmic web around our system. Two inclined cosmic sheets collide at $z\sim 5$ to produce a single massive sheet containing numerous dense filaments along which nearly all halos with $M_{\rm v}>10^{9}{\rm M}_{\odot}$ lie. Following a major merger of one of our two main halos at z\lower 2.15277pt\hbox{;\buildrel<\over{\sim};}3, the filaments merge leaving a single massive filament.

The sheet collision at $z\sim 5$ triggers a strong shock within the resulting sheet, heating the gas to $T\sim 2\times 10^{5}\,{\rm K}$ with its pressure fixed by the ram pressure of infalling gas at $P\sim 50\,{\rm K}\,{\rm cm}^{-3}$. Due to thermal instability, this gas then separates into a volume-filling and a mass-dominating phase with $(n_{\rm H},T)\sim(10^{-4.5}\,{\rm cm}^{-3},10^{6}\,{\rm K})$ and $(10^{-3}\,{\rm cm}^{-3},10^{4.5}\,{\rm K})$ respectively. The cold phase is produced by shattering of the sheet into cloudlets with sizes comparable to the local cooling length, $l_{\rm cool}=c_{\rm s}t_{\rm cool}$ (McCourt et al., 2018). In the post-shock gas, $l_{\rm cool}\sim 10\,{\rm kpc}$ and the typical cell size in our high resolution simulation is $800\,{\rm pc}$. We thus resolve the shattering of the sheet into kpc-scale fragments. While this is comparable to $l_{\rm cool}$ in the cold phase, we caution that the cloud sizes are likely influenced by our numerical resolution. These dense cloudlets result in high-column densities of neutral hydrogen in regions of the sheet which do not intersect any filaments or contain any halos with $M_{\rm v}>10^{9}{\rm M}_{\odot}$. In particular, the covering fraction of pristine LLSs, with $N_{\rm HI}>10^{17.2}\,{\rm cm}^{-2}$, $Z<10^{-4}\,Z_{\odot}$, and dark matter density less than 5 times the universal mean, is $\sim 1\%$ at $z\sim 3$ and 4 and $\sim 35\%$ at $z\sim 5$ when the sheet is viewed close to face-on. Whether individual cloudlets can be observationally disentangled from each other or from the filaments if the sheet were viewed close to edge-on, and whether such viewing angles might yield metal-free Damped Lyman-$\alpha$ Absorbers (DLAs, $N_{\rm HI}>10^{20.3}\,{\rm cm}^{-2}$) will be presented in an upcoming paper (Mandelker et al. in prep.). Preliminary results suggest that the latter seems unlikely.

The large distance of these clouds from any massive galaxy implies that this result is likely robust to the adopted galaxy formation sub-grid physics. The production of metal-free LLSs, as recently observed by Fumagalli, O’Meara & Prochaska (2011); Crighton, O’Meara & Murphy (2016); Robert et al. (2019), thus seems to occur naturally in the IGM due to thermal instabilities induced by the growth of large-scale structure. This supports speculations that a growing fraction of LLSs at $z>3$ can be found in the IGM rather than the CGM around galaxies.

In simulations with resolution comparable to the Illustris TNG100 simulation, the typical cell size in the sheet regions is $\sim 4\,{\rm kpc}$. Therefore, the cooling length in the post-shock gas is only marginally resolved, and the shattering process does not take place. As a result, the covering fraction of pristine LLSs is greatly reduced at $z\sim 5$ and is \lower 2.15277pt\hbox{;\buildrel<\over{\sim};}10^{-4} at z\lower 2.15277pt\hbox{;\buildrel<\over{\sim};}4. This highlights the importance of achieving high-resolution in the IGM, even well outside the CGM of any galaxy.

Another potential application of our results is with regard to the LyAF, and its relation to the underlying dark matter density distribution. The LyAF can be used to probe the mildly nonlinear matter power spectrum, down to scales of $<1$ comoving Mpc. At $z\sim 4-6$, where the physical scales probed are several tens of kpc, the power at these small scales has been used to rule out various warm dark matter (WDM) models (e.g. Viel et al., 2013, and references therein), as a too light WDM particle would suppress power on these small scales compared to observations. These studies are often calibrated against cosmological simulations which have much lower resolution in the IGM than the simulations in this work. However, thermal instabilities such as those identified in this work may lead to excess power in the LyAF on scales of tens of kpc, corresponding to the cooling length of post-shock gas in cosmic sheets (LABEL:phase_diagrams), which is not associated with the underlying dark matter distribution. This may influence the constraints on WDM models from simulations which do not resolve the shattering. In addition to WDM constraints, our results may influence constraints on the temperature and the optical depth of the IGM at high redshift, which are also based on analysis of the LyAF (e.g. Viel et al., 2013; Lidz & Malloy, 2014; Eilers, Davies & Hennawi, 2018).

While the results presented in this paper are likely robust to the galaxy formation sub-grid physics, we caution against drawing too broad conclusions from the single system simulated in this study. In particular, it is currently unclear whether the shattering process discussed here requires a violent collision between sheets, or whether smooth accretion would have the same effect. It is also currently unknown how frequent such collisions between sheets are. Therefore, we cannot confidently extrapolate from the results presented here to the actual number density of metal-free LLSs in the IGM produced by shattering of cosmic sheets. Future simulations which employ similar methods of enhancing the resolution in the IGM will help to shed light on this problem.

Acknowledgments

We thank the anonymous referee for constructive comments which helped improve the clarity of this manuscript. We thank Joe Hennawi, Siang Peng Oh, Drummond Fielding, Michele Fumagalli, and Ruediger Pakmor for helpful discussions. We acknowledge support from the Klauss Tschira Foundation through the HITS Yale Program in Astropysics (HYPA). F.C.v.d.B received additional support from the National Aeronautics and Space Administration through Grant No. 17-ATP17-0028 issued as part of the Astrophysics Theory Program.