Star cluster formation and feedback in different environments of a Milky Way-like galaxy

It remains unclear how galactic environment affects star formation and stellar cluster properties. This is difficult to address in Milky Way-mass galaxy simulations because of limited resolution and less accurate feedback compared to cloud-scale models. We carry out zoom-in simulations to re-simulate 100-300 pc regions of a Milky Way-like galaxy using smoothed particle hydrodynamics, including finer resolution (0.4 Msun per particle), cluster-sink particles, ray-traced photoionization from O stars, H$_2$/CO chemistry, and ISM heating/cooling. We select $10^6$ Msun cloud complexes from a galactic bar, inner spiral arm, outer arm, and inter-arm region (in order of galactocentric radius), retaining the original galactic potentials. The surface densities of star formation rate and neutral gas follow $\Sigma_{SFR} \propto \Sigma_{gas}^{1.3}$, with the bar lying higher up the relation than the other regions. However, the inter-arm region forms stars 2-3x less efficiently than the arm models at the same $\Sigma_{gas}$. The bar produces the most massive cluster, the inner arm the second, and the inter-arm the third. Almost all clusters in the bar and inner arm are small (radii<5 pc), while 30-50 per cent of clusters in the outer arm and inter-arm have larger radii more like associations. Bar and inner arm clusters rotate at least twice as fast, on average, than clusters in the outer arm and inter-arm regions. The degree of spatial clustering also decreases from bar to inter-arm. Our results indicate that young massive clusters, potentially progenitors of globular clusters, may preferentially form near the bar/inner arm compared to outer arm/inter-arm regions.


INTRODUCTION
Star formation takes place in giant molecular clouds (GMCs) with most stars forming in clusters or associations (Lada & Lada 2003). How these clusters/associations form is still an open problem, as is the cause of differences in their properties. Young massive clusters (YMCs; masses >10 4 M ⊙ , radii ∼1 pc) are of particular interest as they may be the progenitors of globular clusters (Portegies Zwart et al. 2010;Longmore et al. 2014). In addition to the cluster properties, the gas itself is influenced by the newly formed stars. Once massive stars (>8 M ⊙ ) form, they become a source of feedback by releasing energy and momentum into the interstellar medium (ISM), changing the gas dynamics while stars are still forming. This can affect star formation rates by dispersing gas reservoirs (Walch et al. 2012) or by compressing them to form new stars (Elmegreen & Lada 1977;Whitworth et al. 1994).
Star formation, feedback, and cluster properties may depend on their birth environment. The observed relation between gas surface density and star formation rate (SFR) surface density depends on galactocentric radius, with SFRs being higher at smaller radii (Bigiel et al. 2008). Numerical simulations have extensively shown that feedback depends on initial cloud conditions such as mass (Dale et al. 2014;Howard et al. 2017 (Fukushima et al. 2020;Ali 2021), structure (Walch et al. 2013;Zamora-Avilés et al. 2019), gravitational boundedness (Howard et al. 2016), and turbulence (Geen et al. 2018;Guszejnov et al. 2022). The most important feedback mechanism on cloud scales appears to be photoionization, which heats gas from ∼10 2 K to 10 4 K, creating a pressure gradient between ionized gas and the neutral ISM. While there are still many uncertainties, photoionization may dominate over other presupernova (SN) mechanisms such as stellar winds (Geen et al. 2021;Ali et al. 2022) and radiation pressure (Kim et al. 2018;Ali 2021). These mechanisms set the structure into which SNe explode, potentially creating low-density channels through which energy can escape (Lucas et al. 2020;Bending et al. 2022).
Results from cloud-scale studies need to be placed in the larger galactic context. GMC evolution is influenced by galaxy-scale potentials, shear, and cloud-cloud tidal forces Jeffreson et al. 2020). Observations in NGC 300 show that feedbackrelated pressure terms exhibit a slight dependence on galactocentric radius, indicating that feedback becomes more powerful at larger radii (McLeod et al. 2021). H ii regions at small galactic radii may be confined by higher ambient pressures (Barnes et al. 2020;Della Bruna et al. 2022), meaning H ii regions in the disc may be able to expand to greater sizes compared to those nearer the centre.
However, it is less clear how specific galactic structures such as bars and spiral arms affect the star formation and feedback processes.
Individual clouds in the Central Molecular Zone of the Milky Way (the innermost 500 pc) are observed to be significantly less starforming than expected given their high densities (Longmore et al. 2013;Kauffmann et al. 2017). This may be due to galaxy-scale potentials causing strong shear (Kruĳssen et al. 2019). Yet galactic centres do contain young massive clusters, including, in our own Galaxy, the Arches and Quintuplet clusters (Longmore et al. 2014). The Galactic bar and its intersections with spiral arms might also host YMCs (Davies et al. 2012;Ramírez Alegría et al. 2014), including W43 (Nguyễn Luong et al. 2011;Carlhoff et al. 2013). While disc clouds may have lower masses or mean densities than clouds in the centre, such regions still manage to form YMCs such as NGC 3603, potentially through cloud-cloud collisions (Fukui et al. 2014;Liow & Dobbs 2020).
Numerical simulations which explore these processes can broadly be split into two types: galaxy-scale simulations which model the evolution of a whole Milky Way-mass galaxy over 100's of Myr and follow the interaction between GMCs and large-scale structures such as spiral arms (see the review by Naab & Ostriker 2017;also, Agertz et al. 2013;Dobbs & Pringle 2013;Smith et al. 2020;Jeffreson et al. 2020;Pettitt et al. 2020a;Keller et al. 2022). The second type follows the cloud-scale evolution over 3-10 Myr and follows star formation and feedback without external influences. In these models, it is computationally feasible to include high-resolution pre-SN feedback methods, but the initial conditions are idealised, usually in the form of turbulent spherical clouds (see the review by Dale 2015;also, Walch et al. 2013;Dale et al. 2014;Grudić et al. 2018;Kim et al. 2018;Ali et al. 2018;Geen et al. 2018).
In this paper, and previous papers (including Bending et al. 2020;Ali et al. 2022;Dobbs et al. 2022a;Herrington et al. 2023), we attempt to bridge this gap by extracting GMC complexes from galaxy simulations (thus starting with more realistic density and velocity distributions; Rey-Raposo et al. 2017), and using feedback methods such as ray-tracing which are usually limited to smaller scales. These zoom-in simulations retain the galactic potentials and include multiple clouds, thus including tidal forces and shear which isolated cloud models would neglect. By extracting GMC complexes (of mass 10 6 M ⊙ and size 100-300 pc) from different parts of a galaxy, we explore the impact of galactic environment on star/cluster formation and feedback. In section 2, we describe the zoom-in method and the implementation of star formation and feedback. We present our results in section 3, discussion in section 4, and conclusions in section 5.

