The escape of LyC photons emitted by massive stars from the dense interstellar medium of galaxies is one of the most significant bottlenecks for cosmological reionization. The escape fraction shows significant scatter between galaxies and anisotropic spatial variation within them motivating further study of the underlying physical factors responsible for these trends. We perform numerical radiation hydrodynamic simulations of idealized clouds with different gas surface densities compactness Sigma sim 102105  Modot rmpc2 meant to emulate star clusterforming clumps ranging from conditions typical of the local Universe to the high ISMpressure conditions more frequently encountered at high redshift. Our results indicate that dense compact star clusters with Sigma gtrsim 104  Modot rmpc2 efficiently leak LyC photons with cloudscale luminosityweighted average escape fractions gtrsim 80 as opposed to lesssim 10 for Sigma sim 100  Modot rmpc2. This occurs due to higher star formation efficiencies and shorter dynamical timescales at higher Sigma the former results in higher intrinsic LyC emission and the latter implies rapid evolution with a burst of star formation followed by rapid gas dispersal permitting high LyC escape well before the intrinsic LyC emission of stellar populations drop sim 4  mathrmMyr. LyC escape in dense clouds is primarily facilitated by highly ionized outflows driven by radiation pressure on dust with velocities  sim 3 times the cloud escape velocity. We also vary the assumed dust abundances Zrmd and find a very mild increase sim 10 in the escape fraction for sim 100 lower Zmathrmd. Our results suggest a scenario in which localized compact bursts of star formation in galaxies are disproportionately productive sites of LyC leakage.
['Astrophysics - Astrophysics of Galaxies']
Bursts of star formation and radiationdriven outflows produce efficient LyC leakage from dense compact star clusters
{'Bursts of star formation and radiation-driven outflows produce efficient LyC leakage from dense compact star clusters': 'Shyam H. Menon , 1, 2 Blakesley Burkhart , 1, 2 Rachel S. Somerville , 2 Todd A. Thompson , 3, 4, 5 and Amiel Sternberg 2, 6, 7  \n1  \nDepartment of Physics and Astronomy, Rutgers University, 136 Frelinghuysen Road, Piscataway, NJ 08854, USA 2 Center for Computational Astrophysics, Flatiron Institute, 162 5th Avenue, New York, NY 10010, USA 3 Center for Cosmology and Astro-Particle Physics, Ohio State University, 140 W. 18th Ave, Columbus, OH 43210, USA 4 Department of Physics, Ohio State University, 191 W. Woodruff Ave, Columbus, OH 43210, USA 5 Department of Astronomy, Ohio State University, 140 W. 18th Ave, Columbus, OH 43210, USA 6 School of Physics and Astronomy, Tel Aviv University, Ramat Aviv 69978, Israel 7 Max-Planck-Institut fur extraterrestrische Physik (MPE), Giessenbachstr., 85748 Garching, Germany', 'ABSTRACT': 'The escape of LyC photons emitted by massive stars from the dense interstellar medium of galaxies is one of the most significant bottlenecks for cosmological reionization. The escape fraction shows significant scatter between galaxies, and anisotropic, spatial variation within them, motivating further study of the underlying physical factors responsible for these trends. We perform numerical radiation hydrodynamic simulations of idealized clouds with different gas surface densities (compactness) Σ ∼ 10 2 -10 5 M ⊙ pc -2 , meant to emulate star cluster-forming clumps ranging from conditions typical of the local Universe to the high ISM-pressure conditions more frequently encountered at high redshift. Our results indicate that dense compact star clusters with Σ ≳ 10 4 M ⊙ pc -2 efficiently leak LyC photons, with cloud-scale luminosity-weighted average escape fractions ≳ 80% as opposed to ≲ 10% for Σ ∼ 100 M ⊙ pc -2 . This occurs due to higher star formation efficiencies and shorter dynamical timescales at higher Σ; the former results in higher intrinsic LyC emission, and the latter implies rapid evolution, with a burst of star formation followed by rapid gas dispersal, permitting high LyC escape well before the intrinsic LyC emission of stellar populations drop ( ∼ 4 Myr). LyC escape in dense clouds is primarily facilitated by highly ionized outflows driven by radiation pressure on dust with velocities ∼ 3 times the cloud escape velocity. We also vary the (assumed) dust abundances ( Z d ) and find a very mild increase ( ∼ 10%) in the escape fraction for ∼ 100 lower Z d . Our results suggest a scenario in which localized compact bursts of star formation in galaxies are disproportionately productive sites of LyC leakage. We briefly discuss possible observational evidence for our predictions and implications for cosmic reionization.  \nKeywords: Stellar feedback (1602) - High-redshift galaxies(734) - Young star clusters(1833) Reionization(1383) - Interstellar medium(847)', '1. INTRODUCTION': "The Epoch of Reionization (EoR) marks the most prominent phase transition of the Universe, when the intergalactic medium (IGM) comprised of neutral Hydrogen (H I ) was (re-)ionized by the first stars, galaxies and quasars (Loeb & Barkana 2001; Gnedin & Madau 2022). The timing of reionization has been well constrained through spectroscopic studies of high redshift (highz ) galaxies and quasars (e.g., McGreer et al. 2015), absorption lines via the Lymanα forest (Fan et al. 2006; McQuinn 2016; Burkhart et al. 2022), and the integrated  \noptical depth of Thomson scattering to the cosmic microwave background (Planck Collaboration et al. 2020), placing the conclusion of the EoR at z ∼ 6 (Robertson et al. 2015; Bosman et al. 2022). However, there is still significant debate on the dominant underlying sources of reionization, with data from the recently launched James Webb Space Telescope proving fertile ground for testable predictions (Robertson 2022).  \nSpecifically, the relative role of the different sources of UVradiation - i.e. low-mass and high-mass star-forming galaxies and the accretion disks of active galactic nuclei  \n(AGN) - is heavily debated. This is primarily quantified through the combination of the relative numbers (UV luminosity functions) of these sources, their intrinsic hydrogen-ionizing Lyman continuum (LyC) photon production efficiencies ( ζ ion ), and their LyC escape fraction ( f esc , LyC ). While AGN have been identified as potential contributors to reionization due to the copious amounts of high-energy radiation they emit and their high escape fractions (Madau & Haardt 2015), their inferred abundances at highz however, have proved difficult to reconcile with the required LyC photon budgets to conclude reionization by z ∼ 6 (Fan et al. 2006; Qin et al. 2017; Mitra et al. 2018; Parsa et al. 2018; Hassan et al. 2018; Dayal et al. 2024) 1 . This suggests that starforming galaxies are the dominant contributors to reionization (Dayal & Ferrara 2018). At highz , low-mass galaxies are more numerous, and harbor weaker potential wells permitting stellar feedback to easily eject gas and leak LyC photons (e.g., Paardekooper et al. 2013; Wise et al. 2014; Xu et al. 2016; Finkelstein et al. 2019; Bera et al. 2023); the feedback, however, would also maintain low star-formation rates and therefore potentially meager LyC photon budgets (Kimm et al. 2017). On the other hand, massive galaxies have higher star formation rates, and thereby a larger flux of LyC photons (Sharma et al. 2016; Naidu et al. 2020); conversely, their large reservoirs of gas and dust, and steep potential wells have been argued as inhibiting factors for LyC escape. Some theoretical studies also find a relatively weak dependence of escape fraction on host galaxy mass (Ma et al. 2015) and/or an evolution of this dependence with redshift (Yeh et al. 2023)  \nWhile galaxy-level studies provide key statistical constraints, young massive stars produce LyC photons while embedded in dense clouds with scales of tens of parsecs. Thus, the most prominent bottleneck to LyC escape is clearing out the dense gas surrounding young stellar populations (Paardekooper et al. 2015; Trebitsch et al. 2017). These scales are unresolved in studies probing indirect tracers of LyC escape at highz , and even in (rare) low-redshift analog LyC leakers (Leitet et al. 2013; Leitherer et al. 2016; Izotov et al. 2016; Flury et al. 2022). State-of-the-art EoR cosmological simulations only marginally resolve these scales at best (e.g., Rosdahl et al. 2022) and typically adopt subgrid recipes for LyC escape, with galaxy-averaged f esc , LyC values highly sensitive to the adopted subgrid escape fraction (e.g., Kostyuk et al. 2023). This is highly rele-  \nvant because LyC escape in resolved galaxy populations and/or lensed galaxies have been found to be highly anisotropic and to show significant variation within the galaxy (e.g., Rivera-Thorsen et al. 2019; Choi et al. 2020; Ramambason et al. 2020). This implies that a galaxyaveraged estimate along a given (projected) line of sight might not be representative of the full 3D LyC escape that takes this spatial variation into account - a subtlety that numerical simulations have been increasingly able to identify with improved numerical resolution (Yeh et al. 2023; Choustikov et al. 2024). Another point to consider is that the timescales over which stellar populations emit LyC photons are short ( ≲ 10 Myr), and the interaction of stellar feedback - especially pre-supernova mechanisms such as radiation and winds - with the ISM over these timescales are likely crucial in the context of the instantaneous value of f esc , LyC . These spatial and temporal variations likely contribute to the significant scatter found in identified correlations of f esc , LyC with host galaxy properties and their (aperture-averaged) emission-line features - such as the [O III ]/[O II ] equivalent width (EW) ratios and Ly α emission, among others (e.g. Flury et al. 2022). All of these factors indicate that probing the physics of massive star (cluster) formation and feedback-regulated LyC escape on small spatial and temporal scales can provide insights relevant for identifying the underlying environmental properties that favor LyC escape.  \nIn this study we conduct idealized numerical simulations of turbulent self-gravitating clouds designed to emulate clumps/regions of star formation within a galaxy. We model a range of clouds with different initial conditions to infer the physical factors that determine the efficiency of LyC escape on cloud scales. Specifically, we probe a range of initial cloud gas surface densities Σ cloud ∼ 10 2 -10 5 M ⊙ pc -2 to mimic ISM conditions typical of clumps ranging from weakly star-forming galaxies in the local Universe, to the high-pressure gasrich conditions more frequently encountered in the high redshift Universe. We extend the upper limits of Σ cloud probed in similarly designed numerical simulations in the literature (Howard et al. 2018; Kim et al. 2019; He et al. 2020; Kimm et al. 2022), by ∼ 2 orders of magnitude motivated by recent findings from JWST revealing young, dense, compact super star clusters (Pascale et al. 2023; Pascale & Dai 2024; Mowla et al. 2024; Adamo et al. 2024) and extremely blue, UV-luminous, compact galaxies (Finkelstein et al. 2023; Casey et al. 2023; McLeod et al. 2024; Robertson et al. 2023; Harikane et al. 2023; Morishita et al. 2023) in the EoR with stellar surface densities ≳ 10 4 -10 5 M ⊙ pc -2 . We show below that these 'extreme' high density star formation envi-  \nhighly productive sites of LyC escape due to their high star formation efficiencies, short dynamical timescales, and radiatively driven outflows that facilitate the formation of density-bounded ionized channels.  \nThe paper is organized as follows: In Section 2 we describe the numerical methods used, the initial conditions of our clouds and the range of parameter space we explore. In Section 3 we present the evolution of our clouds, compare the temporal trends of f esc , LyC we obtain over our parameter space, and discuss the physics driving the escape of LyC photons. In Section 4 we discuss some possible implications of our results, compare with observational findings that seem broadly consistent with our results, and identify the caveats with our findings imposed by our numerical setups. Finally, in Section 5 we conclude with a brief summary of our results.", '2.1. Numerical Methods': 'We solve the equations of self-gravitating radiation hydrodynamics for all our simulations. We use the FLASH magneto-hydrodynamics code (Fryxell et al. 2000; Dubey et al. 2008) for our simulations, with the explicit Godunov method in the split, five-wave HLL5R (approximate) Riemann solver (Waagan et al. 2011) for the hydrodynamics. The Poisson equation for the self-gravity is solved using a multi-grid algorithm implemented in FLASH (Ricker 2008). Sink particles are used to follow the evolution of gas at unresolved scales, the formation of which is triggered when gas properties satisfy a series of conditions to test for collapse and star formation (Federrath et al. 2010a). Specifically, sinks form when the local gas density exceeds ρ sink = πc 2 s /Gλ 2 J - where c s is the thermal sound speed, and lambda J the Jeans length - and the gas is gravitationally bound and collapsing over a zone of radius 2.5 times the minimum cell resolution (i.e. the accretion radius r sink ). Once created, a sink particle can accrete gas, if computational cells within r sink exceed ρ sink , and satisfy the additional conditions above. Gravitational interactions of sink particles with gas and other sinks are included, and a second-order leapfrog integrator is used for sink particle dynamics (Federrath et al. 2010a, 2011). To model the radiative transfer and the associated energy and momentum transfer to gas, we use the Variable Eddington Tensor-closed Transport on Adaptive Meshes method ( VETTAM ; Menon et al. 2022b). VETTAM uses the VET closure obtained with a time-independent ray-trace solution (Buntemeyer et al. 2016) to close the moment equations and an implicit global temporal update for the radiation moment equations for each band,  \naccounting for the coupling between the bands due to UV radiation being reprocessed to the IR by dust. We solve for the (non-equilibrium) ionisation state of hydrogen set by photoionization-recombination balance in an operator-split fashion with subcycling 2 , similar to the approach presented in Kim et al. (2017); the photoionization rate is obtained from the solution for the LyC radiation field from VETTAM , and we use the temperaturedependent case B recombination coefficient reported in Hui & Gnedin (1997). The radiative output from sink particles is included as a smoothed source term in the moment equations, where we have tested convergence in the smoothing parameters (see Menon et al. 2022a). We used a fixed uniform grid resolution of 256 3 for all our simulations. Although VETTAM fully supports AMR, we chose a fixed modest resolution for simplicity and to enable a more controlled experiment across our parameter space. We demonstrate convergence in the SFE (within ≲ 5%) at these resolutions in similar numerical setups in Menon et al. (2023).', '2.2. Simulation setup': "Our simulation setup is very similar to that described in Menon et al. (2023, 2024) with some additional modifications due to the inclusion of photoionizationrecombination balance. We model isolated clouds (primarily) parameterized by their mass surface density Σ cloud = M cloud / ( πR 2 cloud ) M cloud where M cloud is the mass of the cloud and R cloud its radius; this results in a gas density of ρ cloud = M cloud / [(4 / 3) πR 3 cloud ]. We achieve different target values of Σ cloud = 10 2 -10 5 M ⊙ pc -2 by fixing M cloud = 10 5 M ⊙ and varying R cloud . The lower and upper ranges of Σ cloud span the conditions of GMC/star-forming clumps typical of (weakly) star-forming galaxies in the local Universe and high-redshift/starburst conditions respectively. The simulated clouds are placed in an ambient medium with a density ρ = ρ cloud / 100 within a computational domain of size L = 4 R cloud . The fluid is initialized with turbulent velocities following a power spectrum E ( k ) ∝ k -2 , with a natural mix of solenoidal and compressive modes for k/ (2 π/L ) ∈ [2 , 64], generated using methods described in Federrath et al. (2010b) and implemented as in Federrath et al. (2022). The cloud's velocity dispersion σ v is scaled so that the clouds are marginally bound, i.e., α vir = 2, where α vir is defined as  \nα vir = 2 E kin E grav = 5 R cloud σ 2 v 3 GM cloud ,  \nwith E kin = (1 / 2) M cloud σ 2 v and E grav = (3 / 5) GM 2 cloud /R cloud . Diode boundary conditions are applied to the gas quantities, allowing gas to escape through the boundaries but preventing any inflow. We set the initial gas temperature to be equal to 10 K everywhere in the domain implying a turbulent Mach number M = σ v /c s , where c s = √ k B T/m H = 0 . 28 kms -1 is the sound speed of the gas. We also assume that the fluid is a strongly coupled dust-gas mixture with a dustto-gas ratio parameterized by the abundance relative to the solar neighborhood value Z ' d = Z d /Z d , ⊙ , where Z d , ⊙ = 0 . 0081 (Weingartner & Draine 2001).  \nWe model radiation feedback in three wavelength bands: Lyman continuum (LyC; hν ≥ 13 . 6 eV), FarUV (FUV; 6 . 6 ≤ hν < 13 . 6 eV), and the infrared (IR). Photons in the LyC band can be absorbed by, and photoionize, atomic hydrogen in our simulations. Dust absorbs both LyC and FUV photons, heating the dust, which re-emits in the IR; dust can also absorb and reemit the IR photons, i.e. effectively IR scattering, which is treated by our IR band. The sources of LyC and FUV photons are sink particles that form in our simulations that physically represent a population of stars rather than individual stars; their radiative outputs are linearly proportional to the sink particle masses via luminosityto-mass ratios for both bands. We use the analytical fitting functions presented in Kim et al. (2016) 3 for the luminosity-to-mass ratios for a given total star cluster mass, which we compute as the total instantaneous stellar mass. The functions in Kim et al. (2016) were fitted to the median luminosity-to-mass ratios obtained from several realizations of calculations with the SLUG stellar population synthesis code for a range of star cluster masses ( ∼ 10 2 M ⊙ - 10 6 M ⊙ ) assuming a Chabrier IMF; this approach therefore takes into account the effect of the luminosity-to-mass decreasing (on average) due to stochastic sampling of massive stars when the mass of the cluster is low ( ≲ 10 4 M ⊙ ). For cluster masses higher than ≳ 10 4 M ⊙ , the luminosity-per-unit mass saturates to values of Ψ F UV = 943 L ⊙ M -1 ⊙ and Ψ L yC = 380 L ⊙ M -1 ⊙ respectively. However, we do not take into account the temporal evolution of the radiative outputs since most of our simulations are run for relatively short timescales over which these values do not change significantly.  \nWe use a constant (frequency-averaged) cross section of (neutral) hydrogen atoms to LyC photons of σ H = 3 × 10 18 cm 2 . The dust opacities per unit (gas) mass for the LyC and FUV bands are taken to be equal to κ FUV = 1250 Z ' d cm 2 g -1 and κ FUV = 1000 Z ' d cm 2 g -1 respectively. We use the Planck and Rosseland-mean opacities (scaled by Z ' d ) in the IR using the dust temperature-dependent model of Semenov et al. (2003). We assume the dust temperature instantly comes into radiative equilibrium with the IR radiation field; see Menon et al. (2024) for the caveats associated with this assumption, none of which should affect the results of this paper. We initialize our clouds with zero radiation energy/flux in the UV and an IR radiation field corresponding to a dust temperature of T d = 10K. We adopt Marshak boundary conditions (Marshak 1958) for the radiation with the background value set to match the initial conditions.  \nWe use a quasi-isothermal two-temperature equation of state for the thermal evolution of the gas, which by construction assumes instantaneous thermal equilibrium with two distinct (stable) phases that purely depends on the ionization state of the gas - i.e. a warm, fully photoionized phase with temperature T ion = 10 4 K, and a cold, fully neutral phase with temperature T cold = 10K. The temperature is then linearly scaled in between these two limits depending on the ionization fraction for intermediate ionization states. This two-temperature approximation has been adopted by several works in the past (Gritschneder et al. 2009; Kim et al. 2018; Menon et al. 2020), and is a reasonable approximation to first order to model the photoionized/neutral transition since the cooling timescale is very short compared to the dynamical timescale. However, this approach would not be able to produce warmer neutral gas due to other radiative processes that could act on LyC-shielded gas, such as photodissociation and the photoelectric effect. While this warmer gas would likely affect fragmentation in the cloud - which we do not focus on in this paper - it is unlikely to affect the bulk evolution of the cloud and the cloud-scale competition between feedback and gravity to zeroth order; the key physics that is relevant in these contexts are photoionization and radiation pressure on gas/dust, all of which is captured in our simulations. We do not include the feedback pathways of stellar winds, protostellar jets, and supernovae in our simulations. The latter occur on timescales longer than the cloud lifetimes; we discuss the implications of our omission of the former mechanisms in Section 4.5.  \nTable 1. Summary of our simulation suite and their initial condition parameters.  \nNotes : Columns in order indicate - Model: Σ cloud : mass surface density of the cloud given by Σ cloud = M cloud / ( πR 2 cloud ), M cloud : mass of cloud, R cloud : radius of cloud, n cloud : number density of the cloud given by n cloud = 3 M cloud / (4 πR 3 cloud m H ) where m H is the mass of atomic hydrogen, σ v : turbulent velocity dispersion of the cloud, M : turbulent Mach number given by the relation σ v /c s within the cloud where c s is the initial thermal sound speed in the cloud ( ∼ 0 . 2 kms -1 ), v esc : escape velocity of the cloud, t ff : free-fall time of the cloud, Z d : dust-to-gas ratio in units of the corresponding value in the solar neighborhood ( Z d , ⊙ ). All our simulations use a resolution of 256 3 .  \nThe suite of simulations that we present in this study is summarized in Table 1. The primary parameter that we explore in our simulation suite is Σ cloud , where we adopt values = 10 2 M ⊙ pc -2 - 10 5 M ⊙ pc -2 to progressively mimic conditions typical of GMCs in the starforming spiral galaxies that are common in the local Universe to those in the very compact, rapidly star forming galaxies expected (and now observed) in the high-redshift Universe, respectively. The velocity dispersions follow our assumption of initial marginal stability ( α vir = 2), which leads to turbulent Mach numbers between 20 and 110 for a fixed initial gas temperature of T = T cold = 10 K . We note that the escape velocities ( v esc ) are higher and the dynamical timescales ( t ff ) lower for higher Σ cloud - this feature is very important for the results we will show below. All of the above simulations are run with solar-neighbourhood-like dust abundances ( Z ' d = 1); for Σ cloud = 10 4 M ⊙ pc -2 , we also run cases with lower dust abundances of Z ' d = 0 . 1 and Z ' d = 0 . 01. This is meant to mimic the conditions of lower metallicities in the high(er) redshift Universe, especially in low-mass ( M ∗ ∼ 10 7 - 10 9 M ⊙ ) galaxies that primarily contribute to reionization. Our explored range of dust-to-gas ratios spans that occupied by high-redshift LyC leaking galaxies in cosmological simulations (e.g., Choustikov et al. 2024), observed low-redshift LyC leakers (e.g., Flury et al. 2022), and Epoch of Reionization (EoR) galaxies observed with JWST (Curti et al. 2023; Nakajima et al. 2023). We note that the lower dust abundances are only one of several physical effects that manifest differently at lower metallicities; the gas temperatures in thermal balance ( T ion and T cold ) would be  \nhigher due to weaker metal-line cooling 4 , and the UV radiative outputs of lower-metallicity stars (and therefore star clusters) higher. We chose not to include these effects here for a more controlled experiment. In Section 4.5 we discuss the caveats associated with this assumption in interpreting our Z ' d < 1 simulations as lowmetallicity counterparts, and conclude that the change to the dust opacity with metallicity is the principal effect to capture, and that the effects of the gas temperature and stellar spectra are relatively minor in the context of our results.", '3.1. General evolution of the simulations': "We discuss the general time evolution of our simulations in this section. The initial turbulent fluctuations form filamentary structures that become gravitationally unstable, and go on to collapse until they form sink particles (which represent sub-clusters of stars), the first of which forms by a timescale ≲ 0 . 5 t ff , cloud . Star formation continues to occur through both the formation of new sink particles and accretion onto existing sinks. As the sinks become more massive, their UV radiative output proportionately grows, and this starts to photoionize and provide momentum feedback (radiation pressure) on the gas surrounding the sinks. The subsequent evolution differs among runs with different Σ cloud ; lower Σ cloud (Σ cloud ≲ 10 3 M ⊙ pc -2 ) runs photoionize and drive radiation pressure-driven bubbles more successfully, whereas higher Σ cloud runs continue to accrete  \nwith relatively mild dynamical impact from the radiation. Eventually the bulk of the gas in these runs is photoionized and evacuated by a time ≲ 2 t ff , which for the low-Σ cloud run corresponds to ≲ 8 Myr. The higher Σ cloud (Σ cloud > 10 3 M ⊙ pc -2 ) runs on the other hand, still contain a significant amount of dense gas in the central regions, although there are signs of outflows through low-density channels. It is only at much later (dynamical) times ( ≳ 7 t ff ) that the bulk of the gas is driven out in an outflow, consistent with Menon et al. (2023). Broadly, the general evolution shows that at higher Σ cloud , star formation and the dispersal of the remaining gas go on for longer multiples of the cloud free-fall times.  \nThis behavior is reflected in Figure 1, which shows the integrated star formation efficiency ϵ ∗ - the fraction of the initial cloud mass that has been converted into stars - as a function of time. All our simulations show similar trends: ϵ ∗ rises rapidly on a timescale t ff , cloud , and then saturates due to radiative feedback on a timescale that depends on Σ cloud . The principal differences between models are the timescale for saturation and the saturated value of ϵ ∗ which increases sharply with Σ cloud , consistent with previous works (Fall et al. 2010; Grudi'c et al. 2018; Lancaster et al. 2021; Menon et al. 2023). If we compare these values with the results of Menon et al. (2023) where largely similar numerical methods were used, except that photoionization was not included, we find that the runs with Σ cloud ≲ 10 3 M ⊙ pc -2 here have much lower ϵ ∗ , whereas the values for the higher Σ cloud cases are more or less identical. This is consistent with theoretical expectations - the Σ cloud ≲ 10 3 M ⊙ pc -2 cases have v esc ≲ c s , ion , the ionized gas sound speed, where photoionization is expected to be the dominant feedback mechanism (Dale et al. 2012; Kim et al. 2018). The higher Σ cloud cases do not satisfy this condition, and therefore photoionization makes a negligible contribution to the feedback-matter competition, which is rather dominated by radiation pressure on dust and gas. This does not mean that the photoionization does not have any dynamical impact in our simulations; indeed, we will see below that gas needs to be photoionized/photoevaporated to enable escape of LyC photons. Additionally, the ionizing flux also affects the phase structure of radiation pressure-driven outflows.", '3.2. Escape Fractions': 'We now present the escape of UV photons and its dependence on the cloud properties. Since we explicitly solve the radiation moment equations, we can track the local absorption rate of UV photons by dust and gas -˙ E dust and ˙ E gas respectively -and compare this to the  \nnet UV radiative output in the domain to obtain the escape fraction. The net UV radiative output is obtained by summing up the luminosities of the individual sinks ( L ∗ ). We can then define the global escape fraction of the simulation at any time as  \nf esc ,ν = 1 -∑ V ( ˙ E dust ,ν + ˙ E gas ,ν ) dV ∑ ∗ L ∗ ,ν , (1)  \nwhere the subscript ν refers either to the LyC or the FUV bands. We note that gas only absorbs LyC photons 5 and therefore ˙ E gas , FUV = 0, whereas for dust it is non-zero for both bands. We can also define the fraction of radiation that is absorbed by neutral gas ( f gas ) and by dust ( f dust ) as -  \nf gas ,ν = ∑ V ˙ E gas ,ν dV ∑ ∗ L ∗ ,ν , (2)  \nand  \nf dust ,ν = ∑ V ˙ E dust ,ν dV ∑ ∗ L ∗ ,ν , (3)  \nrespectively. This maintains the closure relation f esc ,ν + f gas ,ν + f dust ,ν = 1.  \nIn Figure 2 we show f esc , LyC and f esc , FUV as a function of the time since star formation began ( t ∗ ) for all our simulations. The temporal trend is broadly similar across our simulations - an initial obscured phase followed by a rapid increase in f esc due to dispersal of gas and dust by feedback. We can see, consistent with expectations from the general evolution discussed above, that lower Σ cloud runs achieve higher escape fractions much more rapidly (in units of t ff , cloud ), and that they show a more rapid evolution from obscured ( f esc ∼ 0) to transparent ( f esc ∼ 1) to UV photons.  \nBy comparing the top-left and top-right panels in Figure 2, we can see that the escape fractions of LyC and FUV photons more or less trace each other, albeit with a slightly shorter timescale for FUV escape evident on careful inspection. However, this similar evolution of LyC and FUV escape fractions becomes less pronounced in the simulations with lower dust to gas ratios. Comparing the dashed and dotted lines with the solid lines in Figure 2, we see that there is significantly more FUV photon escape at earlier t/t ff , cloud than LyC escape. The strong effect of the dust abundance on FUV escape relative to LyC escape results from the fact that the only source of opacity for the FUV photons is from dust. On  \n/circledot  \nFigure 1. Time evolution of the integrated star formation efficiency ϵ ∗ - i.e. the fraction of initial cloud gas mass converted to stars - across the simulations in this study, shown as a function of time in units of the cloud free-fall time ( t ff , cloud ; Table 1). ϵ ∗ increases with time until (radiative) feedback from young stars is able to counteract gravity, at which point the rate of star formation ( ˙ ϵ ∗ ) decreases, and eventually ϵ ∗ saturates. The point at which this occurs depends on Σ cloud , with more efficient star formation in denser clouds.  \n<!-- image -->  \nthe other hand, LyC photons can be absorbed by both neutral gas and dust. For the LyC band, Figure 2 shows that different dust abundances yield similar LyC escape fractions at early times, while the curves start to diverge at later times. This suggests that dust absorption is more important at later times. Indeed, in Figure 2 we show that the fraction of LyC photons absorbed by dust, f dust (Equation 3), increases with time, with gas dominating the absorption at earlier times, and dust acting as the main source of opacity later in the evolution. In fact for a brief period - which depends on Σ cloud - most ( ≳ 60%) of the LyC photons emitted by the stellar population are consumed by dust rather than by gas or by escaping the computational volume (see lower panel in Figure 2). This occurs because most of the gas has already been photoionized and/or photoevaporated by ionizing radiation, and it is the surviving dust in the ionized gas that absorbs photons. However, this dusty ionized gas is also gradually evacuated via thermal pressure (at smaller Σ cloud ), or via radiation pressure on dust grains (at higher Σ cloud ), causing f dust  \nto decrease and f esc , LyC to proportionally increase. The former mechanism likely only occurs for the runs with lower Σ cloud , which have v esc ≲ c s , ion , permitting the ionized gas (and the dust along with it) to escape freely; for higher surface density clouds, radiation pressure on dust acts to evacuate the dusty ionized gas. At earlier times when both sources of opacity exist, the factor ∼ 1000 higher cross section of neutral hydrogen to ionizing UV photons over dust ensures that it acts as the dominant opacity. For runs with lower dust abundance, Figure 2 shows that even at later times, dust still plays a relatively minor role as a source of obscuration; for instance for Z d = 0 . 01 Z d , ⊙ , f dust < 0 . 2 at all times.', '3.3. Evacuation of channels by radiative feedback': 'In the previous section we presented the following qualitative scenario for LyC escape: neutral gas is increasingly photoionized by LyC radiation, and the resultant dusty ionized gas subsequently escapes via thermal expansion, or is driven out by radiation pressure on dust in dense clouds. Since the clouds are turbulent,  \n/circledot  \nFigure 2. This figure shows the possible fates of photons (see Section 3.2) emitted by the stellar populations across our simulations, shown as a function of time scaled by t ff , cloud . Top-left: escape fraction of LyC photons ( f esc , LyC ); top-right: escape fraction of FUV photons ( f esc , FUV ); bottom-left: fraction of LyC photons absorbed by neutral HI gas ( f gas , LyC ); bottom-right: fraction of LyC photons absorbed by dust ( f dust , LyC ). We do not show the latter two for the FUV band since f gas , FUV = 0 and f dust , FUV = 1 -f esc , FUV . Most LyC photons are absorbed by ionization-bounded gas at earlier times when f esc , LyC is low, followed by a shorter phase where dust can absorb up to ∼ 50% of photons as radiative feedback clears gas. Denser clouds require significantly more free-fall times to clear their surroundings and leak LyC photons, and have a longer phase where dust absorption is important.  \n<!-- image -->  \nthis likely occurs through increasingly wide channels of UV photon escape. While this qualitative picture holds across our parameter space, there is clearly a distinction in the (dynamical) timescales required for this process to play out. We explain the physics behind this behavior in this section.  \nTo systematically quantify the distinct temporal behavior across our simulation suite, we define four characteristic timescales that we can use to compare the runs:  \n- 1. t ∗ , 50 : The timescale to reach 50% of the final saturated stellar mass.\n- 2. t ∗ , 95 : The timescale to reach 95% of the final saturated stellar mass, therefore representing when star formation more or less ceases.\n- 3. t esc , 0 . 2 : The time when an instantaneous LyC escape fraction of 20% is achieved.\n- 4. t esc , 0 . 7 : The time when an instantaneous LyC escape fraction of 70% is achieved.  \nWe show the projected distribution of gas and stars at these characteristic timescales in Figure 3. We plot these characteristic timescales as a function of Σ cloud in Figure 4. We can see that star formation occurs for a slightly longer duration at higher Σ cloud . More evidently, the timescale for achieving LyC escape is significantly longer at larger Σ cloud . This indicates that the processes of photoevaporating and evacuating channels of dusty gas takes longer (in terms of free-fall timescales) at higher Σ cloud .  \nTo understand this behavior more clearly, we look at the distribution of gas columns as seen by the radiating star clusters. To obtain this, we first make the assumption that the set of distributed radiating sources can be effectively replaced with a single source placed at the centre of mass of the system, with the combined lumi-  \ncloud  \nFigure 3. Projected gas density distributions for the different Σ cloud runs (rows), plotted at the characteristic timescales defined in Section 3.3: i.e. t ∗ , 50 , t ∗ , 95 , t esc , 0 . 2 and t esc , 0 . 7 (left to right). Time in physical units ( t ), star formation efficiencies ( ϵ ∗ ) and instantaneous escape fractions ( f esc ) are annotated. Arrows indicate the velocity field.  \n<!-- image -->  \nnosity of all the sources. We then interpolate the grid quantities from our Cartesian grid onto a spherical grid centered at this point. The spherical grid is discretized with [ N r , N θ , N ϕ ] = [128 , 128 , 256], equally spaced in r ∈ [0 , 2 R cloud ], ϕ ∈ [0 , 2 π ], and cos θ ∈ [ -1 , 1], ensuring that the grid is discretized equally in solid angle. We use trilinear interpolation for this purpose. We then compute the gas column densities along each line of sight by integrating the gas density in the radial direction for a given θ, ϕ . This yields a distribution of the gas column  \nin solid angle space - i.e. f Ω (Σ). We calculate these distributions for both the total gas (Σ g ) and the neutral atomic gas (Σ HI ) separately. We refer to these as the circumcluster gas column densities.  \nTo interpret the significant differences in the LyC escape timescale with Σ cloud , we compare the circumcluster gas column densities between Σ cloud = 100 and 10 4 M ⊙ pc -2 as representative of a relatively diffuse and dense cloud, respectively. In Figure 5 we compare the distributions of Σ HI sampled at the characteristic  \ntimescales defined above. As expected, at t ∗ , 50 all sightlines are obscured by the gas, although the higher Σ cloud case has significantly higher Σ HI values. However, by t ∗ , 95 , a significant fraction of the gas has been photoionized for Σ cloud = 100 M ⊙ pc -2 , such that most sightlines are marginally optically thin. This is not the case for the high surface density model with Σ cloud = 10 4 M ⊙ pc -2 , where most sightlines are still at the mean value, with τ HI , LyC ≳ 10 6 . This indicates that the bulk photoevaporation of the gas has simultaneously been occurring along with star formation at lower Σ cloud , whereas neutral high density regions continue to survive in denser clouds even after star formation more or less ceases. This distinction in outcomes is largely due to the higher densities and higher escape velocities that accompany higher Σ cloud (for fixed mass): higher densities make it harder for the ionizing radiation field to compete with recombinations ( ∝ n 2 ); this is aggravated by the fact that the ionized gas cannot freely expand since the potential well is too deep (i.e. v esc > c s , ion ), and neither can it exert sufficient thermal pressure on the surrounding neutral gas to expand as a D-type front, since the external ram pressure is more significant ( σ v > c s , ion ).  \nInstead, it is the radiation pressure on dust and gas that is responsible for the expansion of HII regions for these denser clouds, consistent with analytic work (Krumholz & Matzner 2009; Murray et al. 2010; Draine 2011) and numerical simulations (e.g., Kim et al. 2018). This expansion lowers the ionized gas density, freeing up LyC photons (which earlier were expended balancing recombinations) to ionize additional material. This radiation pressure-driven expansion is the reason that a significant fraction of material becomes ionized by t esc , 0 . 2 for the high Σ cloud case in Figure 5. That being said, there is still a significant fraction of sightlines with high Σ HI , suggesting that radiation pressure is unable to evacuate these sightlines effectively. We verify in Figure 5 that the total gas column density (Σ g ) also shows signatures of surviving dense columns even at t esc , 0 . 2 with dust optical depths much greater than unity; in fact, inspection of Figure 2 reveals that at these times ( ∼ 8 . 5 t ff ) it is the absorption by dust that dominates LyC absorption rather than the gas, and it is the failure of stellar feedback to eject the existing high-density columns of gas (and dust) that limits the escape fraction of LyC photons. Indeed, semi-analytic models (Thompson & Krumholz 2016) predict that since the radiating sources see a distribution of column densities, there will be a corresponding distribution of outcomes for the competition between gravity and radiation pressure. This competition can be quantified with the Eddington ratio - the ratio of radiation pressure to gravitational (self-  \ntellar) forces - acting through a given column of gas.  \nWe compute the Eddington ratio f Edd for each sightline by cumulatively integrating the total (UV+IR) radiation pressure forces in that sightline, and dividing by the corresponding cumulative gravitational force. In Figure 6 we look at how f Edd varies with the corresponding Σ g among different sightlines to understand why some structures survive. We can see that there is no trend of f Edd with Σ g for Σ g ≲ 5 M ⊙ pc -2 - which corresponds to τ dust = 1 (shown as the dashed line in Figure 6) - following which it falls as Σ -1 g . This is in good agreement with simple back-of-the-envelope estimates of the Eddington ratio for a shell or clump of gas with surface density Σ g (Andrews & Thompson 2011; Thompson et al. 2015; Thompson & Krumholz 2016; Crocker et al. 2018; Blackstone & Thompson 2023). In the limit where the gas is optically thin in the UV and the gravity is dominated by the stellar component ( M ∗ ), f Edd is given by  \nf Edd , OT = κ UV L UV 4 πGM ∗ c ∼ 70 ( Ψ UV 943 L ⊙ M -1 ⊙ ) ( κ UV 1000 cm 2 g -1 ) , (4)  \nwhere we have used Ψ UV = L UV /M ∗ , where Ψ UV is the total UV luminosity-per-unit mass that we have adopted for our simulations. We have expressed the above scaled by the values we have adopted in our simulations. In the single-scattering (SS) regime - where the gas is optically thick in the UV yet optically thin in the IR - the Eddington ratio for a shell/cloud of surface density Σ g is given by (Thompson & Krumholz 2016; Wibking et al. 2018)  \nf Edd , SS = Ψ UV 4 πGc Σ g ∼ 1 ( Ψ UV 943 L ⊙ M -1 ⊙ ) ( Σ g 500 M ⊙ pc -2 ) -1 , (5)  \nwhich is ∝ Σ -1 g . This expression yields the maximum surface density that can be super-Eddington ( f Edd , SS > 1), the Eddington surface density, Σ Edd given by  \nΣ Edd ∼ 500 M ⊙ pc -2 ( Ψ UV 943 L ⊙ M -1 ⊙ ) . (6)  \nWe plot this quantity with dotted lines in Figure 6; we also denote the point where τ dust , UV > 1 to compare with the expected transition from f OT =constant to f SS ∝ Σ -1 g from the above expressions. We can see that both expectations are well satisfied in the obtained trends for our simulations at later times. At earlier times, the trends look different because f Edd , SS is dependent (linearly) on the fraction of mass in stars ( ϵ ∗ ), since gravity has contributions from both the gas and stars,  \nwhereas the radiative force is linearly proportional only to the stellar mass - this lowers Σ Edd at earlier times when ϵ ∗ is less than one. Moreover, at earlier times, the radiation is absorbed very close to the sinks, and Menon et al. (2023) have shown that this leads to cancellation of the forces from individual sinks (see also Socrates & Sironi 2013), such that the net radial force from the stellar population is much lower than the one obtained from a single source with the combined luminosity. These additional complexities become weaker at later times, resulting in the strong agreement.  \nThese figures also reveal that the high column density sightlines ( ≳ Σ Edd ) are mostly sub-Eddington, at least globally. This does not necessarily mean that the gas in those sightlines is being accreted onto the sink particles. While this is certainly true for certain sightlines, it is possible that gas could still be locally superEddington close to the radiating sources, especially considering additional sub-dominant support from the thermal pressure of ionized gas, even if the column as a whole is sub-Eddington - indeed the mass-weighted velocities as a function of Σ g show that the radial velocities are slightly larger than zero, however, insufficient to escape the cloud ( v < v esc ). However, the dynamics of highcolumn density sightlines does reveal why the timescales for evacuating gas are longer for higher Σ cloud - there is a much larger fraction of gas with Σ g > Σ Edd for these cases. The high-density gas in these sightlines is either consumed by star formation, accreted onto existing stars, and/or a small fraction of it could be locally driven outwards by feedback, all of which reduce the Σ g in that sightline. This is unlike the behaviour at lower Σ cloud where most of the sightlines are already superEddington while star formation is still ongoing.', '3.4. Properties of outflows': "The existence of super-Eddington sightlines implies that gas is driven outward by radiative forces in these channels. We show this in Figure 7 for Σ cloud = 100 M ⊙ pc -2 and 10 4 M ⊙ pc -2 for the characteristic times. Specifically, Figure 7 shows the average massweighted radial velocity ( v out ) along sightlines binned by their corresponding Σ g . Comparing the panels, we can see that the lower Σ cloud case already has the most sightlines with v out > 0 by t ∗ , 50 , consistent with the rapid bulk photoevaporation scenario introduced in the previous section, whereas the higher range of Σ g values for Σ cloud = 10 4 M ⊙ pc -2 ensures that outflows are launched only after most of the star formation has ceased. Lower-column density channels have faster outflow velocities, with the lowest columns accelerated to velocities several times the escape speed of the cloud  \nFigure 4. Characteristic timescales (defined in Section 3.3) in units of t ff , cloud for different Σ cloud . t ∗ , 50 and t ∗ , 95 indicate the times at which 50% and 95% of the final saturated stellar mass is achieved. t esc , 0 . 2 and t esc , 0 . 7 are the times at which f esc , LyC are 20% and 70% respectively. Denser clouds form stars for slightly more dynamical times than less dense clouds. They also take significantly more dynamical times to leak LyC photons. Note, however, that the dynamical times of denser clouds are significantly shorter in absolute time (Table 1).  \n<!-- image -->  \nglyph[circledot]  \n( ≳ 5 v esc ). We can also see that this dependence of v out on Σ g follows analytic estimates based on the scaling of the Eddington ratio with Σ g for radiation pressure on dust (Thompson et al. 2015). Since there is a linear anti-correlation of the Eddington ratio ( f Edd , SS ; see Equation 5) with Σ g in the single-scattering regime (Σ g ≳ 5 M ⊙ pc -2 ), the acceleration experienced by the gas, and the resulting outflow velocities, would show a trend of decreasing v out with increasing Σ g in this regime; this dependence is shown in Figure 7. On the other hand, no such dependence is expected in the optically-thin regime (Σ g ≲ 5 M ⊙ pc -2 ), where we expect optically-thin sightlines to be accelerated to a single characteristic velocity (see below).  \nTo quantify this in further detail, we follow Thompson et al. (2015) and estimate the dynamics of a spherical geometrically-thin dusty shell of surface density Σ sh exposed to a radiation source. The radial evolution for the velocity of a thin shell accelerated by radiation pressure, v sh can be written as  \nv sh dv sh dr = -GM ∗ r 2 + Ψ M ∗ 4 πr 2 c Σ sh ( 1 -e -τ UV + τ IR ) , (7)  \nwhere M ∗ is the stellar mass, which we take to dominate the gravitational potential, and where we have expressed its luminosity in terms of the luminosity-per-unit mass Ψ that we adopt. If we make the additional assumption that the shell expands without entraining a significant fraction of its original mass, then the shell column density is ∝ r -2 so that Σ sh = Σ 0 r 2 0 /r 2 , where r 0 is the  \nτ  \nτ  \nFigure 6. The PDF of total gas column densities (Σ g ) in different sightlines as seen by the center of mass of the stellar distribution at the characteristic timescales for Σ cloud = 10 4 M ⊙ pc -2 (left) and the median Eddington ratios in each bin of Σ g with the shaded region representing the 1σ variation. Dashed and dotted gray vertical lines represent the points where τ dust , LyC = 1 and the Eddington surface density respectively. This reveals a scenario where super-Eddington channels are evacuated by radiation pressure on dust, producing the low-column tails in Σ g that leak LyC photons at later times. The fraction of sightlines that are super-Eddington increase with time as gas is consumed by star formation and/or locally driven outward by radiative feedback - both of which systematically lower Σ g .  \n<!-- image -->  \n/circledot  \n/circledot  \nFigure 5. The PDF of neutral HI column densities in different sightlines as seen by the center of mass of the stellar distribution at the characteristic timescales for Σ cloud = 100M ⊙ pc -2 (left) and 10 4 M ⊙ pc -2 (right). The dashed line indicates the transition from optically thin to thick to LyC photons (i.e. τ HI , LyC=1 ). We can see that the higher Σ cloud run contains a significant fraction of high Σ HI sightlines even at later epochs, whereas the lower Σ cloud case has already photoionized the bulk of sightlines when star formation has just ceased (i.e. at t ∗ , 95 ).  \n<!-- image -->  \n/circledot  \n<!-- image -->  \n/circledot  \ninitial radius from which the shell is driven, and Σ 0 the column density at that radius. The right hand side of Equation 7 can now expressed purely as a function of r , and integrated to obtain v sh for a given initial r 0 , Σ 0 . We obtain  \nv 2 sh ( r, Σ 0 ) v 2 esc ( r 0 ) = r 0 r -1 + f Edd , 0 [ r r 0 -1 + e -τ UV , 0 -r r 0 e -τ UV , 0 ( r 0 /r ) 2 + √ πτ UV , 0 ( erf( √ τ UV , 0 ) -erf( r 0 r √ τ UV , 0 ) ) + τ IR , 0 ( 1 -r 0 r )] , (8)  \nwhere τ UV , 0 = κ UV Σ 0 and τ IR , 0 = κ IR Σ 0 are the initial UV and IR optical depths at r 0 , v esc ( r 0 ) = √ 2 GM ∗ /r 0 is the escape velocity measured at r 0 , and f Edd , 0 = Ψ / 4 πGc Σ 0 is the single-scattering Eddington ratio (Equation 5) at r 0 . κ UV and κ IR are the dust UV and IR opacities respectively; we set the former to 1000 cm 2 g -1 consistent with the FUV band value and the latter to 5 cm 2 g -1 . While κ IR is dependent on the dust temperature in our simulations, we use a constant value for simplicity; the adopted value is irrelevant as long as it is within the range of typical dust IR opacities ( ≲ 10 cm 2 g -1 ; Semenov et al. 2003), since the Eddington ratio for IR radiation pressure is less than unity for this range of values and therefore shouldn't affect this calculation (Skinner & Ostriker 2015; Menon et al. 2022a). We can use Equation 8 to compute the analytical expectation for the velocity of shells with different Σ g accelerated by radiation pressure. Following along the lines of Raskutti et al. (2017) we assume that the outflow is launched at the cloud radius ( r 0 = R cloud ), and the shell has reached the edge of the domain ( r = 2 R cloud ) - this choice implies that Σ g = Σ 0 / 4. Using these values in Equation 8 allows us to obtain v sh (Σ g ), which we overplot as the dashed lines in Figure 7. We can see that for Σ cloud = 100 M ⊙ pc -2 , the agreement between the expected velocity from Equation 8 and the simulations is poor, especially at low Σ g . This is not surprising, since photoionization is the dominant source of momentum in this regime (e.g., Kim et al. 2018), which is not considered in the analytic approximation of Equation 8. Moreover, the assumption of a single point source of radiation is less appropriate for lower Σ cloud since the distribution of star particles is more extended due to the weaker potential. However, for Σ cloud = 10 4 M ⊙ pc -2 , for which this comparison is more appropriate, the agreement between the estimate from Equation 8 and the simulation at late times (when the outflow is driven) is striking, as is evident from the lower panel of Figure 7.  \nEquations 7 and 8 also demonstrate why radiation pressure is able to accelerate gas to v out values of several times v esc . The radiation pressure when the dust is optically thick in the UV ( τ UV >> 1) is independent of r , whereas the gravitational force drops as r -2 . This implies that if the shell was super-Eddington at a given radius, the acceleration imposed on the shell is positive and increases with time, until the point where the shell becomes optically thin (Elmegreen & Chiang 1982; Thompson et al. 2015).  \nWhile the outflow velocities are a few times the cloud escape speed ( v esc ) in both Σ cloud cases, the much higher v esc for higher Σ cloud (Table 1) implies that the absolute outflow velocities are much higher - up to 100 km s -1  \nfor the low-column sightlines. To compare the average outflow speeds systematically with Σ cloud , in Figure 8 we plot ⟨ v out ⟩ Ω ,t , the outflow velocity averaged over all sightlines (Ω) and times ( t ) for our different simulations. We can see the trend of increasing ⟨ v out ⟩ Ω ,t with Σ cloud ; inspection of the overplotted v esc shows that ⟨ v out ⟩ Ω ,t ∼ 2-3 v esc for all Σ cloud . The exception are the runs with lower Z d , which have lower ⟨ v out ⟩ Ω ,t . This is because the Eddington ratio for optically thin gas is lower at lower Z d (see Equation 4). In addition, the single-scattering regime operates at Σ g ≳ 50 M ⊙ pc -2 and 500 M ⊙ pc -2 for 0 . 1 Z d , ⊙ and 0 . 01 Z d , ⊙ respectively; therefore the rapidly accelerating phase associated with this regime lasts over a smaller range in surface density, corresponding to a shorter duration in time at lower Z d . This results in the lower asymptotic velocities.", '3.5. Cumulative LyC escape': 'In the preceding sections we have shown that the higher columns in higher Σ cloud cases make it more difficult for feedback to open up channels for LyC escape. In these cases, more dynamical times are required to achieve high f esc , LyC . However, the dynamical timescales for higher Σ cloud are also significantly shorter (see Table 1), which implies that the phase of clearing out feedback-driven channels of LyC escape occurs when the stellar population is younger - so much so that this reverses the trend for the timescales for LyC escape when viewed as a function of the stellar population age. For instance, even though the Σ cloud = 10 4 M ⊙ pc -2 requires ∼ 9 t ff , cloud to achieve f esc , LyC ≳ 50% and Σ cloud = 100 M ⊙ pc -2 only requires ∼ 1 . 5 t ff , cloud (see Figure 2), the value of t ff , cloud is higher in the latter by a factor ∼ 30 (Table 1) - much more than the factor ∼ 6 longer period of cloud dispersal and subsequent LyC escape for the denser cloud. This fact results in the trend shown in Figure 9 (left panel) for f esc , LyC in terms of absolute time. In the context of reionization, it is indeed this trend that is relevant, since the LyC emission rate of a stellar population is very strongly dependent on age; this is because as massive stars die on timescales ≳ 3 Myr, the sources of LyC photons are lost, and the emission rate drops significantly.  \nTo demonstrate this, we show two quantities that take into account the time dependence of LyC emission in a stellar population in Figure 9. In the middle panel, we show the relative escape fraction, f esc , rel , given by the relation  \nf esc , rel = ( f esc , LyC f esc , FUV )( L LyC , ∗ L FUV , ∗ ) , (9)  \nwhere L LyC , ∗ and L FUV , ∗ are the intrinsic luminosities of the stellar population in the LyC and FUV bands, which  \nFigure 8. Time-averaged, mass-weighted outflow velocities for the simulations with different Σ cloud ; the points with (unfilled) circle and star markers are the runs with 0 . 1 Z d , ⊙ and 0 . 01 Z d , ⊙ respectively. The dashed line indicates the escape velocity for a given Σ cloud . The magnitude of the outflow velocities is significantly higher at higher Σ cloud , and smaller for lower Z d for a given Σ cloud . This potentially has implications on whether these outflows can facilitate LyC escape through the ISM of the host galaxies of the star clusters.  \n<!-- image -->  \nglyph[circledot]  \nFigure 7. Mass-weighted average radial velocity in different sightlines ( v out ), binned by the corresponding gas surface densities (Σ g ) at the characteristic times defined in Section 3.3 for Σ cloud = 10 2 M ⊙ pc -2 (top) and 10 4 M ⊙ pc -2 respectively. Dotted gray lines indicate the analytical solution for a thin shell of surface density (Σ g ) accelerated purely by radiation pressure on dusty gas (Equation 7).  \ndepends purely on age (and IMF) for the case of a single burst. This is identical to the commonly used definition in the observational literature for galaxies, without the absorption factor due to the IGM along the line of sight to the galaxy (i.e. exp( τ HI , IGM )). We can see that f esc , rel for high Σ cloud reaches higher values, and also that such clouds do so at much earlier times. In fact, for the lowest Σ cloud case, f esc , rel barely reaches values of 5% in the entirety of its evolution; the highest Σ cloud on the other hand reaches values ≳ 80% in timescales of ≲ 1 Myr. This trend is primarily because high f esc , LyC is achieved at timescales much greater (smaller) than 4 Myr for low (high) Σ cloud .  \n<!-- image -->  \n/circledot  \nThe effect of the relative timescales for the cloud to disperse and the stellar population to evolve are even more dramatic when comparing the cumulative number of LyC photons that escape the star cluster. This depends on the instantaneous emission rate of LyC photons by the cluster ( ˙ n LyC ), instantaneous LyC escape fraction ( f esc , LyC ) and the duration for which LyC escape occurs, and is given by the relation  \nN cum , LyC = ∫ t ∗ 0 ˙ n LyC f esc , LyC dt, (10)  \nwhere t ∗ is the age of the young star cluster. We show how N cum , LyC depends on Σ cloud in the right panel of Figure 9. The temporal trends are similar across Σ cloud with N cum , LyC increasing as f esc , LyC increases, followed by a saturation for t ∗ ≳ 4 Myr. This saturation occurs because the emission rate of LyC photons drops significantly beyond these times. We can see that for higher Σ cloud , the duration of being in a state of high f esc , LyC at young stellar ages ( t ∗ ≲ 4 Myr) is significantly higher. This results in total N cum , LyC over the lifetime of the young stellar population that is orders of magnitude higher than for lower Σ cloud . We can also use N cum , LyC to define a cumulative escape fraction, which is the effective fraction of LyC photons emitted by the stellar population over its lifetime that escape; i.e.  \nf esc , cum = N cum , LyC ( t life ) ∫ t life 0 ˙ n LyC dt , (11)  \nwhere we take t life = 20Myr, since the cumulative LyC emission is well saturated by this point. We show how  \nf esc , cum varies across our simulations in Figure 10. The strong trend of increasing f esc , cum with Σ cloud is evident; the effect of the dust-to-gas ratio is relatively quite minor.', '4.1. Can cloud-scale LyC leaking photons also escape the galaxy?': 'Our results imply that highly dense and compact young star clusters permit the escape of a significantly larger fraction of their intrinsic LyC emission due to their short dynamical timescales, while also emitting a larger number of LyC photons due to their higher star formation efficiencies; see Figures 9 and 10, and Table 2. In other words, the key parameters determining (cumulative) LyC escape are the dynamical timescale ( t ff ∝ n -1 / 2 cloud ) and the star formation efficiency ( ϵ ∗ ∝ Σ cloud ) at the cloud-scale. The former permits higher f esc , LyC while ˙ n LyC is high, and the latter promotes higher ˙ n LyC for a given initial M cloud . However, these findings are in the context of the cloud-scale LyC escape in an isolated cloud numerical setup. In order for LyC photons escaping clusters to contribute to cosmic reionization, the LyC photons must then still be able to escape through the turbulent ISM surrounding the dense star-forming clumps. However, high-resolution cosmological simulations have established that LyC escape in galaxies is determined primarily by escape at the scales of GMCs (Paardekooper et al. 2015; Ma et al. 2015; Trebitsch et al. 2017; Kimm et al. 2017). Therefore, our results already provide insight on the mapping between local ISM properties and galaxy escape fractions. Even so, we caution that the f esc , LyC values we report should be interpreted as upper limits when thinking about the fraction of LyC photons that escape the host galaxy, and it is the trends with cloud properties ( t ff , Σ cloud ) and the physical insights we built (i.e. bursty star formation and radiation-driven outflows) that we would highlight.  \nThere are two primary considerations in using the results presented here to keep in mind while extrapolating our cloud-scale f esc , LyC trends to obtain implications on galaxy scales: i) whether high escape fractions can still be achieved if GMCs accrete from their surroundings and ii) whether radiative feedback can also carve pathways for LyC escape from the galaxy and not just in GMCs.', '4.1.1. Accretion onto GMCs': "GMCs are highly dynamic objects that continually accrete gas from their surroundings, and simultaneously compete with feedback from star formation within them (e.g., Jeffreson et al. 2024). In the context of the re-  \nsults presented in this work, these accreting flows would act to increase the gas columns in sightlines around the young stellar population - i.e. they would affect the PDFs shown in Figure 5 - and thereby lower f esc , LyC . However, we speculate that there are likely two favorable conditions that apply, which would limit the impact of accretion. The first is the expectation that the accretion flow would be highly filamentary - akin to observed hub-filament systems in the Milky Way (Hacar et al. 2023). While these flows would be important in the context of star formation, they would only occupy a small fraction of solid angles - like blinders on a race horse - and preferentially sightlines with higher columns. Inspection of Figure 6 reveals that LyC escape is facilitated by the ejection of low-column sightlines anyways, and as long as additional accretion does not occur along those sightlines, it shouldn't affect f esc , LyC strongly. Secondly, how quickly accretion can fill sightlines would depend on the relative timescales for gas consumption due to star formation/feedback-driven ejection and the accretion rate. The smaller spatial and thereby, faster temporal, scales at which the former occur suggest that the visibility of the cluster shouldn't be strongly affected. Indeed observations in the local Universe do reveal that star clusters are visible in the optical/UV on short timescales after formation (Chevance et al. 2023) - although these may not be suitable analogs for higher density clouds that presumable lie in higher density ISM conditions. Therefore, systematic studies of the impact of accretion on cloud-scale properties, and its interaction with feedback, are required to verify these expectations.", '4.1.2. LyC escape through the ISM': 'An additional question is whether radiative feedback, which we show carves pathways for LyC escape at the cloud-scale can do so at the galactic scale. This depends on the ISM conditions outside the cloud. For a turbulent ISM, it is reasonable to assume that for the surrounding medium this can be represented as a distribution of column densities (sightlines), and some fraction of these are super-Eddington (Thompson & Krumholz 2016). These sightlines would continue to be accelerated by photoionization and radiation pressure at low and high Σ respectively. Whether these sightlines can continue to be driven to carve density-bounded channels in the galaxy depends on how fast gas is driven outward with respect to the circular/escape velocity of the galaxy, and the timescales over which this occurs with respect to the lifetimes of massive stars. We note that these two conditions make the case for the role of dense star clusters stronger. Outflow velocities are much higher for Σ cloud ≳ 10 4 M ⊙ pc -2 , ≳ 50 kms -1  \n/circledot  \nFigure 10. Cumulative LyC escape fractions ( f esc , cum ) averaged over the lifetime of the stellar population (see Equation 11) for the different simulations in our study. Filled markers are at Z d = Z d , ⊙ ; the unfilled circle and star markers indicate the 0 . 1 Z d , ⊙ and 0 . 01 Z d , ⊙ respectively.  \n<!-- image -->  \n∗  \n∗  \n∗  \nFigure 9. Left: Instantaneous value of f esc , LyC shown in terms of absolute time in Myr since star formation started ( t ∗ ) for the different simulations. Middle: f esc , rel , the LyC Escape fraction relative to the escaping bolometric FUV luminosity from a stellar population (see Equation 9). Right: The cumulative number of LyC photons that escape the cloud over the lifetime of the (instantaneous) burst (see Equation 10). The dashed gray vertical line indicates the rough timescale ( ∼ 4 Myr) when massive stars die and the LyC emission rate starts to drop. Denser clouds are able to leak LyC photons well before this timescale due to their rapid evolution governed by their much shorter t ff , cloud . This leads to significantly higher values of N cum , LyC over the lifetime of the stellar population.  \n<!-- image -->  \n/circledot  \non average (Figure 8) and even higher along sightlines with lower columns (Figure 7) - comparable in magnitude to circular velocities in compact dwarf galaxies ( M ∗ ∼ 60 kms -1 ; Simons et al. 2015). More importantly, the rapid dynamical evolution in these clouds implies that there is an important epoch where radiatively driven acceleration of super-Eddington channels can occur before the stellar population ages, the SED reddens, and supernovae take over - a primary limitation  \nin launching radiation-driven outflows of sub-regions in local star-forming galaxies (Murray et al. 2011; Blackstone & Thompson 2023). There is scope for using our insights along with semi-analytic models for radiationdriven winds in a galactic potential (e.g., Thompson et al. 2015) to quantify this possibility in future work. Idealized numerical simulations have demonstrated that radiative feedback can evacuate channels for LyC escape  \nTable 2. Summary of relevant quantities measured across our simulations.  \nNotes : Columns in order indicate - Σ cloud : mass surface density of the cloud; Z d : dust-to-gas ratio in units of the corresponding value in the solar neighborhood ( Z d , ⊙ ); ϵ ∗ ,f : final saturated value of the (integrated) star formation efficiency given by ϵ ∗ ,f = M ∗ ,f /M cloud ; t ∗ , 50 , t ∗ , 95 : timescale to form 50% and 95% of the final stellar mass respectively in Myr and in units of t ff , cloud in brackets; t esc , 0 . 2 , t esc , 0 . 2 : same, but for the times to achieve a LyC escape fraction of 20% and 70% respectively; ⟨ v out ⟩ Ω ,t : Time-averaged outflow velocity; f esc , cum : Luminosity-weighted escape fraction averaged over the lifetime of the stellar population (Equation 11).  \nin the ISM before the onset of supernova feedback (e.g., Kakiichi & Gronke 2021).  \nAn additional point to consider in the high Σ cloud regime is that these clouds likely form in starburst regions - for instance, the compact central regions of dwarf galaxies undergoing mergers or inflows, and/or massive galaxies at high redshift. This implies that there would be effects of feedback from other young star clusters simultaneously acting on the ISM in which clouds are embedded. Combined radiative feedback from clustered star formation (and possibly AGN) can drive stronger outflows and make it easier to carve optically thin channels (Fielding et al. 2018; Orr et al. 2022). Moreover, supernovae from nearby clusters would carve large-scale holes in these regions, which would imply that the only bottleneck for LyC escape is at the cloud-scale in such conditions. This supernovae-assisted high LyC escape in starburst regions has been identified in cosmological simulations (Trebitsch et al. 2017; Rosdahl et al. 2022; Yeh et al. 2023) and observations (Micheva et al. 2018; Martin et al. 2024).', '4.2. LyC leakage from dense compact star clusters: empirical evidence': 'We note that there is some (indirect) empirical evidence for the potentially important role of dense compact star clusters in LyC leakage through the scenario we discussed in this study. For instance, LyC leakers in the nearby Universe ( z ≲ 0 . 5) have been found to be quite compact ( < 1 kpc), and show signs of concentrated star formation (high Σ SFR ; Jaskot & Oey 2013; Nakajima & Ouchi 2014; Borthakur et al. 2014; Schaerer et al. 2016; Izotov et al. 2016, 2018b,a), with even signs of increasing f esc with Σ SFR (Naidu et al. 2020; Flury et al. 2022; Jaskot et al. 2024). They also seem to be characterized by nebular line emission with relatively  \nhigh contribution from high-ionization lines - such as a high [O III ]/[O II ] ratio, among others, and show steep UV spectral slopes with relatively low nebular contamination (Izotov et al. 2017, 2018b). These can be interpreted to be consistent with nebular emission along density-bounded channels around very young ( ≲ 3-4 Myr) stellar populations (Jaskot & Oey 2013; Nakajima & Ouchi 2014). The results of our study seem to suggest the correlation with these spectral features is likely directly related to the compact nature of young super star clusters in these galaxies, whose short dynamical timescales and strong radiative feedback enable LyC leakage while the cluster is still very young. Indeed, nebular signatures seem to suggest very young ages ( ∼ 2-3 Myr) and ionization parameters reminiscent of super star clusters (log U ∼ 2) in local LyC leakers (Izotov et al. 2018b) and z ∼ 2 Lymanα emitters (e.g., Naidu et al. 2022). The combination of compact star formation with weak gravitational potential wells have also been found to be characteristic of Lyα emitters at z ∼ 2 (Pucha et al. 2022). There are also tentative signs of radiatively-driven outflows from highly dense SSCs in local (low-mass) green-pea galaxies (Komarova et al. 2021).  \nMore direct evidence for localized LyC escape from pc-scale compact star formation is found in lensed highredshift galaxies, where lensing has enabled characterization of the physical properties at scales < 100 pc (Bian et al. 2017; Rivera-Thorsen et al. 2019; Vanzella et al. 2020; Kim et al. 2023b; Pascale et al. 2023; Vanzella et al. 2023; Messa et al. 2024b). For instance, in the Sunburst LyC leaker ( z = 2 . 37) with an escape fraction ≳ 50%, the LyC leakage is highly localized in a region tens of parsecs in scale that dominates the UV output of the galaxy. The inferred stellar surface densities are very  \nhigh in the LyC leaking region (Σ ∗ ∼ 10 4 M ⊙ pc -2 ) with their spectra closely resembling SSCs in nearby galaxies (Vanzella et al. 2020). In addition, this region exhibits other indirect signs of density-bounded ionized regions in this localized patch - such as the presence of high ionization nebular lines and steep UV spectral slopes (Messa et al. 2024b) - which are partly averaged out over the entire galaxy (Kim et al. 2023b). This raises the intriguing possibility that LyC leakage is highly anisotropic and dominated by localized patches, and observationally determined escape fractions - when the galaxy is unresolved - are highly dependent on whether these privileged sightlines are accessible in the observations (Nakajima et al. 2020).', '4.3. Implications for reionization by star-forming galaxies': "Our results imply that galaxies that can form a significant population of young, compact, dense star clusters could be proficient LyC leakers. The sensitive dependence on the density of the star cluster (Figure 10) for attaining high (localized) f esc provides a natural explanation for why f esc is typically low in large samples of star-forming galaxies, especially in the local Universe. For instance, stacked samples of intermediate redshift ( z ∼ 2 -4) star-forming galaxies provide upper limits of f esc ≲ 10% (Japelj et al. 2017; Wang et al. 2023; Jung et al. 2024) - much lower than the average f esc ∼ 20% required to match reionization timing constraints (e.g., Robertson et al. 2013). However, there are also a small sample of galaxies that show highly elevated f esc ≳ 40% (Vanzella et al. 2012; Shapley et al. 2016; Izotov et al. 2018b; Rivera-Thorsen et al. 2019; Vanzella et al. 2022; Wang et al. 2023; Marques-Chaves et al. 2024). Apart from line-of-sight anisotropies, it is possible that these galaxies have undergone a recent surge in the formation of these highly dense compact star clusters in high ISM-pressure environments, which may be facilitated via mergers and/or the triggering of central starbursts. The Sunburst LyC leaker does indeed show the presence of a recent bursty event that likely formed the population of young bound star clusters that show evidence of localized LyC leakage (Vanzella et al. 2022; Kim et al. 2023b). This is also consistent with the findings from cosmological simulations on a correlation of high f esc with recent ( ∼ few Myr) bursts of star formation (Rosdahl et al. 2022). Therefore, the low f esc in stacks might simply resemble the small fraction of galaxies that happen to have undergone a recent surge of efficient localized bursts of star formation. Along similar lines, the very low LyC escape fractions in star-forming galaxies of the local Universe may simply reflect the rel-  \nwer ISM-pressures ( P ∼ 10 4 Kcm -3 ) encountered, which are not conducive to the formation of dense stellar systems.  \nIn the early Universe ( z > 6), cosmological evolution of the mean baryon density and sizes ensure that galaxies are compact and have high average densities and short dynamical timescales (Faucher-Gigu'ere 2018). Indeed, JWST has revealed that galaxies at this epoch tend to be compact ( < 0 . 5 pc) and can achieve stellar surface densities Σ ∗ ≳ 10 4 M ⊙ pc -2 over the whole galaxy (Baggen et al. 2023; Casey et al. 2023; Morishita et al. 2023). Such high densities would enable the prodigious formation of several compact clusters with Σ ∗ (at least) of this order, suggesting that this extreme mode of star formation might be relatively commonplace at these redshifts. This is consistent with semianalytic models that have empirically constrained that most ( ∼ 50 -60%) star formation at this epoch would be distributed in massive bound clusters (Belokurov & Kravtsov 2023). Lensed fields have also revealed several young star clusters that achieve these high Σ ∗ values at these redshifts (Adamo et al. 2024; Mowla et al. 2024; Messa et al. 2024a; Fujimoto et al. 2024; Messa et al. 2024b) consistent with these expectations. If this is true - and taking our results at face value - we would expect that LyC escape fractions from these compact galaxies would be systematically higher, and this would be the case when the stellar populations are still young. JWST has provided clues that this may be occurring with the detection of several galaxies that are characterized by high ionizing efficiencies indicative of very young ages and low dust attenuation (e.g., Atek et al. 2024; Simmonds et al. 2024), suggesting that the stellar populations have very quickly cleared out the surrounding gas. In addition, Lymanα emitting galaxies at extremely high redshift ( z > 10) seem to be characterized by very high f esc ∼ 80% (Curtis-Lake et al. 2023; Tacchella et al. 2023; Hainline et al. 2024), consistent with expectations from our results when considering that these galaxies would have their star formation primarily distibuted in highly compact dense star clusters (e.g., Dekel et al. 2023). These considerations make the case that star-forming galaxies are likely more than capable of reionizing the Universe on their own. In fact, the high ionizing production efficiencies and escape fractions may possibly even lead to a 'photon budget crisis' wherein the Universe is reionized too early (Mu˜noz et al. 2024). There is scope to explore these implications by adopting our cloud-scale escape fractions as subgrid prescriptions in semi-analytic galaxy formation models and/or numerical simulations of cosmological volumes.", '4.4. Comparison to other numerical simulations': 'There have been earlier numerical simulations in the literature that made predictions for the escape fraction and its dependence on cloud properties. A comparison to these is not straightforward in practice, since the numerical setups, radiative feedback implementations, physics included and subgrid sink particle algorithms vary amongst the simulations. Therefore we do not attempt to make a one-to-one comparison of escape fractions obtained across the simulations. However, it is still useful to assess whether there is agreement with the trends of escape fractions with cloud properties. For instance, our study suggests that the key physical parameters determining the cumulative LyC escape are t ff and ϵ ∗ , which depend on n cloud and Σ cloud respectively. This is in agreement with the findings of Dale et al. (2012) who found that dense, compact clouds with short t ff ≲ 3 Myr have f esc , cum > 0 . 8, whereas more massive clouds (for a given size) with longer t ff have much lower f esc , cum . Kim et al. (2019) also find a trend of increasing f esc , cum for more compact clouds with shorter t ff . However, we note that Figure 6 in their paper seems to indicate that f esc , cum decreases with Σ cloud for a given cloud mass; we stress that this trend only holds for Σ cloud ≲ 100 M ⊙ pc -2 which is lower than the range we probe; for Σ cloud ≳ 100 M ⊙ pc -2 , f esc , cum shows signs of increasing with Σ cloud consistent with our findings. He et al. (2020) also report a clear trend of increasing f esc , cum with cloud density; they also show that the total number of ionizing photons escaping the cloud (per unit mass) increases with M cloud for compact clouds due to their higher Σ cloud (as ϵ ∗ is higher). On the other hand, Kimm et al. (2022) report that the escape fractions for denser and more efficiently star-forming clouds are lower than diffuse clouds - which is in disagreement with the trend we presented in this work. They attribute their trend to the longer timescales for photoevaporating the cloud and opening up optically-thin channels at higher Σ cloud . While we too find that higher Σ cloud requires more dynamical timescales for LyC photons to start escaping (Figure 4), the factor by which this timescale is longer is smaller than the factor by which the dynamical timescales are shorter at higher Σ cloud , unlike what Kimm et al. (2022) find, indicating that we need to identify possible reasons for these differences. Before doing so, we note that in Kimm et al. (2022), the difference in escape fraction between their denser and diffuse cloud 6 - which vary in t ff by a factor of ∼ 2 - is only  \n10%, and therefore shows a relatively weak difference. Regardless, we note that differences can arise due to different treatments for feedback in their study. For instance, Kimm et al. (2022) assume that most dust (99%) is destroyed in H II regions; this would lower the effect of radiation pressure on dust which we find are important to expand the H II regions in denser clouds, and drive their subsequent photoevaporation. Conversely, for the more diffuse cloud where thermal pressure dominates over radiation pressure, lower dust abundance in ionized gas can permit LyC photons to photoionize material at larger radii. In other words, the impact of destroying dust in ionized gas could work in opposite directions in the regimes where photoionization and radiation pressure dominates, where it aids and hampers LyC escape respectively. Kimm et al. (2022) also include supernovae, which can help disrupt the more diffuse cloud with longer t ff earlier. In addition, the M 1 numerical implementation for radiative feedback in their study, and associated flux cancellations in optically-thin media intrinsic to the approach (see, for instance, Figure 13 in Menon et al. 2022b), could also lower the ability for radiative feedback to evacuate its surroundings. That being said, their higher dynamic range in numerical resolution and star-by-star Poisson sampling approach could also capture more localized effects of LyC absorption although He et al. (2020), who also have these advantages agree with our trends. Any of these factors, or a combination thereof could likely explain the differences in trends, motivating a systematic comparison between different numerical setups and physics implementations and their corresponding impacts on LyC escape.', '4.5. Caveats': 'We have made several choices in our modelling approach, some of which may have implications for our results; we briefly discuss some of them here.', '4.5.1. Binary evolution and radiative outputs': 'In our simulations, we have used the LyC and FUV radiative outputs for stellar populations based on the SLUG stellar population synthesis code, which uses the MIST stellar evolution models (Choi et al. 2016). This model does not account for the presence of massive star binaries and associated physics, or the higher stellar surface temperatures expected at lower metallicities, both of which can increase the intrinsic LyC emission and the timescales up to which LyC photons are emitted. This could possibly have implications on the relative contribution of Σ cloud ≲ 10 3 M ⊙ pc -2 clouds to LyC escape, since we find that the primary reason for their lower f esc , cum is that their intrinsic emission drops by the time they disperse their clouds - a longer timescale for LyC  \nproduction might reduce this effect. Indeed Ma et al. (2016) find that adopting stellar models that include the effects of binary evolution can significantly boost the luminosity-averaged escape fractions in cosmological simulations. In Appendix A we check whether adopting the Binary Population and Spectral Synthesis ( BPASS ) model v2.2 (Stanway & Eldridge 2018) would affect our results on the relative importance of dense clouds on LyC escape. We find that adopting BPASS models can indeed increase N cum , LyC for Σ cloud = 100 M ⊙ pc -2 by factors of up to 10 and decrease the gap in the lifetimeaveraged N cum , LyC with denser clouds. However, this is still insufficient to catch up to the ≳ 100 higher N cum , LyC for Σ cloud ≳ 10 4 M ⊙ pc -2 , providing reassurance that their disproportionate degree of LyC leakage still holds.', '4.5.2. Dust in ionized gas': 'In our simulations we make the assumption that dust grains survive in photoionized gas, and the dust opacity does not change when the gas becomes ionized. However, this is unlikely to be true at UV wavelengths since it is the presence of small carbanaceous grains and polycyclic aromatic hydrocarbons (PAHs) that dominate the extinction at these wavelengths, and these types of grains are likely destroyed by intense UV radiation (Draine & Li 2001) - findings with JWST in H II regions of local galaxies are consistent with this picture (e.g., Egorov et al. 2023). This could impact the absorption coefficients at UV wavelengths, which would lower the amount of dust absorption in photoionized gas, and thereby increase escape fractions. However, the lack of dust absorption in photoionized gas would also reduce the impact of radiation pressure. This is unlikely to affect the evolution of lower Σ cloud cases, since photoionization is the dominant radiative feedback mechanism anyways. This has been shown by Kim et al. (2019) where they considered the two extremes of complete dust survival and destruction for a cloud with Σ cloud ∼ 80 M ⊙ pc -2 , and found that the broad evolution is highly similar; the escape fractions are much higher due to lowered dust extinctions. However, for the denser clouds, it is possible that dust destruction would affect the extent to which radiation pressure can impact the evolution of the clouds. For instance, the simulations with lower dust abundances in our study do have lower outflow velocities, hinting at the diminished impact of radiation pressure. However, this is uncertain since the increased transparency of LyC photons through photoionized gas (as opposed to everywhere) might also shift the momentum deposition due to radiation pressure to larger scales, and make it easier to accelerate gas. Future simulations that systematically  \nquantify the impact of PAH destruction on dust UV opacities, and explore the impact of the modified opacities on dense compact clouds will help to clarify these issues.', '4.5.3. Missing feedback mechanisms': 'We only model early radiative feedback in our simulations, and do not include protostellar jets, stellar winds and supernovae. Jets are unlikely to affect the bulk evolution of our clouds due to their limited momentum injection when averaged over a stellar population; for instance Matzner & Jumper (2015) show that jets can unbind clouds with v esc ≲ 1 kms -1 , reducing their effect to localized clumps 7 . Stellar winds on the other hand have a momentum contribution of the same order of magnitude as (UV) radiation pressure. Therefore, we expect their inclusion to increase escape fractions by making it easier to disperse the gas. However, the relative importance of winds as a feedback mechanism has been shown to be more prominent for denser clouds (Lancaster et al. 2021), which suggests that their inclusion would only amplify the trends with Σ cloud we report in this study. On the other hand, the omission of supernovae would primarily affect the escape fractions of the lowest density cloud (i.e. Σ cloud = 100M ⊙ pc -2 ), since Figure 9 show that these clouds start leaking LyC photons at times comparable to the lifetime of massive stars ( ∼ 4 Myr). Supernovae could quickly remove the surrounding gas when they explode, possibly achieving higher f esc , LyC slightly earlier than we find in our simulations. However, we emphasize that this should not affect our conclusions significantly since the onset of supernovae is precisely the point where the intrinsic LyC emission starts to drop. Therefore the luminosityaveraged escape fraction ( f esc , cum ) and its trends with Σ cloud (Figure 10) is unlikely to be strongly affected.', '5. SUMMARY': "We conduct idealized radiation hydrodynamic simulations of self-gravitating turbulent clouds with different initial cloud gas surface densities (Σ cloud ) and dust abundances ( Z d ) to systematically quantify the escape fraction of LyC (and FUV) photons and its dependence on cloud properties. We use the FLASH MHD code with VETTAM - a highly accurate variable Eddington tensorclosed moment method - to model the transport of LyC,  \nFUV and IR photon bands, and the associated ionization and radiation pressure forces. We explore clouds with Σ cloud ∼ 10 2 -10 5 M ⊙ pc -2 , which resembles conditions found in star cluster-forming clumps typical of star-forming galaxies in the local Universe, to the gasrich, high ISM-pressure conditions that are typical of the epoch of reionization (EoR). We also experiment with Z d = 0 . 01 , 0 . 1 , 1 Z d , ⊙ where Z d , ⊙ is the dust abundance in the solar neighborhood (see Table 1). We summarize relevant quantities measured in our study in Table 2, and draw the following conclusions.  \n- 1. Higher Σ cloud results in much higher (integrated) star formation efficiencies ( ϵ ∗ ) with ϵ ∗ ∼ 15% and 95% for the Σ cloud = 10 2 and 10 5 M ⊙ pc -2 respectively (Figure 1).\n- 2. f esc , LyC increases with time as star formation proceeds, and the radiative feedback increasingly disperses the cloud (Figure 2 and 3). At higher Σ cloud , feedback takes a larger number of dynamical times to start allowing the escape of LyC photons due to higher optical depths (Figure 4).\n- 3. Neutral gas absorbs most LyC photons at early times, and f esc , LyC is high at late times; only for a brief period (when f esc , LyC ∼ 50%) does dust dominate the absorption, and that too only at Z d = Z d , ⊙ . However, this period is longer for higher Σ cloud (Figure 2).\n- 4. The longer duration (in dynamical times) of efficient LyC escape at higher Σ cloud can be explained by the competition between UV radiation pressure and gravity (Eddington limit) and the larger fraction of gas that exceeds this limit at higher Σ cloud (Figures 5 and 6).\n- 5. Radiative feedback drives super-Eddington sightlines in an ionized outflow, producing densitybounded channels for LyC escape, with faster outflow velocities at lower columns (Figure 7), and average velocities that are higher for denser clouds (Figure 8).\n- 6. While denser clouds (Σ cloud ≳ 10 4 M ⊙ pc -2 ) take a larger number of dynamical times to permit LyC escape, their significantly shorter dynamical times cause this trend to reverse when viewed in terms of the stellar population age in absolute time (Figure 9), with f esc , LyC > 50% in ≲ 2 Myr. This ensures that the phases of high f esc , LyC and high intrinsic LyC emission from stellar populations ( ≲ 3 Myr) are occurring contemporaneously, such  \nthat the cumulative number of LyC photons that escape are significantly higher.  \n- 7. This results in a luminosity-weighted average escape fraction f esc , cum that is much higher for denser clouds: ≳ 70% and ≲ 10% for Σ cloud ≳ 10 4 M ⊙ pc -2 and Σ cloud ∼ 10 2 M ⊙ pc -2 respectively (Figure 10).  \nWe provide a brief discussion on what implications the cloud-scale LyC escape would have on its counterpart at the galaxy scale (Section 4.1), compile some signs of possible evidence in observational studies of the aforementioned scenario playing out (Section 4.2), and discuss its implications on cosmic reionization (Section 4.3). We compare with results of other numerical simulations in the literature and discuss agreements/disagreements with them (Section 4.4) and discuss some caveats associated with our simulations and their possible implications on our findings (Section 4.5). Overall, our results seem to suggest a scenario where dense compact bursts of star formation forming bound clusters (Σ > 10 4 M ⊙ pc -2 ) in a galaxy are highly efficient at leaking LyC photons due to their rapid bursts of star formation and radiationdriven outflows, which carve density-bounded channels before massive stars in the cluster die and the intrinsic LyC emission drops. Increasing evidence for clumpy star clusters in (lensed) galaxies at the EoR, approaching (and exceeding) these 'extreme' stellar surface densities (Vanzella et al. 2023; Adamo et al. 2024; Messa et al. 2024a; Mowla et al. 2024; Fujimoto et al. 2024; Messa et al. 2024b) hints at their important role as the drivers of reionization of the Universe.  \nS. H. M would like to thank Nicholas Choustikov for extensive discussions on escape fractions in galaxies. S. H. M also thanks Eve Ostriker, Rahul Kannan, Aaron Smith, Ulrich Steinwandel, Lachlan Lancaster, Oleg Gnedin and Crystal Martin for useful discussions in the context of this paper. We acknowledge high-performance computing resources provided by the provided by the Simons Foundation as part of the CCA. B.B. and S.H.M. acknowledge the support of NASA grant No. 80NSSC20K0500 and NSF grant AST2009679. B.B. also thanks the Alfred P. Sloan Foundation and the Packard Foundation for support. This research was supported in part by grant NSF PHY2309135 to the Kavli Institute for Theoretical Physics (KITP). T.A.T. is supported in part by NASA grant No. 80NSSC23K1480. This research has made use of NASA's Astrophysics Data System (ADS) Bibliographic 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17  \n- Services. 18  \nSoftware: SLUG (da Silva et al. 2012; Krumholz et al. 2015), PETSc (Balay et al. 1997, 2021), NumPy (Harris et al. 2020), SciPy (Virtanen et al. 2020), Matplotlib  \n(Hunter 2007), yt (Turk et al. 2010), TurbGen (Federrath et al. 2010b, 2022), CFpack .", 'A. EFFECT OF BINARY STELLAR EVOLUTION': 'To check whether binary evolution/metallicity would have an impact on our conclusions, we (re)compute and compare the time evolution of N cum , LyC (Equation 10) for Σ cloud = 100M ⊙ pc -2 and 10 4 M ⊙ pc -2 , but with the intrinsic emission ( ˙ n ion ) taken from the BPASS model v2.2 (Stanway & Eldridge 2018). We use the (pre-computed) LyC emission rates for a 10 6 M ⊙ stellar population with a Chabrier (2003) IMF, an upper mass limit of 100M ⊙ , and stellar metallcities of 0 . 1 Z ⊙ and 0 . 01 Z ⊙ . We use this to compare the LyC emission rate per unit-mass (Ψ LyC ) in these models with the adopted MIST v1.0 values (Choi et al. 2016) in our simulations; this is shown in Figure 11. We can clearly see the enhanced LyC outputs at young ages and more importantly, the longer timescales of LyC emission in the BPASS models. We use these to (re)compute the time evolution of N cum , LyC . We do not rerun the simulations, so f esc , LyC ( t ) and M ∗ ( t ) are identical to those presented in the paper. While the BPASS luminosities would also change this in principle - permitting quicker cloud disruption due to enhanced UV outputs - we do not expect it to change significantly since the enhancement of the UV output is only a factor ≲ 3. Moreover, this elevated UV outputs would also slightly decrease the star formation efficiency (by ≲ 20%; Menon et al. 2024), so the combined effect on N cum , LyC remains likely small. As for the longer timescale of intrinsic LyC emission, we did not account for the temporal drop in intrinsic UV emission in our simulations anyways, so are already overestimating the impact of radiative feedback at late times ( t ≳ 4 Myr) at low Σ cloud and therefore f esc , LyC . These considerations suggest that sticking with the f esc , LyC ( t ) we have obtained is reasonable.  \nIn Figure 12 we compare the time evolution of N cum , LyC obtained with the BPASS and MIST models for Σ cloud = 100 M ⊙ pc -2 and 10 4 M ⊙ pc -2 to see its effect on their relative (leaking) LyC outputs. We can see that the impact of BPASS models is more pronounced for the lower Σ cloud case - increasing N cum , LyC by a factor ∼ 10; this can be easily understood as arising from the longer intrinsic LyC emission in these models. The impact on the denser cloud is milder since the timescale of LyC leakage is very short anyways, and the additional impact at longer timescales is modest. Regardless, even with the factor ∼ 10 increase, N cum , LyC for the denser cloud is still higher by a factor ∼ 10. Therefore, adopting BPASS models do not qualitatively change the picture of disproportionate LyC leakage from dense clusters, but only quantitatively change the degree to which the denser clouds are important.', 'REFERENCES': "Adamo, A., Bradley, L. D., Vanzella, E., et al. 2024, arXiv e-prints, arXiv:2401.03224,  \ndoi: 10.48550/arXiv.2401.03224  \nAndrews, B. H., & Thompson, T. A. 2011, ApJ, 727, 97,  \ndoi: 10.1088/0004-637X/727/2/97  \nAtek, H., Labb'e, I., Furtak, L. J., et al. 2024, Nature, 626, 975, doi: 10.1038/s41586-024-07043-6  \nBaggen, J. F. W., van Dokkum, P., Labb'e, I., et al. 2023, ApJL, 955, L12, doi: 10.3847/2041-8213/acf5ef  \nBalay, S., Gropp, W. D., McInnes, L. C., & Smith, B. F. 1997, in Modern Software Tools in Scientific Computing, ed. E. Arge, A. M. Bruaset, & H. P. Langtangen (Birkhauser Press), 163-202  \nBalay, S., Abhyankar, S., Adams, M. F., et al. 2021, PETSc/TAO Users Manual, Tech. Rep. ANL-21/39 Revision 3.16, Argonne National Laboratory  \nBelokurov, V., & Kravtsov, A. 2023, MNRAS, 525, 4456, doi: 10.1093/mnras/stad2241 Bera, A., Hassan, S., Smith, A., et al. 2023, ApJ, 959, 2, doi: 10.3847/1538-4357/ad05c0 Bian, F., Fan, X., McGreer, I., Cai, Z., & Jiang, L. 2017, ApJL, 837, L12, doi: 10.3847/2041-8213/aa5ff7 Blackstone, I., & Thompson, T. A. 2023, arXiv e-prints, arXiv:2302.10136, doi: 10.48550/arXiv.2302.10136 Borthakur, S., Heckman, T. M., Leitherer, C., & Overzier, R. A. 2014, Science, 346, 216, doi: 10.1126/science.1254214 Bosman, S. E. I., Davies, F. B., Becker, G. D., et al. 2022, MNRAS, 514, 55, doi: 10.1093/mnras/stac1046 Buntemeyer, L., Banerjee, R., Peters, T., Klassen, M., & Pudritz, R. E. 2016, New Astronomy, 43, 49, doi: 10.1016/j.newast.2015.07.002  \nFigure 11. LyC photon emission rate per unit mass (Ψ LyC ) emitted by (fully IMF-sampled) stellar populations using the ( MIST ) v1.0 isochrones (Choi et al. 2016) used in SLUG - which we adopt in this study - compared with that obtained with the BPASS v2.2 binary stellar population models (Stanway & Eldridge 2018) at metallicities of 0 . 1 Z ⊙ and 0 . 01 Z ⊙ . We can see that Ψ LyC is higher and remains higher for longer in the BPASS models.  \n<!-- image -->  \n∗  \nBurkhart, B., Tillman, M., Gurvich, A. B., et al. 2022, ApJL, 933, L46, doi: 10.3847/2041-8213/ac7e49 Casey, C. M., Akins, H. B., Shuntov, M., et al. 2023, arXiv e-prints, arXiv:2308.10932, doi: 10.48550/arXiv.2308.10932 Chabrier, G. 2003, PASP, 115, 763, doi: 10.1086/376392 Chevance, M., Krumholz, M. R., McLeod, A. F., et al. 2023, in Astronomical Society of the Pacific Conference Series, Vol. 534, Protostars and Planets VII, ed. S. Inutsuka, Y. Aikawa, T. Muto, K. Tomida, & M. Tamura, 1, doi: 10.48550/arXiv.2203.09570 Choi, J., Dotter, A., Conroy, C., et al. 2016, ApJ, 823, 102, doi: 10.3847/0004-637X/823/2/102 Choi, Y., Dalcanton, J. J., Williams, B. F., et al. 2020, ApJ, 902, 54, doi: 10.3847/1538-4357/abb467 Choustikov, N., Katz, H., Saxena, A., et al. 2024, MNRAS, 529, 3751, doi: 10.1093/mnras/stae776  \nCrocker, R. M., Krumholz, M. R., Thompson, T. A., Baumgardt, H., & Mackey, D. 2018, MNRAS, 481, 4895, doi: 10.1093/mnras/sty2659  \nCurti, M., Maiolino, R., Curtis-Lake, E., et al. 2023, arXiv e-prints, arXiv:2304.08516,  \ndoi: 10.48550/arXiv.2304.08516  \nCurtis-Lake, E., Carniani, S., Cameron, A., et al. 2023, Nature Astronomy, 7, 622,  \ndoi: 10.1038/s41550-023-01918-w  \nda Silva, R. L., Fumagalli, M., & Krumholz, M. 2012, ApJ,  \n745, 145, doi: 10.1088/0004-637X/745/2/145  \nDale, J. E., Ercolano, B., & Bonnell, I. A. 2012, MNRAS, 424, 377, doi: 10.1111/j.1365-2966.2012.21205.x  \nDayal, P., & Ferrara, A. 2018, PhR, 780, 1,  \ndoi: 10.1016/j.physrep.2018.10.002  \nDayal, P., Volonteri, M., Greene, J. E., et al. 2024, arXiv e-prints, arXiv:2401.11242,  \ndoi: 10.48550/arXiv.2401.11242  \nDekel, A., Sarkar, K. C., Birnboim, Y., Mandelker, N., & Li, Z. 2023, MNRAS, 523, 3201,  \ndoi: 10.1093/mnras/stad1557  \nDraine, B. T. 2011, ApJ, 732, 100,  \ndoi: 10.1088/0004-637X/732/2/100  \n/circledot  \nFigure 12. Comparison of the cumulative number of escaping LyC photons ( N cum , LyC ; Equation 10) obtained with our adopted MIST stellar tracks for the LyC emission rate, with the corresponding curves obtained with the BPASS models. We show only two cases, Σ cloud = 100 M ⊙ pc -2 and 10 4 M ⊙ pc -2 , to check if the trends with Σ cloud reported in our paper are robust to these changes. To produce these curves we still use f esc , LyC obtained in our simulations (which use MIST ), and only use the BPASS for ˙ n LyC in Equation 10. We can see that the relative effect of higher and more extended Ψ LyC in the BPASS models is more prominent at lower Σ cloud , yet is insufficient to strongly affect the trend of N cum , LyC with Σ cloud that we have reported in this work.  \n<!-- image -->  \n∗  \nDraine, B. T., & Li, A. 2001, ApJ, 551, 807,  \ndoi: 10.1086/320227  \nDubey, A., Reid, L. B., & Fisher, R. 2008, Physica Scripta Volume T, 132, 014046,  \ndoi: 10.1088/0031-8949/2008/T132/014046  \nEgorov, O. V., Kreckel, K., Sandstrom, K. M., et al. 2023, ApJL, 944, L16, doi: 10.3847/2041-8213/acac92  \nElmegreen, B. G., & Chiang, W. H. 1982, ApJ, 253, 666, doi: 10.1086/159667  \nFall, S. M., Krumholz, M. R., & Matzner, C. D. 2010, Astrophysical Journal Letters, 710, L142,  \ndoi: 10.1088/2041-8205/710/2/L142  \nFan, X., Carilli, C. L., & Keating, B. 2006, ARA&A, 44, 415, doi: 10.1146/annurev.astro.44.051905.092514 Faucher-Gigu'ere, C.-A. 2018, MNRAS, 473, 3717,  \ndoi: 10.1093/mnras/stx2595  \nFederrath, C., Banerjee, R., Clark, P. C., & Klessen, R. S. 2010a, ApJ, 713, 269, doi: 10.1088/0004-637X/713/1/269  \nFederrath, C., Banerjee, R., Seifried, D., Clark, P. C., & Klessen, R. S. 2011, in IAU Symposium, Vol. 270, Computational Star Formation, ed. J. Alves, B. G. Elmegreen, J. M. Girart, & V. Trimble, 425-428, doi: 10.1017/S1743921311000755 Federrath, C., Roman-Duval, J., Klessen, R. S., Schmidt, W., & Mac Low, M. M. 2010b, A&A, 512, A81, doi: 10.1051/0004-6361/200912437 -. 2022, TG: Turbulence Generator, Astrophysics Source Code Library, record ascl:2204.001. http://ascl.net/2204.001 Fielding, D., Quataert, E., & Martizzi, D. 2018, MNRAS, 481, 3325, doi: 10.1093/mnras/sty2466 Finkelstein, S. L., D'Aloisio, A., Paardekooper, J.-P., et al. 2019, ApJ, 879, 36, doi: 10.3847/1538-4357/ab1ea8 Finkelstein, S. L., Bagley, M. B., Ferguson, H. C., et al. 2023, ApJL, 946, L13, doi: 10.3847/2041-8213/acade4 Flury, S. R., Jaskot, A. E., Ferguson, H. C., et al. 2022, ApJ, 930, 126, doi: 10.3847/1538-4357/ac61e4"}