Impact of Radiation Feedback on the Formation of Globular Cluster Candidates during Cloud-Cloud Collisions

To understand the impact of radiation feedback during the formation of a globular cluster (GC), we simulate a head-on collision of two turbulent giant molecular clouds (GMCs). A series of idealized radiation-hydrodynamic simulations is performed, with and without stellar radiation or Type II supernovae. We find that a gravitationally bound, compact star cluster of mass $M_{\rm GC} \sim 10^5\,M_\odot$ forms within $\approx 3\,{\rm Myr}$ when two GMCs with mass $M_{\rm GMC}=3.6\times10^5\,M_\odot$ collide. The GC candidate does not form during a single collapsing event but emerges due to the mergers of local dense gas clumps and gas accretion. The momentum transfer due to the absorption of the ionizing radiation is the dominant feedback process that suppresses the gas collapse and photoionization becomes efficient once a sufficient number of stars form. The cluster mass is larger by a factor of $\sim 2$ when the radiation feedback is neglected, and the difference is slightly more pronounced (16\%) when extreme $\rm Ly\alpha$ feedback is considered in the fiducial run. In the simulations with radiation feedback, supernovae explode after the star-forming clouds are dispersed, and their metal ejecta are not instantaneously recycled to form stars.


INTRODUCTION
The formation of a globular cluster (GC) has been a long-standing puzzle in galactic astronomy. Peebles (1984) first put forward the idea that the primordial gas in a dark matter halo may rapidly condense into a dense structure due to Lyman alpha cooling. Cosmological radiation-hydrodynamic (RHD) simulations verified that halos with mass M h ∼ 10 8 M can form a dense cluster when purely atomic gas collapses during the transition from a molecular-cooling to atomiccooling regime (Kimm et al. 2016); however, the predicted internal metallicity spread in each GC candidate was significantly larger than that observed. Moreover, Trenti et al. (2015) proposed a similar scenario wherein the mergers of star-free mini-halos cause shock-induced gas collapse but the halo gas needs to be somehow predaniel.han@yonsei.ac.kr tkimm@yonsei.ac.kr enriched by external sources to match the observed GC metallicity. Alternatively, a rapid collapse may occur inside a massive giant molecular cloud (GMC) in starforming galaxies at high redshift (Kravtsov & Gnedin 2005;Howard et al. 2019). This idea forms the basis of tagging methods wherein a fraction of star particles are assumed to represent potential GC candidates (e.g. Renaud et al. 2017;Pfeffer et al. 2018). Despite their simplicity, the models can successfully reproduce the color bimodality in GCs, which makes the formation channel appealing. The most relevant observational evidence in support of this scenario stems from the detection of young super star clusters in the Antennae galaxies (Whitmore & Schweizer 1995). Additionally, compressive tides and/or shock compression due to galaxy mergers are likely to have played an important role in GC formation (Renaud et al. 2015;Kim et al. 2018b;Lahén et al. 2019). Improving upon tagging methods, Grudić et al. (2022b) applied a cluster formation model derived from their GMC simulations (Grudić et al. 2021) to a Milky Way-like galaxy, and showed that the clus-ter mass-size relation observed in nearby star-forming galaxies can be explained if denser GMCs form clusters with high efficiency. However, the detailed formation processes have not yet been rigorously investigated.
In this paper, we focus on another scenario in which GCs form through the collision of GMCs. GMCs are known to be birthplaces of present-day star clusters (e.g. Cambrésy et al. 2002;Heyer et al. 2006). Although the star clusters tend to be less massive and compact than GCs, observational analyses of DR 21 (Dobashi et al. 2019), R136 (Fukui et al. 2017), and the Antennae galaxies (Tsuge et al. 2021) suggest that cloud collisions likely trigger star formation. Indeed, previous studies based on hydrodynamic simulations demonstrated that dense cores form in converging flows (Heitsch et al. 2009) and perhaps more effectively via cloud-cloud collisions (e.g., Takahira et al. 2018;Sakre et al. 2021). Furthermore, Tanvir & Dale (2020) showed that considerably more stars are created when colliding clouds are gravitationally bound, than when bound-unbound or unboundunbound clouds collide. Moreover, Dobbs et al. (2020) demonstrated that for the formation of young massive clusters, high gas densities (≥ 100 cm −3 ) and velocities greater than 20 km s −1 are required when two supervirial clouds collide.
Studies have also shown that the effects of radiation need to be understood to comprehend the evolution of GMCs. Krumholz et al. (2007) performed AU-scale simulations of turbulent protostellar cores with mass of ∼ 100 M and showed that the initial stellar mass function is significantly affected by the heating due to photoionization. Additionally, Geen et al. (2017) argued that the observed star formation efficiency of ∼ 10 % in the Milky Way molecular clouds could be reproduced by including ionizing radiation in magnetized GMC simulations. Kim et al. (2018a) showed that photoionization feedback from stars effectively controls star formation efficiency by disrupting the GMCs within ∼ 5 Myr, as supported by the CO-to-Hα analysis of GMCs in starforming galaxies (e.g., Chevance et al. 2022). In contrast, stellar winds, which are another form of early feedback, seemingly play a minor role iRn setting the properties of star clusters (Dale et al. 2013;Grudić et al. 2021).
Recently, Dobbs et al. (2022) reported the impact of photoionization in a dense region of a spiral galaxy with converging flows. More specifically, they simulated cluster formation in strongly and moderately converging regions along a spiral arm in the disks of galaxies similar to the Milky Way. They found that the star formation rate decrease by ∼ 20% in the isolated filamentary structures when radiation feedback was included, whereas the ef-fect was negligible in regions with multiple converging filaments.
Building upon previous studies, this work aims to investigate the impact of radiation feedback on the formation process of a massive compact cluster during cloudcloud collisions, which may be seen as a progenitor of present-day GCs. Although such collisions between massive clouds seem to be rare in the local Universe, these may occur more frequently at high redshifts during the formation of dwarf galaxies (Kim et al. 2018b;Sameie et al. 2022). Thus, we employ RHD simulations with various forms of feedback, including photoionization, direct radiation pressure, non-thermal pressure due to multiple scattering of IR and Lyα photons, and Type II SN explosions.
This paper is organized as follows. Section 2 describes the initial conditions, numerical methods, and physical parameters of our simulations. Section 3 examines the formation process of a massive compact cluster and effects of stellar feedback on it. Section 4 discusses the relative importance of several forms of stellar feedback, the compactness of the massive cluster, and caveats in our simulations. Finally, Section 5 summarizes our results.

