Cosmic Reionization on Computers: Statistics, Physical Properties, and Environments of Lyman Limit Systems at z ∼ 6

Lyman limit systems (LLSs) are dense hydrogen clouds with high enough H i column densities to absorb Lyman continuum photons emitted from distant quasars. Their high column densities imply an origin in dense environments; however, the statistics and distribution of LLSs at high redshifts still remain uncertain. In this paper, we use self-consistent radiative transfer cosmological simulations from the Cosmic Reionization on Computers (CROC) project to study the physical properties of LLSs at the tail end of cosmic reionization at z ∼ 6. We generate 3000 synthetic quasar sight lines to obtain a large number of LLS samples in the simulations. In addition, with the high physical fidelity and resolution of CROC, we are able to quantify the association between these LLS samples and nearby galaxies. Our results show that the fraction of LLSs spatially associated with nearby galaxies increases with H i column density. Moreover, we find that LLSs that are not near any galaxy typically reside in filamentary structures connecting neighboring galaxies in the intergalactic medium (IGM). This quantification of the distribution and association of LLSs to large-scale structure informs our understanding of the IGM–galaxy connection during the “Epoch of Reionization,” and provides a theoretical basis for interpreting future observations.


INTRODUCTION
As the first stars and galaxies form, photons emitted from these luminous sources travel through the Universe.These photons encounter hydrogen atoms residing in a variety of environments, from the low-density, diffuse intergalactic Medium (IGM) to denser filaments of cosmic large-scale structures to even denser gas in the interstellar Medium (ISM) of galaxies.The first ionizing photons rapidly ionize the neutral hydrogen in the diffuse IGM and turn it transparent to ionizing radiation (sometimes called Lyman Continuum, or LyC).On the other hand, the denser regions remain progressively more neutral and more opaque to LyC.In ob-Corresponding author: Jiawen Fan johnfan@umich.eduserved absorption spectra of quasars and galaxies, the higher density and relatively more neutral regions manifest as Lyman Limit Systems (LLSs) (Sargent et al. 1989;Lanzetta 1991).LLSs are absorbers with the value of the hydrogen column density along the line-of-sight N HI > 1.6 × 10 17 cm −2 (which corresponds to the unit optical depth at the Lyman limit) and less than about 2 × 10 20 cm −2 (Crighton et al. 2019).The higher column density absorbers N HI > 2 × 10 20 cm −2 are traditionally called "Damped Lyα Systems", or DLAs.
Over the last decade, several groups have begun to explore how different physical processes and environmental factors affect the properties and distribution of LLSs.For example, some studies examining the effects of radiation from stars on the distribution of LLSs at z = 3 found minimal dependence on stellar feedback implementation (Yajima et al. 2012;Altay et al. 2013).Other studies focused on understanding LLS properties in the context of their association with galaxies and streams at these lower redshifts (Fumagalli et al. 2011;Rahmati & Schaye 2014;Erkal 2015).However, studies of LLSs at higher redshifts are sparse in the literature (Songaila & Cowie 2010;Crighton et al. 2019), likely due to the previous limitations in available observational data for comparison.In particular, observational detection of LLSs at the highest redshifts requires more quasar spectra, particularly with complete spectral coverage of both the Lyα and the Lyman break (Shukla et al. 2016).
It is important to understand LLSs at high z.It is well established that LLSs contribute more to the ionizing photon mean free path (MFP) than the Lyα forest (i.e.systems with N HI ≲ 1.6 × 10 17 cm −2 ; Miralda-Escudé 2003), and LLSs are much more frequently seen than the denser and rarer DLAs.Therefore, LLSs are expected to play a particularly important role during cosmic reionization.As the ionized bubbles grow and merge, the embedded LLSs restrict the propagation of ionizing photons inside the ionized bubbles.LLSs thereby limit the MFP of the LyC photons, eventually causing local reionization along that line of sight to stop.An understanding of the nature of LLSs requires a combination of improved observations and simulation models.
The recent launch of JWST has opened a new window to study galaxy formation during the Epoch of Reionization.JWST near-infrared and mid-infrared bands enable observers to study the earliest galaxies and to probe galaxy formation to fainter magnitudes (Endsley et al. 2023;Saxena et al. 2023;Jones et al. 2023;Atek et al. 2023;Yang et al. 2023;Wang et al. 2023;Treu et al. 2022;Prieto-Lyon et al. 2023;Finkelstein et al. 2023;Jung et al. 2023).This new accessibility of galaxy formation data during the Epoch of Reionization is a major advancement compared with the historic paucity of high redshift LLS data.
Theoretical models, particularly cosmological simulations, are necessary complements to the deluge of observational data.Fortunately, the rapid progress in developing fast and accurate methods for modeling radiative transfer in cosmological simulations allows us to study the LLSs at high redshifts at an unprecedented level of detail.
Realistic models are important to our understanding of high-density absorbers (e.g.Erkal 2015).With galaxy formation physics, variations in the balance between heating and cooling processes can lead to a range of effects on the distribution of absorbers.For example, more efficient cooling in gas clumps can lead to more compact dense gas distributions that background quasars are less likely to probe; these objects will fill the high-density tail of the true distribution of absorbers in the Universe but will suffer from an observational bias.On the other hand, energy feedback can expel high-density gas to larger radii, creating sufficiently dense gas for star formation in a "positive feedback" scenario (Ishibashi et al. 2013;Luisi et al. 2021).This kind of mechanism might potentially boost the number of intermediate-column-density absorbers at larger distances from galaxies.However, some studies on lowerredshift LLSs suggest LLS robustness to a realistic range of physical properties; star-forming gas in a galaxy will adjust such that outflows driven by feedback balance inflows due to accretion (Altay et al. 2013).
The structure of this paper is as follows.We describe the Cosmic Reionization on Computers (CROC) simulations and our methods to identify and quantify LLSs in Section 2, present the statistics and properties of LLSs in Section 3, and summarize our results in Section 4.

