We present a theoretical framework in which bound stellar clusters arise naturally at the high-density end of the hierarchy of the interstellar medium (ISM). Due to short free-fall times, these high-density regions achieve high local star formation efficiencies, enabling them to form bound clusters. Star-forming regions of lower density remain substructured and gas-rich, ending up unbound when the residual gas is expelled. Additionally, the tidal perturbation of star-forming regions by nearby, dense giant molecular clouds imposes a minimum density contrast required for the collapse to a bound cluster. The fraction of all star formation that occurs in bound stellar clusters (the cluster formation efficiency or CFE) follows by integration of these local clustering and survival properties over the full density spectrum of the ISM, and hence is set by galaxy-scale physics. We derive the CFE as a function of observable galaxy properties, and find that it increases with the gas surface density, from ~1% in low-density galaxies to a peak value of ~70% at densities of ~10^3 Msun pc^-2. This explains the observation that the CFE increases with the star formation rate density in nearby dwarf, spiral, and starburst galaxies. Indeed, comparing our model results with observed galaxies yields excellent agreement. The model is applied further by calculating the spatial variation of the CFE within single galaxies. We also consider the variation of the CFE with cosmic time and show that it increases with redshift, peaking in high-redshift, gas-rich disc galaxies. It is estimated that up to 30-35% of all stars in the Universe once formed in bound stellar clusters. We discuss how our theory can be verified with Gaia and ALMA, and provide implementations for future theoretical work and for simulations of galaxy formation and evolution.

On the fraction of star formation occurring in bound stellar clusters