RADIATION-HYDRODYNAMIC SIMULATIONS
We perform RHD simulations of cloud-cloud collisions with the adaptive mesh refinement code ramsesrt (Teyssier 2002;Rosdahl & Teyssier 2015). For each spherical cloud, an initial number density of n H = 100 cm −3 and temperature of T = 100 K are adopted. The initial radius of the clouds is set as 30 pc, based on the mass-size relations of observed GMCs Heyer et al. 2009;Oka et al. 2001). We place two clouds with mass M GMC = 3.6 × 10 5 M 100 pc apart in a uniform medium with density n H = 0.54 cm −3 and temperature 10 4 K. Turbulent velocity fields following the Kolmogorov power spectrum are developed over half of the free-fall time (0.5 t ff ≈ 2.2 Myr) of the cloud, at which point the virial parameter of each cloud reaches α vir = 0.7. We generate the random force fields via an Ornstein-Uhlenbeck process in a Fourier domain with a compressive to solenoidal ratio of 1:3 (see e.g., Federrath et al. 2010). To avoid a bulk cloud-scale motion due to the random forces, we limit the wave numbers to be greater than 6. The simulated volume of (512 pc) 3 is covered by 128 3 coarse cells with 4 pc resolution, which are further refined up to 0.25 pc (fiducial) or 0.125 pc (high-resolution case) when the gas mass in a cell exceeds 0.1 M . To better capture turbulent structures, the thermal Jeans length is resolved using at least 32 cells until the maximum refinement level Table 1. Input parameters and physical ingredients of the simulations. From left to right, each column indicates the initial gas mass of a cloud (MGMC), minimum cell size (∆xmin), mass resolution of a star particle (m res star ), inclusion of radiation feedback (RF), SN, and Lyα pressure, time at the end of simulation (t final ), total mass of star particles (Mstar), the half-mass radius (R c,h ), and stellar mass of the GC candidate (MGC). The length of the simulated box is 512 pc. The simulations with RF include photoionization, direct radiation pressure, and non-thermal pressure due to the multiple scattering of IR photons. Run-woRF-woSN 3.6 × 10 5 0.25 100 ---13.9 4.095 × 10 5 0.36 2.716 × 10 5 is reached (Federrath & Klessen 2012). Additionally, cells are refined if the change in pressure with respect to the neighboring cells is greater than 100%. To reduce the computational cost, we limit the refinement to regions representing a GMC using two passive scalars. Specifically, the outside of the clouds is minimally refined unless ≥ 1% of the cell mass stems from either of the clouds. These yield ∼ 10 8 − 2 × 10 8 computational cells, of which ∼ 1/3 are resolved at the maximum level.
To simulate star formation, the sink particle algorithm developed by Bleuler & Teyssier (2014) is employed. First, we identify gravitationally bound gas clumps using the Parallel HiErarchical Watershed (phew) algorithm (Bleuler et al. 2015) with a density of 0.1 ρ th , where ρ th is the threshold density for sink particle formation. Following Bleuler & Teyssier (2014), we adopt ρ th = πc 2 s /16G∆x 2 min , where c s is the speed of sound, G is the gravitational constant, and ∆x min is the size of the finest cell. The corresponding ρ th values for the runs with ∆x min = 0.25 pc and 0.125 pc resolutions are 3.0 × 10 −20 g cm −3 and 1.2 × 10 −19 g cm −3 , respectively. The accretion onto the sink particle is modeled using a flux accretion scheme. If the mass of a sink particle exceeds m res star = 100 M , a star particle of mass m res star is created, which generates ultraviolet and optical photons, as detailed below. The numerical methods and physical models of gas cooling and stellar feedback employed herein are identical to those employed in Kimm et al. (2019). The hydrodynamic equations are solved using the MUSCL scheme with a multi-dimensional MonCen limiter and HLLC Riemann solver (Toro et al. 1994). For radiative transfer, a GLF solver is adopted with a reduced speed of light approximation (10 −3 c, Rosdahl et al. 2013 [24.59, 54.42), and [54.42,∞). The explosion of Type II SNe is modeled using a mechanical scheme developed by Kimm et al. (2015). To bracket the impact of strong radiation feedback, we also examine an extreme scenario where additional pressure due to Lyα photons is included, assuming that the local velocity structures of the scattering medium are static (Kimm et al. 2018). The spectral energy distribution of each star particle is derived from the Binary Population And Spectral Synthesis model (Stanway et al. 2016, BPASS v2.0) assuming a Kroupa initial mass function (IMF) with lower and upper cutoff masses of 0.1 and 100 M , respectively. The initial metallicity of the clouds is set as Z ini = 0.1 Z , where Z = 0.02 is the solar metallicity. We adopt the metal yield of 0.075 for SN explosions. Gas cooling due to primordial species (HI, HII, HeI, HeII, HeIII, H 2 , and e − ) is included by solving the non-equilibrium chemistry (Rosdahl et al. 2013;Katz et al. 2017), and the cooling due to metallic species is included based on the Cloudy calculations (Ferland et al. 1998) and fine-structure line cooling by Rosen & Bregman (1995, 'cc07' model in ramses).
To model the head-on collision of two massive clouds, the velocity v col is added to the cloud in the lower left at t = 0.5 t ff (Figure 1, panel-(a)), and the other is kept stationary to minimize the numerical diffusion caused by the relative motion of clouds across the hydrodynamic grid (e.g. Robertson et al. 2010). In this study, v col is set to 30 km s −1 , motivated by the velocity range of collisions that induce locally enhanced star formation according to the recent idealized simulations of a disk galaxy embedded in a 10 11 M dark matter halo (Yoo log column density (g/cm 2 ) Figure 1. Formation of a massive compact cluster (i.e. potential GC) in the fiducial run (Run-Fid). Each panel displays the projected distributions of gas densities in the (128 pc) 3 (panels (a)-(c)) or (32 pc) 3 cube (panels (d)-(f)). The blue and red dots represent the stars formed during the first and second star formation events, respectively, whereas the green dots represent stars that do not belong to the massive cluster. Panel (a) shows the initial turbulent structures of the two GMCs with mass 3.6 × 10 5 M . Different panels correspond to the moment at which the GMCs coalesce (b), the moment at which half of the population A (i.e., premerger of the gas clumps) (d) or population B (i.e., postmerger) emerges (e), and the moment at which the GMCs are destroyed by the radiation feedback (f). The gray colors indicate the gas column density, as shown in the legend.
et al. 2020, priv.comm.) 1 . We also perform a simula-tion with no cloud collision (Run-Static, v col = 0) as a control sample. Each simulation is run for 11-14 Myr, at which point no additional stars are formed during the last ≥ 1 Myr. The star clusters are identified using the DBSCAN method (Ester et al. 1996). The final cluster mass M GC and radius R GC,h are measured 1 Myr after the end of star formation in the GC candidates, except in the Run-woRF-woSN run where star formation continues. These correspond to 10.7, 10.1, 12.3, 10.6, 11.0, and 12.6 Myr for Run-Fid, Run-Lya, Run-woRF, Run-HR, Run-Fid-woSN, and Run-Lya-woSN run, respectively. In the case of Run-Static where no GC candidate is formed, we measure the stellar mass at t = 10.7 Myr, which is 1 Myr after the end of star formation in the entire simulated volume. Table 1 summarizes the setup of our simulations and physical properties of the GC candidate.

Formation of a massive compact cluster
We begin by showing how a compact, massive star cluster forms when two GMCs, each of mass M GMC = 3.6 × 10 5 M , collide head-on with v col = 30 km s −1 (Run-Fid). As shown in Figure 1, the two clouds first coalesce at t = 5.3 Myr with a mild compression at their interface (panel-b). Because of the turbulent structures of the GMCs, star formation is not immediately triggered at the interface (c.f., Hunter et al. 2021), but the gas continues to collapse locally. At this stage, stars with the total stellar mass of 1.5 × 10 4 M are formed; they account for only ∼ 10 % of the total stellar mass produced by the end of the simulation. Note that none of these stars belong to the massive compact star cluster that emerges later. Most of the potential GC population starts appearing 1 Myr after the coalescence (t > 6.5 Myr). Several dense gas clumps that are located within a small region ( 10 pc) first form stars individually and then start to merge around t = 7.8 Myr (not shown, but see the blue dots in Figure 1). The star cluster formed from the multiple clumps is similar to the GC formation picture presented in the hydrodynamic simulations of dwarf galaxy mergers (Lahén et al. 2019). Neighboring clouds and diffuse gas continue to accrete onto the cluster and form stars (red dots) until the radiation feedback entirely disrupts the gas clump at t ≈ 10 Myr. This is corroborated in Figure 2  using blue and red lines (A and B, respectively), indicate that the compact cluster does not form in a single event but emerges due to the local collapse and merger, followed by gas accretion. The formation process occurs over a short timescale of ≈ 2-4 Myr, resulting in a final cluster of mass M GC = 1.08 × 10 5 M with a half-mass radius of 0.33 pc in the fiducial case. We note that the star formation efficiency per freefall time is considerably higher (∼ 20%) than in ob-served molecular clouds (1-2%) (Krumholz & Tan 2007;Utomo et al. 2018). It is generally argued that this observed low efficiency arises from turbulence-regulated star formation in clouds with relatively low column densities ∼ 100 M pc −2 and marginally high virial parameters α vir ∼ 1-2 (e.g., Krumholz et al. 2006Krumholz et al. , 2009. In contrast, our simulated GMCs have lower virial parameters (α vir = 0.7) and higher gas column densities (≈ 300 M pc −2 ), compared to GMCs located in the Local Group galaxies (Bolatto et al. 2008) or the Milky Way (Heyer et al. 2009). Indeed, in model Run-static where no collision is introduced, ≈ 8% of gas is converted into stars per free-fall time (Figure 2 (panel-a), grey dashed line), although no GC candidate is formed in Run-static. Such a high star formation efficiency is often obtained in the simulations of GMCs with high gas surface densities (e.g., Grudić et al. 2021;Kimm et al. 2022). We surmise that these more extreme physical conditions, typical of high-redshift galaxies, increase the likelihood of forming dense and compact clumps during cloud-cloud collisions, thus providing a favorable environment for viable GC candidate formation by boosting star formation efficiency.
In addition to the main cluster, the simulated GMCs form small compact star clusters with radius 1 pc. In the fiducial run, two stellar aggregations with > 40 particles are identified; their masses and half-mass radii are 9.6 × 10 3 M , 4.9 × 10 3 M and ≈ 0.2 pc, 0.1 pc, respectively. They form either in a single clump (former) or during the merging of gas clumps (latter); however, no subsequent gas accretion occurs in these small systems, unlike the formation process of massive GC candidates. Such stellar aggregations are also formed in the simulation without cloud collision (Run-Static). We note that these small clusters do not form in the higher resolution run (Run-HR), as gas fragmentation is better captured, yielding even smaller clumps in size and mass (e.g., Price & Federrath 2010;Federrath & Klessen 2012). Stellar feedback then efficiently disperses gas clumps before they merge and form a massive clump that would otherwise end up as a star cluster of mass 9.6 × 10 3 M . In contrast, the GC candidate is still present in Run-HR, denoting that the formation of the massive compact cluster is less affected by resolution.

Effects of radiation feedback
With the potential GC, we study the impact of radiation feedback on the stellar mass, gas inflow/outflow, and metal enrichment of the compact clusters. . Formation of a massive compact cluster in the run without radiation feedback (Run-woRF; top panels), with Lyα feedback (Run-Lya; middle panels), and with higher resolution (Run-HR; bottom panels). Each panel exhibits the projected distributions of gas densities in the (128 pc) 3 (the first column) and (32 pc) 3 cubes (the other three columns). From left to right, the panels show the coalescence of two GMCs, the first or second star formation episodes, and the massive gas expulsion by radiation (panels (b4) and (c4)) or supernova feedback (panel (a4)). The colored dots represent star particles formed in different phases, as in Figure 1.

Star formation history
ation feedback (Run-woRF) where the star particles do not produce photons. Thus, both photoionization and radiation pressure due to UV and the multiple scattering of IR photons are ignored, while still considering SN explosions. We find that an approximately twice more massive GC candidate forms in the Run-woRF run than in the fiducial run. In Figure 2, the stellar mass of the population A (stars formed before the one final massive clump appears) is slightly reduced, but this is because the mergers between gas clumps take place earlier than those in the fiducial run. The growth of stellar masses is not controlled in Run-woRF, and stars continue forming even after t 9 Myr at which point the radiation feedback blows the gas away in Run-Fid (see also Figure 3). Consequently, 31% of the two GMCs is converted into the GC candidate within ≈ 6 Myr in Run-woRF, while only 15% forms stars in the Run-Fid run. Additionally, the mass of the second or third massive stellar aggregations increases by a factor of 2-3 (3 × 10 4 or 1 × 10 4 M , respectively) in the absence of radiation feedback. This clearly demonstrates the importance of radiation in the formation of compact clusters. Thus, any simulation that does not account for radiation feedback will likely overestimate the stellar mass of a GC by a factor of ∼ 2.
The predicted stellar mass of a potential GC may decrease even further if additional nonthermal pressure operates on star-forming sites (see the second row in Figure 3). As demonstrated by Dijkstra & Loeb (2008, see also Smith et al. 2017Tomaselli & Ferrara 2021), the resonant scattering of Lyα photons can accelerate dustpoor gas to several tens of km s −1 , possibly suppressing the formation of metal-poor star clusters (e.g., Kimm et al. 2018;Abe & Yajima 2018). Since the column density of neutral hydrogen in our metal-poor star-forming clumps tends to be high (N HI 10 22 cm −2 ), Lyα photons can easily scatter 10 7 -10 8 times. Accordingly, the typical force multiplier of Lyα photons (M F ), which is a measure of the enhancement in momentum transfer compared to the single scattering limit, can be large (M F ∼ 10-50) (Kimm et al. 2018, see their Fig. 1). Thus, Lyα pressure may play an important role in GCforming sites. Strictly speaking, the number of resonant scatterings sensitively depends on the velocity structure of the ISM, but here the scattering medium is assumed to be static to maximize the impact of Lyα feedback. Therefore, the Run-Lya run should be considered as an extreme model where early feedback due to radiation is enhanced, disrupting star-forming clumps. More realistic simulations will likely yield results between the fiducial and Run-Lya model.
We find that the additional pressure prevents the gas from accreting onto the formation site and suppresses the total stellar mass growth. The stellar mass of population A is little affected because the gas clumps rapidly collapse before the Lyα pressure becomes effective. However, Lyα feedback from this population reduces the gas mass within the radius of 6 pc (10 pc) from the cluster center at t = 8 Myr from 8.9 × 10 4 M (11.4 × 10 4 M ) in the fiducial run to 6.4 × 10 4 M (8.8 × 10 4 M ), respectively. Moreover, as shown in Figure 2 (panel-(c)), Lyα pressure intervenes in the star formation of population B in ∆t ≈ 0.5 Myr (see the red line), which is significantly faster than that for the fiducial case, wherein it takes ∆t ≈ 1.0 Myr to quench star formation. The mass of population B is reduce by 33% in Run-Lya, and the total stellar mass of the potential GC decreases by 16% (M GC = 9.08 × 10 4 M ). Figure 4 presents the gas outflow and inflow rates measured at the spherical surface of radius 128 pc from the center of the two GMCs (gas plus stars). The surface is taken as large as possible in order to avoid any confusion due to the colliding motion of the two GMCs. By doing so, we confirm that the net inflow plus outflow rate is close to zero in Run-woRF-woSN. Outflows start appearing at t ∼ 7 Myr even before the massive compact cluster forms and destroys the gas clumps. This happens because the radiation from the underlying population in the neighboring region also drives the winds (see the green lines in Figure 2). After t ∼ 9 Myr, SN explosions from the underlying population also play a role in enhancing the gas outflow, as can be inferred from the comparison with the simulations without SNe (dashed lines in Figure 4). The dense gas in the potential GC-forming site is dispersed only at t ≈ 9.5 Myr, and SN explosions from the GC candidate are mainly responsible for the enhanced outflows at t 10.5 Myr, with velocities up to v ∼ 600 km s −1 at t = 13.0 Myr. Given that the outflow rates areṀ out ∼ 0.1 M yr −1 at t 10 Myr, the bulk mass loading factor may be estimated asṀ out / Ṁ star ∼ 3.5, where Ṁ star is simply taken as the total stellar mass divided by 5 Myr. In contrast, if we take the stellar mass of the underlying population, the mass loading factor would increase by a factor of four (Ṁ out / Ṁ star ∼ 13.8).

Gas inflow and outflow
The outflow rates in the fiducial run with radiation feedback are systematically larger than those of Run-woRF at t 12 Myr. The typical density and velocity of the outflows are n H ∼ 1-2 cm −3 and v out ∼ 30-100 km s −1 , respectively. Note that even in this case with a high star formation efficiency of ≈ 20%, the density of the gas outflows on the GMC scale is significantly smaller than the electron density (n e ∼ 300 cm −3 ) inferred from the metal emission lines in star-forming galaxies at high redshift (e.g., Sanders et al. 2016;Schreiber et al. 2019). We attribute this to the fact that photoionization heating overpressurizes the medium, efficiently lowering the density (e.g., Krumholz et al. 2007); in contrast shock compression due to SNe may play a role on the galactic scale (Schreiber et al. 2019). The outflowing density and velocity in the Run-Lya run further increase due to the destruction of the dense starforming clumps because of the strong radiation pressure. This increases the outflows by ∼ 20%. Furthermore, gas inflow rates are suppressed in the runs with radiation feedback, and even become almost zero within ∼ 3 Myr from the onset of star formation; this indicates that the impact of radiation is not limited to a small volume of star-forming clumps but extends to a large solid angle. The outflows become weak rapidly at t 12 Myr, as the majority of the GMC gas is blown out and the average density in the sphere with radius of 128 pc is reduced (n H ≈ 2 cm −3 at t = 13 Myr).
Note that the formation of a potential GC ends before their massive stars explode as SNe. This means that radiation can drastically alter the environment into which SNe propagate (e.g., Chevance et al. 2022). For example, SNe around t ≈ 11 Myr emerge at n H ∼ 10 4 cm −3 in the absence of radiation feedback (Run-woRF), but they explode at n H ∼ 0.1 cm −3 in the fiducial run (see yellow lines in Figure 6). Given that the final momentum from SNe (p SN ) is proportional to n −2/17 H (e.g., Thornton et al. 1998), the high-density environments denote that p SN from SNe that emerge early in Run-woRF (i.e., before they destroy the star-forming clouds) will likely be underestimated by ∼ 0.1/10 4 −2/17 ≈ 4, mitigating the effect of SNe on galactic scales.

Internal metal enrichment
As is well established, young GCs are composed of simple stellar populations with a small internal spread in [Fe/H] (e.g., Hollyhead et al. 2019). In general, the observed metallicity spread in Milky Way GCs is ∼ 0.045 dex (Carretta et al. 2009;Bailin 2019). We find that the presence of radiation feedback suppresses the internal metal enrichment due to SNe. In our simulations, Type II SN explosions occur between 3.5 and 40 Myr when stars more massive than 8 M evolve off the main sequence, but the formation of the star clusters . Metallicity distribution of the stars in the GC candidate. Different color-codes correspond to different runs, as indicated in the legend. Note that the simulation without radiation feedback over-predicts the internal metal enrichment due to SN explosions. We do not include the distribution from Run-HR, as it is very similar to that from Run-Fid. ends when their age is less than 3 Myr (t ≈ 9.0 Myr, Figure 2). Thus, no stellar particle is significantly selfenriched by the explosion of SNe inside the GC candidate. Although some of the stars formed before the collision explode at t 8.0 Myr, only a tiny amount of the ejecta is accreted onto the sink and all of the cluster stars exhibit metallicity close to the initial metallicity (Z ini ), i.e., ∆Z /Z ini < 0.01 ( Figure 5).
In contrast, the run without radiation feedback comprises a non-negligible fraction (16%) of the GC stars that are self-enriched (∆Z /Z ini > 0.01). Owing to the explosions of massive stars and efficient gas consumption, star formation is still quenched quite early (∼ 6 Myr) but the metallicity of the cluster stars increases up to 0.16 Z . Although this metallicity is not significantly large, care must be taken especially when modeling chemical mixing in multiple populations of GCs using hydrodynamic simulations without radiative transfer.
We note that the presence of radiation feedback does not necessarily indicate zero metallicity spread in GCs. First, any potential variation in the intrinsic metallicity of the colliding GMCs is neglected in this work (e.g., De Cia et al. 2021). If the initial metallicities of the GMCs are different, the cluster formed by their mergers will inevitably have an internal metallicity spread (e.g., McKenzie & Bekki 2021). Second, the star formation timescale could be longer than the SN timescale, if the GMCs are significantly more massive. For example, in a cosmological RHD simulation of a dark matter halo with mass ∼ 5 × 10 7 M (Kimm et al. 2016), the stellar metallicities of GC candidates increase up to 0.3 Z from zero metallicity. This is because the gas mass in the halo is an order of magnitude larger than the mass of the colliding GMCs and the star formation continues for ∼ 10 Myr, allowing the SN ejecta to be recycled for ∆t ∼ 6 Myr. Finally, as suggested in the literature (e.g., Kim & Lee 2018), multiple star formation episodes followed by metal enrichment on several hundreds of Myrto-Gyr timescales may well explain the complex chemical abundance patterns observed in GCs (Bastian & Lardo 2018;Carretta 2019). These diverse conditions will be the subject of future studies and our simulations should be regarded as controlled experiments for better understanding the detailed effects of radiation feedback in GC (candidate) formation.

Dominant radiation feedback process
Several types of radiation feedback occur in GMCs, such as photoionization heating or direct radiation pressure. To determine their relative importance in the suppression of star formation in GC-forming sites, we compute the maximum extent to which photoionization heating (r PH ), direct radiation pressure (r DP ), and Lyα pressure (r α ) balance the ISM pressure in Figure 6. These scales are computed as (Rosdahl & Teyssier 2015;Kimm et al. 2017) r PH = 3N ph,ion 4πα B 1/3 k B T ion P ther + P grav + P turb 2/3 (1) r DP = L abs 4πc (P ther + P grav + P turb ) (2) where N ph,ion is the number of ionizing photons radiated per unit time, T ion = 10 4 K is the temperature of the ionized medium, and α B is the recombination coefficient at T ion (Hui & Gnedin 1997). The thermal and turbulent pressures are calculated within a sphere of radius 2 pc from the cluster center, as P ther = (ρ gas /µm H ) k B T dV /V and P turb = ρ gas | v − v | 2 dV /V , respectively, where µ is the mean molecular weight. The pressure due to gravity is taken as P grav = − ρ g · rdV /V . To evaluate the effect of direct radiation pressure due to UV and optical photons, the luminosity absorbed by hydrogen, helium, and dust (L abs ) is computed. Additionally, the total Lyα luminosity emitted by gas within the same volume is computed from the recombinative and collisional radiation, as described in Kimm et al. (2018). The effect of nonthermal IR pressure is not presented here, as its contribution to the total pressure is negligible in our simulation. This is because the typical optical depth of the IR photons by dust is small (τ d 0.1) and they simply escape without being trapped in such metal-poor environments. Figure 6 shows that in the fiducial run, direct radiation pressure more effectively controls gas collapse than photoionization heating while forming the potential GC population (t 9 Myr). This is also evident in Run-woRF if r PH and r DP are measured, based on the luminosities expected from the stellar population (panel (b1)). In the run with Lyα, r α is even larger than r DP , indicating that Lyα pressure is the strongest form of radiation feedback early on. However, r α plunges at t ≥ 9.2 Myr once the cloud is fully ionized and Lyα freely escapes the system.
Specifically, in the early stage of cluster formation (6.5 t 7.5 Myr), the stellar mass grows and the volume-weighted gas density within the central 2 pc sphere decreases (panels (a2) and (c2) in Figure 6). Subsequently, the pressure due to gravity and turbulence drops, which increases r DP and r PH in the runs with radiation feedback. Note that the turbulent pressure and the pressure due to gravity are similar at t 9 Myr (before the dense clumps are disrupted), as the initial virial parameters of the GMCs are set to α vir = 0.7. Despite the decrease in density, the thermal pressure tends to build up or remain constant because the radiation feedback heats up the gas; however, P ther is still subdominant compared to P grav or P turb . In contrast, P ther in Run-woRF is smaller than that in the other runs by more than an order of magnitude, as the cloud temperature is kept cold (T ∼ 100 K). At t ≈ 8 Myr, the gas density increases again in Run-Fid and Run-Lya due to the merger with other gas clumps, causing temporary reductions in both r DP and r PH . Once the direct radiation pressure removes the gas from the central region, photoionization takes over and r PH rapidly expands at t ≈ 9 Myr 2 . Then, the gas density in the cluster-forming region reduces to a level similar to the ISM (n H 0.1 cm −3 ), Run-Lya r DP r PH r 6 7 8 9 10 11 12 (c2) Figure 6. The top panels display the maximum extent measured from the center of the potential GC to which photoionization heating (rPH, blue), direct radiation pressure (rDP, black), and Lyα pressure (rα, red) can overcome the ISM pressure in three different runs. In Run-woRF, where the radiation feedback is absent, rPH and rDP are computed assuming that stellar particles radiate and are presented by dashed lines. The bottom panels exhibit volume-weighted pressure due to gravity (brown), thermal motion (magenta), and turbulence (green). The orange lines indicate the volume-weighted (solid) or mass-weighted (dashed) densities. All quantities are measured inside a sphere of 2 pc radius from the cluster center. and the absorption of optical and UV light becomes insignificant (t ≈ 10 Myr). Figures 1 and 3 show that the formed clusters are generally very compact; the half-mass radius of the Run-Fid cluster is R c,h = 0.33 pc with a stellar mass of M GC = 1.083 × 10 5 M at the end of the simulation. The final mass of the cluster decreases by 16% in Run-Lya or increases by a factor of 2.07 in Run-woRF, but their R c,h are still comparable to that of the fiducial case. Furthermore, the cluster size is even more compact (R c,h = 0.23 pc) in the run with the higher resolution (Run-HR), demonstrating that the size may not yet have converged. This is likely because the turbulent structure of the cluster-forming region is not fully resolved. As noted in previous studies, the Jeans length needs to be resolved using more than 32 cells to capture the turbulent structure in the ISM (e.g., Federrath & Klessen 2012). However, herein the Jeans length in the cluster forming region is often resolved using only 8 or so cells (λ J /∆x min ∼ 1 pc/0.125 pc). Nevertheless, be-cause the cluster is formed by mergers of gas clumps and stars bearing angular momentum, we believe that the resultant stellar structure will not be arbitrarily compact. Moreover, since the massive star cluster is better resolved in Run-HR than in the fiducial case, the cluster size is unlikely to dramatically change. However, we note that the stellar mass of the cluster might be overestimated, given that an intermediate-mass black hole could form via stellar collisions in dense clusters (e.g., Katz et al. 2015). Figure 7 compares the compactness of the simulated clusters (filled circles) with that of the present-day Milky Way GCs (gray points). We employ the mass-to-light ratio of M/L V = 2.0 M /L V, , typical of GCs with ∼ 10 Gyr (Baumgardt 2017), to estimate the V-band luminosity of the simulated clusters. We remind the reader that such cloud-cloud collisions may occur at high redshift (z 5) during the formation of dwarf galaxies (e.g., Kim et al. 2018b;Sameie et al. 2022). Figure 7 shows that at these luminosities, the simulated clusters are more compact than the Milky Way GCs (0.2 R c,h 5 pc). How- ever, Gieles et al. (2010Gieles et al. ( , 2011 argued that the GCs are likely to lose up to ∼ 90% of their initial stellar mass due to stellar evolution, tidal disruption, and evaporation due to two-body relaxation, while their size increases by a factor of up to two (see also Larsen et al. 2012;Madrid et al. 2012;Kruijssen 2015;Rodriguez et al. 2022). If these effects are considered, as shown by dashed boxes in Figure 7, it is possible that the simulated clusters evolve into normal present-day GCs. Although a direct N-body simulation is needed to accurately follow the evolution of the GC candidate and to determine whether the structure can survive from tidal disruption, our results suggest that the cloud-cloud collision may be a viable route for the formation of GC.

Caveats
We present one of the first studies that attempt to understand the impact of radiation feedback on the formation of massive compact clusters using RHD (see also, Dobbs et al. 2022). However, several caveats should be considered when interpreting our results. First, although a high resolution of up to ∆x min = 0.125 pc is employed herein, along with an aggressive refinement scheme to capture turbulent interactions, the size of the simulated clusters does not fully converge. As discussed in the previous section, the adoption of resolution elements that can resolve the Jeans length using more than 32 cells may produce a slightly more compact cluster. Second, our star particle of mass 100 M represents a coeval cluster, but this should be modelled in a more realistic way by sampling the individual stellar mass for a given IMF (e.g., Geen et al. 2018). To mimic the uncertainty in the IMF sampling, additional simulations are performed with a larger stellar mass resolution of m res star = 1000 M , but the stellar mass of the potential GC is only changed by ∼ 20% and our conclusions on the formation process remain largely unchanged. Third, several physical processes that may be relevant for GC formation are neglected. Grudić et al. (2022a) showed that although they are not responsible for disrupting clouds, protostellar jets can play an important role in regulating massive star formation. Momentum-driven stellar winds may also disperse very dense gas in the vicinity of metal-rich O stars, although they are less effective than photoionization heating in destroying gas clouds with n H ∼ 10 4 -10 5 cm −3 (e.g., Dale et al. 2013). These imply that the stellar mass in our simulations may be overestimated. In contrast, Sakre et al. (2021) argued that strong magnetic fields delay gas collapse, increasing the number of massive dense cores when the field lines are aligned with the cloud collision axis. Finally, the determination of the stellar mass formed during GMC collisions may be affected by numerical diffusion. As pointed out by Pontzen et al. (2021), numerical diffusion present in Eulerian codes tend to delay gas collapse. This may increase the probability of the formation of a compact cluster due to the suppression of star formation in the moving GMC prior to the collision, although our higher resolution run, where diffusion is reduced somewhat demonstrates that this is not a dominant effect. The inclusion of magnetic fields, which are believed to play such a role (Sakre et al. 2021) might also alleviate this numerical issue.

CONCLUSIONS
Motivated by recent numerical and observational findings (e.g., Kim et al. 2018b;Tsuge et al. 2021), we simulated the formation of a massive compact cluster via the collision of two GMCs, each weighing 3.6 × 10 5 M , using RHD simulations. In particular, we focused on the impact of radiation feedback on the formation process of GC candidates and found that it can play an important role in determining the star formation histories, stellar masses, and metallicity distributions (see also Dobbs et al. 2022). Our results are summarized as follows.
1. A massive star cluster (≈ 10 5 M ) is formed when two GMCs collide. Roughly half of the cluster stars are formed in several gas clumps that collapse 1 Myr after the coalescence of the GMCs. The other half is produced in the massive clump that is formed by the merging of the early starforming clumps and the accretion of neighboring gas (Figure 1-3). In total, ≈ 15% of the gas clouds turns into a GC candidate in 2-4 Myr.
2. Radiation feedback efficiently controls the growth of stellar mass in cluster-forming gas clumps. Star formation is completed in ≈ 3 Myr in the presence of radiation feedback (Run-Fid), while the clusterforming site keeps accreting the neighboring gas clumps and forming stars for ≈ 6 Myr when the radiation feedback is neglected (Run-woRF). Thus, the cluster mass in Run-woRF is over-estimated by a factor of ≈ 2, compared to that in the fiducial case ( Figure 2).
3. Among the various radiation feedback processes considered throughout the formation of the potential GC, the dominant mechanism for the suppression of the gas collapse is the direct radiation pressure, i.e., momentum transfer via the absorption of LyC radiation. However, once the gas density drops to n H 10 3 cm −3 due to outflows driven by direct radiation pressure and gas consumption by star formation, photoionization is more dominant than direct radiation pressure in Run-Fid ( Figure 6). The impact of non-thermal radiation pressure due to IR photons is always negligible due to the low dust optical depth.
4. The additional pressure due to the multiple scatterings of Lyα photons partially inhibits the gas accretion onto the cluster-forming sites and suppresses the stellar mass by 16% in Run-Lya. The effect of Lyα pressure is only noticeable after a sufficient number of stars have formed and it progressively weakens as the star-forming region becomes ionized ( Figure 6).
5. The simulated clusters are much more compact than the present-day Milky Way GCs. The halfmass radius of the cluster is ≈ 0.3 pc at the end of the simulations, which is about three times smaller than the typical size of Milky Way GCs.
6. Radiation feedback suppresses the internal metal enrichment due to SN explosions, which is consistent with the observations of a simple stellar population in the star clusters younger than ∼ 2 Gyr (Bastian & Lardo 2018). However, 15% of the cluster stars are self-enriched by SNe in the run without radiation feedback.
Our numerical experiments focused on one way of forming a compact, massive star cluster during a headon collision between two GMCs. However, any process that triggers rapid gas collapse and star formation warrants investigation, including galaxy mergers (e.g., Whitmore & Schweizer 1995;Renaud et al. 2015;Kim et al. 2018b;Lee et al. 2021), to fully understand the GC formation process. Although the star cluster mass may be overestimated herein due to the absence of other forms of feedback, such as protostellar jets or stellar winds, the impact of direct radiation pressure and photoionization is likely to be persistent in other GC formation scenarios. As shown in Section 3.2.3, the early feedback affects the local environment into which SNe propagate, and hence the metal enrichment in the starforming clouds and neighboring medium. This implies that chemical mixing is complex in star-forming sites and warrants further studies on the detailed chemical evolution of GCs (e.g., Tenorio- Tagle