The studies of disordered heterogeneous media and galaxy cosmology share a common goal: analyzing the disordered distribution of particles and/or building blocks at microscales to predict physical properties of the medium at macroscales, whether it be a liquid, colloidal suspension, composite material, galaxy cluster, or entire Universe. The theory of disordered heterogeneous media provides an array of theoretical and computational techniques to characterize a wide class of complex material microstructures. In this work, we apply them to describe the disordered distributions of galaxies obtained from recent suites of dark matter simulations. We focus on the determination of lower-order correlation functions, void and particle nearest-neighbor functions, certain cluster statistics, pair-connectedness functions, percolation properties, and a scalar order metric to quantify the degree of order. Compared to analogous homogeneous Poisson and typical disordered systems, the cosmological simulations exhibit enhanced large-scale clustering and longer tails in the void and particle nearest-neighbor functions, due to the presence of quasi-long-range correlations imprinted by early Universe physics, with a minimum particle separation far below the mean nearest-neighbor distance. On large scales, the system appears hyperuniform, as a result of primordial density fluctuations, while on the smallest scales, the system becomes almost antihyperuniform, as evidenced by its number variance. Additionally, via a finite-scaling analysis, we compute the percolation threshold of the galaxy catalogs, finding this to be significantly lower than for Poisson realizations (at reduced density ηc=0.25 in our fiducial analysis compared to ηc=0.34), with strong dependence on the mean density; this is consistent with the observation that the galaxy distribution contains voids of up to 50% larger radius. However, the two sets of simulations appear to share the same fractal dimension on scales much larger than the average intergalaxy separation, implying that they lie in the same universality class. We also show that the distribution of galaxies is a highly correlated disordered system (relative to the uncorrelated Poisson distribution), as measured by the τ order metric. Finally, we consider the ability of large-scale clustering statistics to constrain cosmological parameters, such as the Universe’s expansion rate, using simulation-based inference. Both the nearest-neighbor distribution and pair-connectedness function (which includes contributions from correlation functions of all order) are found to considerably tighten bounds on the amplitude of quantum-mechanical fluctuations from inflation at a level equivalent to observing 25 times more galaxies. The pair-connectedness function in particular provides a useful alternative to the standard three-particle correlation, since it contains similar large-scale information to the three-point function, can be computed highly efficiently, and can be straightforwardly extended to small scales (though likely requires simulation-based modeling). This work provides the first application of such techniques to cosmology, providing both a novel system to test heterogeneous media descriptors and a tranche of new tools for cosmological analyses. A range of extensions are possible, including implementation on observational data; this will require further study on various observational effects, necessitating high-resolution simulations.6 MoreReceived 16 July 2022Revised 2 January 2023Accepted 10 January 2023DOI:https://doi.org/10.1103/PhysRevX.13.011038Published by the American Physical Society under the terms of the Creative Commons Attribution 4.0 International license. Further distribution of this work must maintain attribution to the author(s) and the published article’s title, journal citation, and DOI.Published by the American Physical SocietyPhysics Subject Headings (PhySH)Research AreasLarge scale structure of the UniverseGravitation, Cosmology & AstrophysicsInterdisciplinary Physics