"Cosmic Reionization on Computers": CROC Simulation
CROC simulations use the Adaptive Refinement Tree code (Kravtsov et al. 1997(Kravtsov et al. , 2002;;Rudd et al. 2008).To account for the physics necessary to self-consistently model cosmic reionization, CROC simulations include gravity, gas dynamics, fully coupled (to gas dynamics) radiative transfer, atomic cooling and heating processes, molecular hydrogen formation, star formation, and stellar feedback.Full details of the simulation are described in the CROC method paper (Gnedin 2014).
For our analysis, we use three simulation runs with box length L box = 40 h −1 cMpc.The three runs are performed with exactly the same code and only differ in their random realizations of cosmological initial conditions.The simulations have a maximum spatial resolution of 100 pc in proper units.We generate 1000 lines of sight in each of the simulation boxes.These lines of sight start at random locations uniformly sampling the simulation volume, and they are oriented along random directions that sample the unit sphere uniformly.Each line of sight is sampled nonuniformly, fully tracing the underlying high resolution of the adaptively refined grid.We choose 100 h −1 cMpc for the length of the line of sight since it is a reasonable distance to use for the simulation boxes with L box = 40 h −1 cMpc, while providing a sufficient total length to allows us to find even rare absorbers.

Identifying Lyman limit systems
The observational signature of an LLS is a break at the Lyman limit (a drop in transmitted flux at 912 Å rest frame).In numerical simulations, we can simply look for regions with high enough neutral fraction (X HI ). Figure 1 shows the neutral fraction along one line-of-sight in our simulation.Along each line of sight we search for regions with X HI values larger than a given threshold, e.g. the region indicated by the red rectangle in the lower panel of Figure 1.We define the edges of these highly neutral regions to be at the intersection points between the threshold and the line of sight, e.g. the left and right sides of the red rectangle, which also corresponds to the "piece" of the LLS identified.The column density of the LLS is then simply the integral of the neutral hydrogen fraction along the line of sight between the two intersection points.However, we note that such a method for determining the HI column density would not work reliably in the Lyα forest.
We use X HI = 0.001 as our fiducial threshold for highly neutral regions.This choice is somewhat arbitrary and needs to be tested.Using a different threshold may change the column density of each LLS.In order to explore the artifacts induced by our choice of this threshold, we consider three different values for X HI (0.0001, 0.001, and 0.003) and compare column density distributions for the three cases.
To calculate the column densities of LLSs, we count the LLSs in bins of natural log of the column densities, normalized by the total comoving path length sampled in the simulation (∆l) and the bin size.Historically, the observational data has been reported per unit "absorption distance" ∆X.It is mostly a historical artifact and is defined as whereas we measure the comoving path length ∆l in the simulation.The conversion between ∆l and ∆X is straightforward: We illustrate how the LLS column density distribution changes with varying the threshold X HI in Figure 2. The primary takeaway from the Figure is that the number of LLSs with N HI > 3×10 17 cm −2 is robust to the threshold choice.

Galaxy Associations and Visualization of Lyman limit systems
A primary goal of this paper is to have a clearer understanding of where LLSs are located with respect to surrounding galaxies.We use the Rockstar halo finder (Behroozi et al. 2013) to identify dark matter halos in our simulations.The well-defined center of a dark mat-ter halo (i.e. the minimum of the gravitational potential) provides an accurate proxy for the location of each hosted galaxy, enabling spatial associations with LLSs.
We define the "association" between LLSs and nearby galaxies using both the spatial locations and virial radii (R vir ) of the dark matter halos hosting these galaxies.We quantify the strength of associations using the quantity f r<n×Rvir (n|n ∈ {1, 2, 3, 4}), where r is the distance between the location of the LLS neutral fraction peak and the center of the host dark matter halo of the nearby galaxy.For example, f r<1×Rvir (1) corresponds to the fraction of LLSs' that sit inside 1 × R vir of any galaxy.Higher n values indicate weaker associations with any galaxy, since the LLS sits at larger distances with respect to the R vir of any galaxy.
Our definition for the LLS-galaxy association differs from some in the literature, e.g.Rahmati & Schaye (2014), who defined the associated galaxy as the closest detected galaxy to the LLS.(Kohler & Gnedin 2007) showed that the latter approach is not as robust: for a given pair of a LLS and a nearby galaxy, one can frequently find a galaxy with a lower mass that is closer to the LLS, since lower mass galaxies are more numerous.And, as long as the numerical resolution of a simulation is sufficient to resolve such lower mass galaxies, these objects will be identified as the associated object to the LLS.As the result, the nearest neighbor identification is dependent on numerical resolution.In contrast, defining the LLS association with respect to the virial radius of the candidate galaxy automatically excludes all smaller galaxies within the virial radius of a larger one (i.e.satellite galaxies).
In addition, we note that the association of LLSs with satellite galaxies can be rather ambiguous.First, satellites can be tidally truncated, so their virial radius is not always meaningful.Second, the satellites may have tidal tails that contain material from the satellite galaxy but are bound to the host; it is ambiguous whether the material is associated with a satellite or the host.We choose to associate LLSs with the host galaxy through the virial radius association criterion in order to avoid this ambiguity.Additionally, we consider this approach more (although not completely) numerically robust and less dependent on numerical resolution.Granted, if a given LLS is not associated with any galaxy for a given value of n, this does not imply that there is no association for a lower mass galaxy that is not resolved in the simulation.
To help visualize the LLSs, we generate gas density projection plots around some LLSs in simulation box snapshots using the "yt" simulation analysis package (Turk et al. 2011).In these projection plots, one can visually identify structures surrounding the LLS including galaxies and connective filaments.

Statistics of LLSs
We use the column density distribution of LLSs defined in Section 2.2 to compare column density distribution of LLSs for CROC simulations.In Figure 3, the left panel shows the individual distributions from each of the three different realizations of the CROC simulations as an estimate of the sample variance in a finite simulation volume.Here, we see that the column density distributions of LLSs in CROC exhibit an overall trend where the LLS number density decreases with increasing column density.Overall, the three simulation boxes produce similar column density distributions within the level of the numerical noise with which we are able to measure them.
The right hand panel shows the column density distribution of LLSs from the combined three CROC boxes, alongside the fits to observational data from Crighton et al. ( 2019) at z = 2.4 and z = 4.4.We see that our simulation result and the observational fits have similar shapes, with CROC simulations having more LLSs with 16 < log 10 (N HI ) < 19.The differences between the column density distributions of Crighton et al. (2019) data and the average distribution across all CROC boxes are reasonable given the limited observational data and modeling uncertainties.First, we expect potential variations up to an order of magnitude for the distribution at the lowest column density bins (log 10 N HI < 17.5) that depend on the neutral fraction threshold criterion for identifying LLSs in our simulations (see Figure 2).Next, we see that the observational data in the highest column density bins (log 10 N HI ≳ 20) do not exhibit any strong evolutionary trend.The measurements from CROC are consistent with this behavior.The overall larger amplitude in column density distribution of LLSs in CROC at z ∼ 6 is consistent with the suggested evolution from observations indicating an increasing normalization with redshift.Some increase in normalization with redshift is plausible since we expect more dense absorbers to exist during reionization, with a larger fraction of neutral hydrogen present compared with lower redshift.It is sufficient for our purpose that our simulation measurements capture the dominant trends in the data.

Association of LLSs with Galaxies
A primary goal of this paper is to provide theoretical predictions of how LLSs at z ∼ 6, near the end of cosmic reionization, are associated with galaxies and nearby structures.In order to quantify the association of Lyman limit systems and galaxies in our simulation, we show in Figure 4 4.The fraction of Lyman limit systems within n×Rvir of any galaxies, with n ∈ 1, 2, 3, 4. Blue, red, orange, and green indicate respective fractions of Lyman limit systems within that distance of any galaxy.We see that higher column density systems have progressively stronger associations with galaxies.
est range, N HI < 10 17 cm −2 , only 10% of these systems are within 1 × R vir of a galaxy center and only 40% are within 2 × R vir .
Of the population of highest column density LLSs (i.e.those with 10 19 cm −2 < N HI < 10 20 cm −2 ), 80% are within 2 × R vir of some galaxy.Figure 5 shows the projected neutral hydrogen density around a LLS with the column density of N HI = 10 19.4 cm −2 , located within 1 × R vir of the virial radius of a galaxy.The left panel shows the projection in the y-z plane, and the right in the x-z plane.The depth of projection corresponds to the estimated size of the LLS (∼10 kpc).
We indicate the LLS location with a green star, and surrounding dark matter halo centers with their virial radii respectively in black crosses and black circles.The galaxy associated with the LLS has a dark matter halo mass and virial radius of M = 2.8 × 10 9 M ⊙ and R vir = 6.8 kpc respectively; we show it in red.
The remaining 20% of the highest column density LLSs are outside of 2 × R vir of any galaxy.While these are not closely associated with any galaxy, we find that they are associated with overdense filamentary structures in the large-scale surrounding environments.Figure 6 shows the projected neutral hydrogen map surrounding such a system, with similar projections and color mapping as Figure 5. Here, we see that the LLS sits in a filamentary bridge connecting galaxies, but lies well outside twice the virial radius of nearby galaxies.

Physical Properties of LLSs
In this section, we investigate the relationship between the optical depth of LLSs and their various physical properties in CROC simulations.The LLS physical properties of interest here are: optical depth, characteristic size, neutral fraction, and density.We measure Figure 5. Projected neutral hydrogen maps centered around a Lyman limit system with NHI = 10 19.4 cm −2 .The projection depth of the graph is ∼ 10 kpc, corresponding to the characteristic size of the Lyman limit system.Left: zy-projection plane.Right: xz-projection plane.The green star indicates the position of the Lyman limit system; the black circles center on galaxies and indicate their corresponding dark matter halo virial radii; the red circle corresponds to the galaxy that contains the Lyman limit system within its virial radius.Most of the high column density Lyman limit systems reside within one virial radius of a galaxy.
Figure 6.Projected neutral hydrogen maps centered around an example Lyman limit system weakly associated with any galaxy.The projection depth is ∼ 10 kpc.Markers indicate the same properties as in Figure 5.This Lyman limit system is well outside of 2 times the virial radius of any galaxy.For illustration purposes, we annotate galaxies up to 5 virial radii away from the Lyman limit system.all physical properties for each LLS at the peak neutral fraction location of the LLS along the line of sight.
Optical depth describes opacity of the LLSs and is directly proportional to the column density.Optical depth is defined as: where σ LL = 6.3 × 10 −18 cm 2 is the hydrogen photoionization cross section at the Lyman limit.
We define the characteristic size of a LLS as: The blue, red, orange, and green contours respectively correspond to systems within 1, 2, 3, and 4 times the virial radius of any galaxy.For each color, the contours correspond to 1-σ and 2-σ regions of the distribution.
the simulated LLSs.The color bar on top of each plot represents the number of data points inside each bin.The horizontal dashed lines show the τ range corresponding to the LLS column density range, although we note that the numbers are arbitrarily defined.We see from the middle and right panels that τ increases with the peak neutral fraction and the peak gas density of the LLSs.This means that higher column density LLSs are more neutral and denser at the center.However, there is not much of a trend between τ and L eff -i.e., at a fixed column density (or, equivalently, τ ), some LLSs are compact dense gas clumps while others are more spatially extended but less dense.Most LLSs, however, have sizes between 1 and ∼5 physical kpc.Note, that for LLSs with τ ≥ 100, their neutral fractions are close to unity.We also know that these high column density LLSs are strongly associated with nearby galaxies from Figure 4.This means that despite their proximity to ionizing sources, high column den-sity LLSs still have highly neutral centers that are not affected by the ionizing photons.
It is also instructive to separate LLSs based on the strength of their association to galaxies, and plot the physical quantities in the bottom panel of Figure 7.The colors we use here are the same as in Figure 4. Blue, red, orange, and green correspond to systems within 1, 2, 3, and 4 times the virial radius of any galaxy.The contour regions correspond to the 1-σ and 2-σ regions of the distribution for that subset of the LLSs.We chose to fill in the blue and red contours to visually clarify the distribution at each association strength.The strength of association with galaxies does not appear to have any significant effect on the physical properties of LLSs at z ∼ 6.

Shapes and Orientation of Lyman limit systems
LLSs at a fixed column density can come from both (a) systems of lower average column density seen along a particular direction with higher gas column density .Distribution of shape axis ratios measured from our Lyman limit system samples.Both b/c and a/c peak near 1.3, indicating that most of the systems have gas distributions that deviate from spherical.The ratio of the semimajor axis to the semi-minor axis, a/c exhibits a longer tail and skews large, corresponding to the particularly elongated systems.
and (b) systems of higher average column density seen along a particular direction with low gas column density.Such interdependencies imply that the LLS column density distribution quantified in Figure 3 reflects both the distribution of LLSs orientations in space and the distribution of their shapes.This motivates us to study the shapes and orientations of LLSs.
In order to explore the distribution of shapes of LLSs and any potential bias associated with their nonsphericity, we calculate the moment of inertia tensors of all simulated LLSs within 5 proper kpc from the location of the density peak in the LOS.We then compute the axis ratios of equivalent ellipsoids (i.e. an ellipsoid with the same moment of inertia tensor).We choose 5 kpc based on the left panel of Figure 7, which shows that L eff ≲ 5 pkpc for the vast majority of all LLSs.We show the resulting distribution of axis ratios of equivalent ellipsoids for our simulated LLS sample in Figure 8.
Since LLSs are not completely spherical, lines of sight through an LLS may only capture a threshold density peak with preferred orientations to the LLS principle axes.In particular, if a line of sight aligns with the longest principal axis of the LLS (which corresponds to the smallest eigenvalue of the moment of inertia tensor), the observed column density will bias high compared to the average for any possible line of sight.Specifically, ⟨N HI ⟩ ∝ X HI (1 + δ)L eff defined by the physical quantities shown in Figure 7.In Figure 9, we show the distribution of (cosines of) angles, cos θ, between the longest principal axis of the moments of inertia tensor and the LOS direction for all LLSs and plot this against the col- umn density.There is a mild trend of more LLSs with smaller angles, resulting in the mean cos θ being slightly lower than 0.5.Figure 10 shows the distribution of cos θ for all LLSs with N HI > 10 16 cm −2 (a more quantitative representation of the y-axis histogram from Figure 9).This is, however, what is expected for approximately ellipsoidal shapes.Namely, for a homogeneous ellipsoid with axes a and c (here, for simplicity, we adopt b = c) a column density along an angle θ with respect to the major axis is N = n × l, where n is the density and l is determined from the equation, or If the distribution of LLSs over the column densities and angles θ (shown in Figure 9) is f (N, µ) (where we use a short-hand notation µ ≡ cos θ), then, where N a ≡ n × a is the column density along the major axis, f (N a ) is the column density distribution with respect to N a , and Ñ (N a , µ) is the function that gives a value of N for a given values of N a and µ.Here we also assume that f (N a ) is independent of µ, since f (N a ) is intrinsic to the system and, hence, independent of the location of the observer, i.e. f (N a , µ) = Cf (N a ), with C = 1 since µ goes from 0 to 1.This derivation serves as a mere illustration, so here we also assume that the ratio c/a is constant, so that there is no additional distribution over c for a fixed a.
From Equation (7) it immediately follows that with α ≡ c/a.Here Ña (N, µ) is the function that gives a value of N a for a given values of N and µ, i.e. the inverse function to Ñ (N a , µ).We consider this model distribution over µ for N HI > 10 16 cm −21 and the value of α = 0.54, the average value from Figure 8, and overplot it in Figure 10 (blue dashed line) to compare against the distribution in alignment from our simulated sample (black solid line).The agreement between the simulation results and the simple model suggests that there is no significant orientation-dependent bias present.The small deviations between the two lines cannot be con-sidered significant given the simplicity of the analytical model.

SUMMARY AND DISCUSSION
We use the state-of-the-art reionization simulation, CROC, to investigate the statistics, environment and physical properties of Lyman limit systems at z ∼ 6.We cast a total of 3000 lines-of-sight to simulate synthetic quasar sight lines.Along these sight lines, we are able to identify LLSs in the simulated IGM.
Our main results are: 1. Galaxy distribution largely dictates where we will find most of the LLSs.The highest column density self-shielded absorbers, such as damped Lyα systems, have the strongest association with galaxies.These are typically located within one virial radius of a galaxy.
2. Lower column density LLSs may be located further from the nearest galaxy, but even the systems with column densities as low as 10 16 cm −2 are mostly (80%) located within four virial radii from the nearest galaxy, with the majority (50%) of LLSs being located within just two virial radii from the nearest galaxy.
3. Rare LLSs that are not located close to any galaxy are preferentially found along the filament connecting more distant galaxies.
4. The column density of LLSs strongly correlates with the gas density and the hydrogen neutral fraction.Their characteristic sizes, however, show little correlation.
5. The gas distributions around LLSs tend toward a slightly elongated gas distribution, with the peak distribution in axis ratio being a/c ∼ 1.4.Despite this slight preference, our mock observed selection of LLSs remain unbiased: they are consistent with every LLS being sampled fairly along all possible directions.
One of the physical properties of the LLSs that we measure is the effective size.The effective size of the LLSs provides a reference point.The effective size of the LLSs is comparable to the sizes of galactic disks and smaller than virial radii of galactic halos at this redshift.Therefore, it is likely that LLSs are not largely ISM structures as DLAs are, which would make the effective size comparable to the disk scale height.But, LLSs could be tidal tails from disrupted satellites or intergalactic filaments, as Figures 5 and 6 seem to indicate.
In our fourth point above, we note that the strong correlation with density is the result of the effective size being weakly correlated with the column density.In general, column densities and physical densities do not need to be correlated.For example, if all LLSs were coming from pancakes or filaments, there would be no or weak n 1/3 correlation with density.
We emphasize that the last point is nontrivial.For a random line of sight to exhibit a density peak above our selection threshold, we would expect some dependence on the exact shape of the column density distribution.In our mock observed distribution, it is just as likely for the more numerous lower column density systems to be seen along their major axis directions as for rarer higher column density systems to be seen along a direction more aligned with their minor axes.

Figure 1 .Figure 2 .
Figure 1.Neutral fraction along an example line of sight.The top panel shows two identified "highly neutral" regions in a single line-of-sight (red peaks), and the bottom panel zooms into the second peak with the red rectangle indicating the region of the line-of-sight above a threshold neutral fraction of XHI = 0.001 used to identify such regions.

Figure 3 .
Figure 3. Left: column density distribution of LLSs for three different 40 h −1 cMpc boxes.Right: comparison between the observational fit at z = 2.4, z = 4.4 (Crighton et al. 2019) and our simulations.

Figure 7 .
Figure 7. Lyman limit system optical depth vs various physical properties of the Lyman limit systems in CROC.From left to right: characteristic size, neutral fraction, and gas overdensity, respectively.The top panels are two-d histogram plots that visualize the number density distribution of data points for each LLS.The color bar above indicates the number of data points in each square bin.The bottom panels show the same distribution of systems in subsets of association strength with any galaxy.The blue, red, orange, and green contours respectively correspond to systems within 1, 2, 3, and 4 times the virial radius of any galaxy.For each color, the contours correspond to 1-σ and 2-σ regions of the distribution.
Figure8.Distribution of shape axis ratios measured from our Lyman limit system samples.Both b/c and a/c peak near 1.3, indicating that most of the systems have gas distributions that deviate from spherical.The ratio of the semimajor axis to the semi-minor axis, a/c exhibits a longer tail and skews large, corresponding to the particularly elongated systems.

Figure 9 .
Figure9.Blue points: Alignment between the line of sight and the principal axis of each LLS, cos θ, as a function of the measured column density for that system, log 10 NHI .Black line: Regression line with shaded regions corresponding to the confidence interval.Blue histograms on axes: Corresponding marginal histograms for each quantity.There is no apparent trend in alignment with the column density.
Figure10.Probability density function of cos(θ) for our simulated sample of LLSs (NHI > 10 16 cm −2 , the black line -note that it is this cut that makes this distribution µdependent).The blue line shows a simple model from Equation (7).This figure serves as a more quantitative version of the y-axis histogram from Fig.9.