NUMERICAL METHODS
We use the smoothed particle hydrodynamics code sphNG, which originated with Benz et al. (1990) and Benz (1990), with substantial modifications by Bate et al. (1995) and Price & Monaghan (2007) such as the inclusion of sink particles and magnetic fields (although we do not include the latter here). Our models include self-gravity, ISM heating/cooling, H 2 and CO chemistry (Glover & Mac Low 2007;Dobbs et al. 2008), and the galactic potentials described in section 2.1. We use the methods introduced by Bending et al. (2020) for cluster-sink particles and photoionizing radiation via ray-tracing. We also include supernovae as described by Bending et al. (2022). We summarise the sink method in section 2.2 and the feedback processes in section 2.3.

Initial conditions
We set up our initial conditions by extracting a region from a galaxy evolution model, enhancing the resolution, then re-running the region with the higher accuracy methods for sinks and feedback (i.e. a zoom-in simulation). The galaxy is based on the Milky Way and includes analytic potentials for a bar and four spiral arms. We use a modified version of the BrSp4 model from Pettitt et al. (2020a) taken at 340 Myr post-initialisation. This simulation was carried out using the SPH code gasoline2 (Wadsley et al. 2004. It included stellar feedback in the form of winds from evolved low-mass stars (<8 M ⊙ ) and supernovae (type II and type Ia; Stinson et al. 2006;Pettitt et al. 2017). Main sequence feedback from high-mass stars was not included. We retain the bar and arm potentials in our zoom-in simulations using sphNG. The bar (Wada & Koda 2001) has a scale length of √ 2 kpc and a pattern speed of 60 km s −1 kpc −1 . The spiral arms (Cox & Gómez 2002) have a pitch angle of 15 • and pattern speed 20 km s −1 kpc −1 . There is also an axisymmetric potential for the combined disc + bulge + halo (Pettitt et al. 2014). This model differs to that used in Pettitt et al. (2020a) in that it includes a greater gas-mass resolution of 600 M ⊙ per SPH particle and a gravitational softening scale of 5 pc (compared to 1500 M ⊙ and 50 pc presented in the aforementioned paper). The global dynamics and structure are essentially identical to that of the model studied in Pettitt et al. (2020a) and are not discussed here.
We select gas particles in a region of the galaxy and increase the resolution using the particle-splitting method of Bending et al. (2020). In the following, we use sphNG and retain the original galactic potentials, but we do not take any of the original star particles. For each model, we do the resolution enhancement in two stages. This reduces the effect of grid artefacts in the final particle setup. The first enhancement takes the region (size of the order of ∼1 kpc) and increases the resolution from the base resolution of the galaxy simulation to ∼13 M ⊙ per particle. The galactic position of each region at this stage of the process is shown in the left-most panel of Fig. 1 overlaid on the original galaxy. We evolve this region (without feedback) for approximately 0.5 Myr to allow the particles to settle. Next, we select a sub-region (size 100-300 pc) and enhance to ∼0.43 M ⊙ per particle. This is used as the initial condition to evolve with cluster-sinks and stellar feedback.
The parameters of each initial condition are listed in Table 1, sorted by galactocentric distance. One region explored is in the bar, two are in spiral arms (one inner arm and one outer arm), and finally one is in an inter-arm region. The initial mean density is largest for the bar and decreases with galactocentric distance, with the bar being 10 times denser than the outermost region (inter-arm model). The panels to the right of the galaxy in Fig. 1 show the models 2-4 Myr after the onset of ionization.

Cluster-sink particles
The zoom-in simulations form cluster-sink particles which represent (sub-)clusters of stars. Sink formation is based on the criteria of Bate et al. (1995). Gas particles above a density threshold of 1.2 × 10 4 cm −3 are tested to see if the neighbourhood of ∼ 50 particles is collapsing and converging. If so, the particle and its neighbours are converted to a sink particle. Sink formation is forced for densities above 1.2 × 10 6 cm −3 . The sink accretion radius is 0.1 pc and the sink merger radius is 0.03 pc.
The method for converting sink masses to a stellar population is based on a method by Geen et al. (2018)   simulation, we create a list of massive stars. This is done by sampling from a Kroupa (2001) initial mass function (IMF) up to a total mass of 3 × 10 6 M ⊙ and grouping stars into mass bins, keeping a list of the bin indices for masses above 18 M ⊙ . The same list and ordering is used for all models presented here, as well as in Bending et al. (2020) and Ali et al. (2022). We keep track of the mass accreted by each sink. When the total mass accreted over all sinks reaches 305 M ⊙ , we take the next massive star from the list and assign its properties (e.g. ionizing flux) to a chosen sink -this is the sink with the most mass not made up of massive stars. This is only done if the sink is massive enough to accept the star; otherwise, we wait until the next time step. The fraction of the sink mass that is available for star formation is 50 per cent -the rest is assumed to be a gas reservoir. Stars are in bins according to spectral type with representative masses and ionizing fluxes listed in table 2 of Bending et al. (2020).

Stellar feedback
We use a ray-tracing method for calculating photoionization equilibrium along lines of sight (LOS) between gas particles and sinks. This is a similar method to Dale et al. (2007). A full description is available in Bending et al. (2020). For each gas particle, we calculate the ionizing flux it receives from all ionizing sources, taking into account the reduction along the LOS. All particles with a smoothing length which overlaps with the LOS are included, with quantities interpolated at the position on the LOS (rather than the particle position, as Dale et al. do). The photoionization rate is balanced with the recombination rate at the gas particle density. We use the on-the-spot approximation with case B recombination coefficient B = 2.7 × 10 −13 cm 3 s −1 and ionization temperature of 10 4 K. If a gas particle receives ionizing radiation from multiple sources, the column density contributions from all lines of sight are divided by the number of sources contributing the flux. This is described in Herrington et al. (2023).
To reduce the computational time, we set the maximum LOS distance to 100 pc -this is tested by Bending et al. (2020) in a 500 pc region. In that case, the longer range ionization only had a small effect on top of the limited LOS model, sweeping up shells near the boundaries where massive stars did not form. Here we investigate smaller regions (100-300 pc, closer to the LOS radius), and the sink particles are more evenly spread out, meaning we do not expect this limit to change our results significantly. We also group ionizing sinks into nodes if they are close together, which reduces the number of ionizing 'sources' that a gas particle needs to loop over. For the most ionizing sink, we find the minimum radius at which the average ionization fraction of enclosed gas drops below 90 per cent. Any other sink that sits within half this radius is grouped with that sink. We repeat with the next most ionizing sink that has not already been grouped, until all sinks have been grouped or tested. The summed flux of each sink group propagates from its centre of flux. We find that this reduces the total number of ionizing sources in the latter stages of our simulations by a factor of between 2 and 3. Testing with and without this optimisation shows only a small effect on the total amount of ionized gas (<1 per cent) and a negligible effect on H ii region morphology.
We also include supernovae (SNe) using the method by Dobbs et al. (2011), which was updated for cluster-sink particles by Bending et al. (2022). When a star above 18 M ⊙ becomes old enough to explode as a SN, we insert energy around the host sink inside a radius which encompasses its 80 nearest particles. This radius is used to calculate the age, temperature, and velocity of a SN bubble in the snowplough phase, assuming each SN contributes 10 51 erg. This solution provides the energy to be inserted inside the radius as a combination of thermal and kinetic energy. We do not include SNe from less massive stars as their life times are longer than the simulation run times.
We do not include stellar winds in this paper, as their impact on the gas is negligible compared with photoionization (Ali et al. 2022).
Their impact on cluster properties may be marginally more important compared with gas properties, so this is planned for future models. Fig. 2 shows, for each region, the time evolution of column density. Four snapshots are presented between 0.5 Myr after the onset of feedback and the end of the simulation run time. The end time is typically limited by the small time steps caused by sink dynamics (in particular the bar model) or heating from SNe. The effect of rotation around the galactic centre can be seen most clearly for the bar region, where the gas structure starts almost vertical and rotates almost 90 degrees over the next 3.2 Myr. The inner arm also shows some rotation, but this is marginal compared to the bar. The outer arm and inter-arm regions do not show rotation over the time scales modelled here. The bar has a burst of star formation in the first Myr and is dominated by a large, dense cluster which forms from the densest gas. Star formation in the inner arm and outer arm models occurs along the length of the central arm structure, with clusters forming in a chain-like pattern. Star formation in the inter-arm region is more sparsely distributedthis has three main sites of star formation, with two close together on the left and one ∼ 100 pc away to the right. The inter-arm regions bear the closest resemblance to isolated cloud models. Fig. 3 shows the time-evolution of star formation efficiency (SFE), star formation rate surface density (Σ SFR ), and the cumulative star mass. The SFE is defined in terms of the total sink and gas mass,

Star formation
where the stellar mass is half the sink mass according to the clustersink prescription described in section 2.2. Σ SFR is the star formation rate divided by the rectangular area which contains 99 per cent of the neutral gas mass, with the origin at the centre of mass, as viewed in the -(top-down) plane, i.e.
This is calculated over time intervals Δ = 0.047 Myr and the dimensions are also re-calculated at every SFR measurement time. Δ sinks is the change in sink mass over the time interval Δ . Zero time in the plot is defined to be when ionizing radiation is first emitted. The star symbols show when supernovae explode, the first of which occurs after 4.5 Myr.
The final SFE decreases with galactocentric distance, with the bar model having the highest SFE of 25 per cent by the end of its runtime of 3.2 Myr. For the models which undergo at least one supernova event, the SFEs just before the first SN (in per cent) are 16 (inner arm), 12 (outer arm), and 6 (inter-arm). This shows that inner regions form proportionately more stars in the same amount of time. For most of the evolution, Σ SFR is also ordered by galactocentric distance as rates get lower with larger distance. However, this is slightly different in the first Myr, where the two arm regions fluctuate over each other. Similarly, from 4 Myr, the two outermost regions (the outer arm and inter-arm) overlap before the latter overtakes the former. By 5 Myr, there is not much variation between the three non-bar regions, with the range in Σ SFR being within a factor of 2. The peak values of Σ SFR are reached at 0.8 Myr (bar), 1.8 Myr (inner arm), 1.0 Myr (outer arm), and 4.2 Myr (inter-arm); the latter model however shows two peaks of star formation, with the first peak occurring around 1.5 Myr. The other three models experience a burst of star formation early on, before declining -star formation is spread more evenly in position across these regions, while the inter-arm has two distinct sites of star formation separated by ∼ 100 pc. All four models reach total stellar masses above 10 5 M ⊙ , with the bar being the first model to achieve this at 0.85 Myr and the inter-arm region being the last at 4.1 Myr.
The top panel of Fig. 4 shows the relation between Σ SFR and Σ gas , the gas mass surface density of neutral gas (equivalent to H i + H 2 ). For each model, we plot a data point at every simulation dump time, starting from 0.5 Myr after the first ionizing source starts radiating; we exclude earlier points to avoid the initial large scatter in SFR before this time due to the small number of sinks (c.f. the middle panel of Fig. 3). The black dashed line shows the power-law fit to the data, which is Σ SFR ∝ Σ 1.3 gas . This is found by applying a least-squares method to the data in log-space. The index agrees with the standard Kennicutt-Schmidt index of 1.4 ± 0.15 (Kennicutt et al. 2007), especially at later times (lower Σ gas ) for the bar, inner arm and outer arm models when the points evolve to follow the power law line. The Ali et al. (2022) spiral arm zoom-in is also included here, and is shown in green. However, the inter-arm model is shifted down compared to the other arm regions modelled in this paper -at early times, Σ SFR is lower by a factor of 2-3. This model shows a double bump in SFR without tailing downward at later time (lower Σ gas ) like the arm and bar models; this is shown more clearly in Fig. 3 as a function of time. The bar region is located higher up the power-law line than the arm regions, which generally all have lower densities and lie at similar points in the figure.
For comparison, observed regions from Bigiel et al. (2008) are shown as orange diamonds -these are from 7 nearby spiral galaxies with 750 pc resolution. Observations by Pessa et al. (2021) from 18 galaxies at a resolution of 100 pc are also plotted (assuming Σ mol,gas ≈ Σ neu,gas ). The models lie above the observational dataour star formation rate surface densities are a factor of ∼100 higher. Lines of constant gas depletion time scale (∝ Σ gas /Σ SFR ) are plotted as dotted lines. The models have gas depletion time scales between 10-30 Myr, while the Bigiel et al. (2008) and Pessa et al. (2021) data show time scales of 1-10 Gyr, with some of the latter regions approaching 100-300 Myr. Similarly, 1.5 kpc resolution observations by Sun et al. (2023) have depletion time scales above 1 Gyr. Milky Way regions with higher resolution from Heiderman et al. (2010) are shown in magenta. Their values of Σ SFR are closer to the simulations than the Bigiel et al. (2008) results, but have higher Σ gas as they are individual star-forming clouds (sizes smaller than Orion) rather than cloud complexes. We test how the pixel size and sink parameters affect the results in appendix A.
The bottom panel of Fig. 4 shows the star formation relation in terms of volume density (dividing by ) instead of surface density -this now includes the height of the region in the -dimension as well (above/below the galactic plane). This is SFR against gas . In observations, the main difference between the volume density law and the surface density law is that the latter is well known to have a break at low Σ gas , while the volumetric law shows indications of being the same for all gas (for the same method). The simulations lie along the dashed line which shows the power law fit to the models, SFR ∝ 1.3 gas . The index of 1.3 agrees with one of the values calculated by Bacchini et al. (2019a,b), who use surface density measurements in 12 nearby galaxies and turn this to a volumetric quantity by calculating the radius-dependent scale height of a disc in vertical hydrostatic equilibrium. They measure two power law indices, 1.3 (when using a constant SFR scale height of 100 pc) and 1.9 (when using a radially varying SFR scale height).

Cluster identification
We study the properties of the clusters formed in our simulations using two different methods, INDICATE (Buckner et al. 2019) and HDBSCAN (Campello et al. 2013). Both methods suggest the same trends in the properties of clusters with galactic region, namely that smaller denser clusters occur preferentially in the inner arm and bar regions, and larger clusters or associations in the outer arm and inter-arm regions. For this section, we analyse the clusters in each simulation at the time when the first supernova occurs (4.5 Myr after the onset of feedback), except for the bar region which is analysed at 3 Myr.

INDICATE (INdex to Define Inherent Clustering And TEndencies;
Buckner et al. 2019) is a local indicator of spatial association. This means that rather than finding discrete clusters like (H)DBSCAN (see section 3.2.2), it assigns an index to each point in a dataset that describes the spatial distribution in its local neighbourhood. The index has a range of 0 ⩽ 5 ⩽ tot −1 5 , where tot is the total number of points in the dataset and higher values represent greater degrees of association. INDICATE calibrates 5 against random distributions to identify the minimum value, sig , which denotes a point is spatially clustered (rather than randomly distributed). This value is defined as three standard deviations greater than¯r andom 5 i.e. sig =¯r andom 5 + 3 . Statistical testing by the authors has shown INDICATE to be independent of dataset size, shape, and density; robust against edge effects and outliers; and valid for sample sizes ⩾ 50 that are up to    gions in Fig. 5. The algorithm is applied to the 3D sink positions, with the figure showing the results projected onto the -plane. The value of 5 (colour scale) denotes the degree of association for sink particles that have been identified as spatially clustered, while randomly distributed sinks are plotted in grey. Sinks with higher values of 5 are more clustered. Fig. 5 shows that the highest values occur in the bar, followed by the inner arm, the outer arm and inter-arm. For sinks identified as being clustered, the median value of 5 in each model is 87 (bar), 92 (inner arm), 59 (outer arm), and 60 (inter-arm). Thus we see that denser, tighter clusters, are found in the bar and inner arm compared to the inter-arm and outer arm.

(H)DBSCAN
DBSCAN (Ester et al. 1996) is an algorithm commonly adopted observationally to identify clusters (  . The colour bar denotes the index value, 5 , of sinks found to be spatially clustered. Sinks with higher 5 are more clustered. Points in grey denote sinks found to have a random distribution. Sinks in the bar region are the most clustered, followed by the inner arm, outer arm, and inter-arm regions. Figure 6. Clusters of sink particles identified for each model using HDBSCAN (top-down projection). The different colours denote discrete clusters found, with members of the same cluster sharing colours. Particles which are not associated with any cluster are plotted in black. The clusters tend to be more spatially extended in the outer arm and inter-arm models, compared to the bar and inner arm models. core/border members of a cluster, or noise. is the radius around each star that is searched for neighbouring stars, and MinPts-1 is the minimum number of neighbours needed for to be a core member of a cluster. If has less than this minimum number of neighbours but is within the search radius of a core star it is considered a border star, else noise. We chose MinPts=30 simply as measure of the size of a cluster which is well resolved in the simulation. We originally tried using DBSCAN to identify clusters, ideally keeping MinPts and the same so that clusters could be compared consistently across different models, but we could not find uniform values of the parameter across all models. Using the method proposed by Ester et al. (1996) to determine , with MinPts=30, gave values of 0.39, 0.55, 0.67 and 0.88 pc for the bar, inner arm, inter-arm region and outer arm respectively, indicating that the bar and inner arm contain smaller clusters compared to the other regions. With these values of , the bar and inner arm preferentially produce smaller, denser clusters whereas the outer arm and inter-arm regions produce spatially larger, lower density clusters.
Although DBSCAN is commonly used in the observational literature, a disadvantage is that it is designed to find clusters of similar densities. This is especially true for datasets such as ours, in which clusters of different densities exist side-by-side in the same region, and different regions have different densities as well. HDBSCAN (Campello et al. 2013) is a successor to DBSCAN which utilises a hierarchical clustering approach to find clusters in different density regions. This increased sensitivity makes HDBSCAN a more effective algorithm for recovering clusters in observational datasets (Hunt & Reffert 2021). HDBSCAN does not require the parameter to be user-defined as it is essentially an implementation of DBSCAN which varies this value. Instead the main input to HDBSCAN is min_cluster_size, indicative of the minimum number of points in a cluster. A second input is min_samples, which is a measure of how strict the cluster assignment is. We found HDBSCAN produced satisfactory clusters without the need to alter the parameter, and unlike DBSCAN, we use the same parameters for HDBSCAN to find clusters across all the models.
To identify the best choice of input parameters for HDBSCAN, we determined the median INDICATE values of the clusters identified across the different models, giving us a range of values for min_cluster_size and min_samples which produced similar 5 values (note that there is some degeneracy between the two values, so for example increasing min_cluster_size and min_samples both tend to produce larger clusters). We also checked the clusters identified with the different parameters by eye, and compared the peaks identified with INDICATE (see Fig. 5) with the resultant clusters. We found that the inner arm, inter-arm, and outer arm favoured smaller values of min_cluster_size and min_samples, whilst the bar favours larger values. This again indicates that the bar model produces denser clusters with more particles (similar to the INDI-CATE results). We selected the largest values of min_cluster_size and min_samples within our optimal range which did not spuriously group particles together which were not clusters by eye in the inter-arm, inner arm and outer arm models. These values then overlapped with the lower range of optimal values for the bar model. Overall this approach gave values of min_cluster_size= 55 and min_samples= 40.
We show the clusters picked out with the HDBSCAN algorithm using these parameters for the different models in Fig. 6. As with INDICATE, HDBSCAN uses the 3D sink positions, with the figure showing a 2D projection. Although the scales vary slightly between the different panels, the figure indicates that more spatially extended clusters are found in the inter-arm model, and to some extent the   (2021) as a 2D histogram. The solid line shows the fit to the observations. Bottom -fraction of clusters in different radius bins. Clusters with radii > 5 pc are predominantly found in the outer arm and inter-arm regions, whilst no spatially larger clusters are found in the inner arm region. outer arm model, whereas more compact clusters are found in the inner arm and bar models.

Cluster masses and radii
Having identified the clusters, we show their effective radii and masses in the upper panel of Fig. 7, where the data from our models is plotted over observational data from Brown & Gnedin (2021). The effective radius we use is the half-mass radius, which is comparable to the observed radius. Similarly to Dobbs et al. (2022b), the points collectively indicate very similar trends to the observed data, and a similar increase in radius with mass. As also found in Dobbs et al. (2022b) the points have a slightly larger spread compared to the observational data. As mentioned in section 3.1, the star formation rates are higher than would be expected, so the clusters tend to be at the high mass end of the observational data. We find that the clusters with the largest radii form in the inter-arm and outer arm models -this is also indicated by eye in Fig. 6. With radii of around 10 pc, these objects are more comparable to local observed associations (Portegies Zwart et al. 2010) than clusters (note that smaller associations may be part of much larger ∼100 pc regions or associations (Wright 2020) but here we compare with the ∼ 5 pc size associations discussed in Portegies Zwart et al. (2010)). The clusters in the bar and inner arm models appear to follow fairly well the observed distribution. The clusters in the outer arm however exhibit relatively large radii, and taken in isolation would exhibit a much steeper radius mass relation compared to the observational dataset.
In the lower panel of Fig. 7, we plot the frequency of clusters with different radii for the different models. Again the bar and inner arm models contain spatially smaller clusters compared to the outer arm and inter-arm regions. The inner arm in particular contains no clusters of radii > 5 pc, whereas nearly half the clusters in the outer arm, and over a third in the inter-arm region, have radii > 5 pc.
We note that the tendency of the spatially largest clusters, for a given mass, to occur in the outer arm and inter-arm models is independent of our choice of algorithm or input parameters for HDB-SCAN. Using DBSCAN with the above values (see section 3.2.2), or choosing lower values for min_cluster_size and min_samples with HDBSCAN, tends to break up the clusters more, and the points are shifted to lower masses; in which case the larger clusters for the outer arm and inter-arm models are shifted to the top left area of Fig. 7 (upper panel).
We also see from Fig. 7 (upper panel) that the most massive cluster is formed in the bar, then the second most massive in the inner arm, followed by the inter-arm region, then the outer arm. Again this trend is fairly robust to the choice of algorithm and input parameters. We can compare the cluster masses and radii with clusters in the Milky Way. Although there is not a complete map of clusters in our Galaxy, Portegies Zwart et al. (2010) list the main clusters (or YMCs) and associations in our quadrant. At the end of the bar lies RSGC02 which is the most massive cluster (using phot from table 2 of Portegies Zwart et al. (2010)), at 4 × 10 4 M ⊙ (and RSGC01 and RSGC03 are close by with similar masses). The next most massive is Westerlund 1, in an inner spiral arm. Between arms, on a minor spiral arm, lies Orion which according to the table has a combined mass of 2 × 10 4 M ⊙ , whilst NGC 3603 just outside the solar circle has a mass of 1.2 × 10 4 M ⊙ . Therefore the trend in cluster mass with region seen in the simulations is the same as seen in our Galaxy.

Cluster rotation and expansion
We measure the bulk rotation of the clusters identified in section 3.2. We use a similar method to Ballone et al. (2020) and Verliat et al. (2022). For each cluster, we identify the centre of (sink) mass and calculate the angular momentum of each sink L = r × v, where r and v are its position and velocity, respectively, relative to the centre of mass. We then calculate the mean angular momentum ⟨L⟩ and rotate the reference frame so that the new ′ axis is parallel to ⟨L⟩, meaning this is defined to be the rotation axis of the cluster. For each sink, we calculate the azimuthal velocity component = v ′ ·ˆ′, whereˆ′ is the azimuthal unit vector around the ′ axis. The angular velocity is then = / , where is the distance from the ′ -axis. The radial velocity component is = v ′ ·r ′ . In Fig. 8 Note that while clusters in the inter-arm may have higher on average compared to the inner arm, they are also larger in size (see Fig. 7), hence the division by results in a smaller . The median radial velocities show that the majority of clusters in each region are expanding, except in the inner arm where most (56 per cent) are contracting. Radial velocities rarely exceed ±2 km s −1 , which is consistent with the observed clusters listed in table 2 of Kuhn et al. (2019). The only cluster with statistically significant rotation in the Kuhn et al. (2019) dataset is Tr 15 in the Carina Nebula, which lies in a spiral arm about 2.4 kpc from the Sun (Shull et al. 2021). It has a median = (1.7 ± 0.5) km s −1 , which is consistent with the velocities we find in our models, but less so with the outer arm region where most of the clusters rotate more slowly than this. R136 in the Large Magellanic Cloud has a mean = (3 ± 1) km s −1 and = (0.75 ± 0.22) Myr −1 (Hénault-Brunet et al. 2012), which fits with the bar/inner arm values, but is on the extreme end of the outer arm/inter-arm results.

Mass flows
We investigate the dispersal or replenishment of gas under the influence of feedback by measuring the radial mass flux near clusters. First, we identify the cluster with the highest total ionizing flux, and then the sink in that cluster with the highest ionizing flux -we do this at the same snapshots as in section 3.2, just before the first supernova for the inner arm, outer arm, and inter-arm models (4.5 Myr after feedback starts), and at 3 Myr for the bar region. We define the position of this sink as the origin of a sphere of radius and locate particles at the surface, for which we calculate the instantaneous mass flux through the surface (radially outward and radially inward). To calculate this numerically, we discretize the surface using HEALPix (Górski et al. 2005), creating equal-area cells (Δ ) at a defined radius. We locate particles between and ( − 1 pc), and sort them according to HEALPix cell. In each cell, we calculate the mean value of the mass flux per unit area , where is mass volume density and = v ·r, i.e. the velocity (relative to the origin sink) in the direction pointing radially away from the origin sink. We then integrate this value over cells to calculate the mass flux in units of M ⊙ yr −1 , separating out the positive (outward) components and the negative (inward) components. We do this for two radii, = 10 pc and 30 pc representing small and large scales respectively. For reference, if gas moves radially outward from the origin at the ionized sound speed (∼ 10 km s −1 ), it will take 1 Myr to reach 10 pc and 3 Myr to reach 30 pc. We track the origin sink backwards through time and repeat the procedure every 0.047 Myr. In Fig. 9, we plot the total mass inside the sphere as a function of time. We also plot the mass flux flowing radially outward through the surface, + , and inward, − . For context, Fig. 10 shows column densities at the time highlighted by the blue dashed line of Fig. 9 -this shows the differences in gas morphology across the region at 1.5 Myr, enough time for radiation to propagate and ionize gas, and for gas to reach a surface.
While feedback is occurring, the inner arm and outer arm models shows the most significant change in total mass inside 10 pc (top row, leftmost panel). In the inner arm, the mass increases by more than a factor of 2, from 0.5 to 1.2 × 10 5 M ⊙ at 2.3 Myr. The streams of dense gas are more collimated compared with the other regions where dense gas is spread out more evenly over the surfaces. After the peak mass is reached, the total mass then declines slowly for more than 3 Myr as gas is gradually dispersed from the system by ionization (while most of the sinks actually fall into the 10 pc sphere over this time).
On the other hand, the outer arm shows a decline in mass over the whole runtime, decreasing by a factor of 2 within the first 2 Myr of feedback. This is partly due to the gas dynamics set by the original galaxy simulation, which shows up as outflow when clouds move across the region, separate to the effect of stellar feedback. As the sinks form and feedback progresses (from this origin sink as well as nearby sinks), the decrease is due to the sinks dispersing in addition to the gas itself. This model is already the lowest density of the four sub-regions (Fig. 10), meaning feedback does not have to be strong to disperse it. The outer arm has low inward fluxes over the whole runtime, but does exhibit large spikes as clumps of material pass through the surface; however, these stop after about 1.3 Myr as much of the material has left by this point -the outward mass flux declines from 10 −2 M ⊙ yr −1 around 1 Myr to 10 −4 M ⊙ yr −1 by 2 Myr.
The inter-arm follows a similar pattern as the inner arm despite the different morphology, albeit over a shorter timespan. It increases its mass by a factor of 2 by the time the feedback starts, at which point this plateaus and the total mass slowly declines (c.f. the inner arm where this takes 2 Myr). This region has high outward mass fluxes (top row, middle panel), with values often exceeding 10 −1 M ⊙ yr −1 and becomes the model with the second-highest flux (behind the bar). The inward mass flux (right panel) drops by almost two orders of magnitude around 1.5 Myr, meaning this region disperses material effectively while not maintaining any infall -the infall here then matches the outer arm, reaching the smallest values of a few 10 −4 M ⊙ yr −1 .
The bar shows the smallest change to the total mass inside 10 pc, simply increasing its mass slowly over the course of the simulation. The material in this model rotates around, causing mass to enter the surface and a comparable amount of mass to leave it, with high mass fluxes of 1-10 M ⊙ yr −1 being reached in this model. Indeed, in the first 1.5 Myr at 10 pc, the bar has the highest inward mass fluxes, followed by the inner arm, inter-arm, then outer arm.
At 30 pc (bottom row), the inward fluxes are generally similar for all models except the bar, which is 1-2 orders of magnitude higher than the other models, especially in the first 1.5 Myr. The other models are all similar to each other, with any morphological differences averaging out over this larger radius. On the other hand, the outward mass fluxes differ more strongly between the models, with the outer arm being an order of magnitude greater than the inner arm and interarm in the first Myr. The fluxes for the latter two slowly increase to match the outer arm by 2-2.5 Myr, as the dispersive effect of feedback increases in these regions (which are initially denser than the outer arm). Overall, the models differ from each other more noticeably when considering the local gas flows (10 pc), with less difference between regions at larger distances (30 pc); the exception is the bar, which has higher fluxes than the other regions even at the larger radius.

DISCUSSION
There have now been multiple studies showing the variation of giant molecular cloud (GMC) properties with environment, initially with individual galaxies (Colombo et al. 2014;Pan & Kuno 2017), but now across larger samples with the PHANGS survey (Sun et al. 2018). The variation of cluster properties is less established, but there is some observational evidence of variation with environment. Messa et al. (2018) find that the cluster mass function is steeper in the interarm region of M51, whilst higher mass clusters are present in the inner parts of M83 (Adamo et al. 2015;Della Bruna et al. 2022). In our models, we also find that the more massive clusters form in the inner regions compared to outer regions, and also in spiral arms compared to inter-arm regions. Our cluster masses also vary with environment similarly to the Milky Way.
We see that the bar region produces the most massive and dense clusters, and the star formation rate is highest here. Observationally, there are often particularly massive clusters at the ends of bars. For example in the Milky Way, the possible super-star cluster W43 appears to be forming at the end of the bar (Nguyễn Luong et al. 2011;Carlhoff et al. 2013). High density gas can be collected together at the bar ends (as seen in the galaxy NGC 3627; Beuther et al. 2017), and strong tidal fields and turbulence can create higher density clouds in the Galactic Centre generally compared to the disc (Oka et al. 2001;Henshaw et al. 2016;Kruĳssen et al. 2019) -these extreme initial conditions could produce high mass clusters. Massive clusters are observed along the dust lanes of the disc-shaped region surrounding the bar in NGC 1365 (Elmegreen et al. 2009;Schinnerer et al. 2023). This is similar in morphology (although larger) to the ring structure produced by the bar in our galaxy simulation. This seems to contrast with some surveys (e.g. Sheth et al. 2000;Momose et al. 2010), which find that there is a lower star formation rate in the bar. However, if we consider the wider bar region of the original galaxy simulation (Fig. 1, first panel), we see that there are also large areas with very little gas. So the wider bar region contains both low density areas and very dense ring or disc-like structures (see also Renaud et al. 2015;Shimizu et al. 2019;Querejeta et al. 2021;Iles et al. 2022;Maeda et al. 2023).
Generally, the trends we see with environment are equivalently trends with initial density, but it is not necessarily possible to readily distinguish between the two. For example, as we suggest above, the particular dynamics of the bar lead to very high densities, and similarly the spiral arms are regions of more strongly convergent flows. We do see some indication of a difference in star formation rate surface density between the arm and inter-arm region for a given gas surface density. The spiral arms may gather gas together, forming large complexes with increased star formation rates compared to other clouds , though it is not clear whether the arms make a large difference to the galactic star formation rate (Dobbs et al. 2011;Eden et al. 2013Eden et al. , 2015Pettitt et al. 2020b;Urquhart et al. 2021).
We find that our star formation rates are high compared with observed extragalactic regions (Fig. 4). One explanation for the discrepancy between our results and the Bigiel et al. (2008) observations is the latter have coarser spatial resolution -see also appendix A and the discussion in Heiderman et al. (2010), comparing the extragalactic results to Milky Way regions and discussing the effect of pixel size. This means regions where little or no star formation is taking place is also included in the observed SFR measure, driving down the spatially averaged Σ SFR . The simulated boxes, meanwhile, encompass the star forming regions without including too much of the non-star forming material -and indeed, there is a selection bias to the measurements as the initial conditions were chosen in part for their potential to form stars. This does not totally explain the discrepancy, however, as our star formation rates/efficiencies are still higher than most of the 100 pc resolution data from e.g. Chevance et al. (2020), Pessa et al. (2021), and Kim et al. (2022). Given the computational difficulty of resolving star formation self-consistently on these scales, it is necessary to use a sub-grid model instead. Together, our tests in appendix A imply the input parameter for the star formation efficiency per cluster-sink is too high (currently 50 per cent) and that the sink accretion radii may need reducing to avoid too much material being accreted. It is also possible that our star formation method is triggered at lower densities than they should be, as our simulation results are slightly shifted to the left (lower gas density) compared to resolved Milky Way clouds. Our clouds can still be compared reliably to each other as they have the same simulation parameters, and in terms of observations, represent clouds with active, high star formation.
Measurements of the rotation of young clusters are rare (Kuhn et al. 2019), but our results predict that rotation is highest in the bar and inner arm, which could simply reflect the higher angular velocity of the galaxy at smaller radii. Finally, our simulations also suggest that differences in the gas flux from the wider environment (>30 pc) are minor -however, there are larger differences between gas inflows/outflows to/from clusters on 10 pc scales. The exception is the bar area, with gas inflow over larger scales in this region, again likely because of the high rotation of the bar compared to the arms. The inflow rates for the bar region are comparable to the values observed flowing along the dust lanes towards the Central Molecular Zone of the Milky Way (Sormani & Barnes 2019) and the nuclear ring of NGC 1097 (Sormani et al. 2023).

SUMMARY AND CONCLUSIONS
We present zoom-in simulations of cloud complexes extracted from a galaxy evolution model similar to the Milky Way, which contains a bar and four spiral arms (Pettitt et al. 2020a). Clouds have been taken from the bar, inner and outer spiral arms, and an inter-arm region, with masses 2 × 10 6 M ⊙ and sizes 100-300 pc. The zoom-ins include ray-traced photoionization from cluster-sink particles (Bending et al. 2020). The new resolution is 0.4 M ⊙ per particle, compared to 600 M ⊙ per particle in the original galaxy run. We have calculated star formation measures and cluster properties as a function of galactic environment. Clusters have been identified with HDBSCAN (Campello et al. 2013) and the degree of clustering measured with INDICATE (Buckner et al. 2019). Our key results are: (i) Denser regions form stars at a higher rate, following the relation Σ SFR ∝ Σ 1.3 gas , which is consistent with the Kennicutt-Schmidt index (Kennicutt et al. 2007). However, the inter-arm model forms stars less efficiently than the spiral arm regions for the same Σ gas , as Σ SFR is a factor of 2-3 below the arms. The bar is always the most star-forming model.
(ii) Almost all the clusters in the bar and inner arm are smaller than 5 pc. Half the clusters in the outer arm and a third in the inter-arm are larger than 5 pc, with radii more similar to associations.
(iii) Similarly, applying INDICATE shows that the degree of clustering is highest in the bar and decreases sequentially down to the inter-arm.
(iv) The bar and inner arm regions are able to form faster rotating clusters, while the outer arm and inter-arm regions tend to produce slower rotators on average. The representative angular velocities /Myr −1 = 0.57 (bar), 0.54 (inner arm), 0.18 (outer arm), and 0.22 (inter-arm).
(v) The dispersive effect of feedback is shown through radially outward mass fluxes measured at spherical surfaces around the most ionizing cluster. Gas streams away from massive stars, and dense clumps show up as bursts of high fluxes. Radially inward fluxes can still be maintained for the bar and inner arm. Regions differ from each other the most at the smaller scale (10 pc), whereas they are more similar at the larger scale (30 pc).
These models do not include stellar winds. In previous zoom-in models of a different galaxy, we have shown this can affect the sink and cluster properties, typically producing smaller clusters (Ali et al. 2022). However, winds only have a minor effect on the gas dynamics and morphology compared to photoionization (see also Gatto et al. 2017;Rathjen et al. 2021). Finally, we do not discuss the role of supernovae as the bar region has not yet evolved to the point of the first supernova, making a comparison between regions difficult. We expect to investigate the role of supernovae in different galactic environments in a future paper.

APPENDIX A: EFFECT OF PIXEL SIZE, ACCRETION RADIUS, AND SFE ON THE KENNICUTT-SCHMIDT RELATION
We recalculate Σ SFR and Σ gas in a fixed area with dimensions = = 100 pc, instead of recalculating the actual area each time as we did for Fig. 4. We repeat this for 300 pc. We calculate the mass of sinks and neutral gas in this area around the centre of mass, and repeat the procedure in section 3.1 (with the exception that here we average the SFR over a longer Δ ≈ 0.094 Myr, to reduce the number of points in the plot for clarity). The results are shown in Fig. A1 with open circles for the smaller area and filled circles for the larger area. Points for the smaller area are shifted towards the top right compared to the larger area, and are close to the Milky Way results of Heiderman et al. (2010). There is also more scatter, as sinks move in and out of the area between time steps, and this changes Δ sinks in equation (2). For the 300 pc area, the bar results now have systematically higher Σ SFR for the same Σ gas than the other regions (whereas for the 100 pc area, and in Fig. 4, they have higher Σ gas too) -i.e. the bar points are vertically higher than the arms and inter-arm, not horizontally. This is because the bar is smaller than 300 pc and hence is not resolved with this pixel size, whereas it is resolved with the 100 pc pixel.
We also test the effect of sink accretion radius, acc , in different runs of the spiral arm region from Bending et al. (2020). One run   Fig. 4, but now the models show the spiral arm region from Bending et al. (2020 with different sink accretion radii ( acc ) and different efficiencies for converting sinks to stars (SFE). has acc = 0.78 pc and one has acc = 0.1 pc. The results are shown in Fig. A2 with green circles showing the larger acc and purple circles the smaller acc . Smaller acc results in lower SFR and brings the depletion time scales into agreement with the Milky Way regions. As before, the pixel size determines the diagonal position.
Lastly, we test the star formation efficiency (SFE) imposed on sink particles, i.e. the proportion of each cluster-sink which is available for conversion to stars as described in section 2.2. The black circles in Fig. A2 have SFE=1 and can be compared with the green points where SFE=0.5. The latter points again have lower SFR.
Combined, these tests show that the K-S relation is sensitive to pixel size (diagonal position), SFE per sink (lower SFE gives lower Σ SFR ), and sink accretion radius (lower acc gives lower Σ SFR ). This paper has been typeset from a T E X/L A T E X file prepared by the author.