Modelling of the Effects of Stellar Feedback during Star Cluster Formation Using a Hybrid Gas and N-Body Method

Understanding the formation of stellar clusters requires following the interplay between gas and newly formed stars accurately. We therefore couple the magnetohydrodynamics code FLASH to the N-body code ph4 and the stellar evolution code SeBa using the Astrophysical Multipurpose Software Environment (AMUSE) to model stellar dynamics, evolution, and collisional N-body dynamics and the formation of binary and higher-order multiple systems, while implementing stellar feedback in the form of radiation, stellar winds and supernovae in FLASH. We use this novel numerical method to simulate the formation and early evolution of open clusters of $\sim 1000$ stars formed from clouds with a mass range of $10^3$-$10^5$ M$_\odot$. Analyzing the effects of stellar feedback on the gas and stars of the natal cluster, we find that our clusters are resilient to disruption, even in the presence of intense feedback. This can even slightly increase the amount of dense, Jeans unstable gas by sweeping up shells; thus, a stellar wind strong enough to trap its own H II region shows modest triggering of star formation. Our clusters are born moderately mass segregated, an effect enhanced by feedback, and retained after the ejection of their natal gas, in agreement with observations.


INTRODUCTION
Star cluster formation is a nonlinear, attenuated, feedback problem: initially cold and dense gas starts forming stars, and the more dense gas is available, the more vigorous the star formation becomes (Mac Low & Klessen 2004). However as the number of stars increase, so do the chances of producing multiple OB stars that can prevent further star formation by injection of kinetic energy, momentum, and radiation in surrounding regions extending out to parsecs. This feedback can not only make the gas Jeans stable, thereby halting star formation, but can also eject the gas from the natal cluster altogether. If the gas dominates the gravitational potential of the cluster at the time of ejection, the cluster may be disrupted by the gas expulsion entirely (Baumgardt & Kroupa 2007), perhaps providing an explanation for the 90% destruction rate of all young clusters found by Lada & Lada (2003).
Stellar feedback in the context of entire clusters has been studied by several authors. Early work included radiation (Peters et al. 2010b;Krumholz et al. 2011) and winds (Krumholz et al. 2012), as well as studying the effects of clustered supernovae (Joung & Mac Low 2006). In a seminal series of papers Dale and coauthors used a smooth particle hydrodynamics code (Monaghan 1992) to study radiation (Dale et al. 2012a(Dale et al. , 2013a and winds (Dale et al. 2013b) independently as well as in combination (Dale et al. 2014(Dale et al. , 2015a for their effects on natal gas clouds. Subsequent studies have further investigated feedback in the form of radiation (Rosen et al. 2016;Peters et al. 2017), winds (Gatto et al. 2017) and supernovae (Kim & Ostriker 2015;Simpson et al. 2015;Ibáñez-Mejía et al. 2016;Girichidis et al. 2016) with increasing accuracy. However none of these works has included at the same time ionizing radiation, winds, supernovae, stellar evolution, and collisional N-body dy-namics capable of following the formation and dynamical evolution of wide binaries and multiple systems.
Our goal is to predict the initial conditions of a newly born star cluster: A cluster formed star by star from magnetized gas and remaining gravitationally bound during the gas expulsion process driven by stellar feedback from evolving massive stars. The gravitationally bound state of young star clusters has been supported by both observations (Tobin et al. 2009;Karnath et al. 2019) and previous simulations of young clusters (Offner et al. 2009), although more recent studies have suggested many young stellar groups may actually be supervirial or unbound (Gouliermis 2018;Kuhn et al. 2019). We are cautious with regard to observational claims of supervirial ratios, though, as mass segregation has been shown to bias virial measures towards smaller group masses and more unbound configurations (Fleck et al. 2006).
Our simulations bridge the gap between gasdominated, star-formation simulations and gas-free, Nbody, star cluster simulations. In a previous paper (Wall et al. 2019, hereafter Paper I), we explained how we use the AMUSE framework (Portegies  to couple the adaptive mesh refinement magnetohydrodynamics (MHD) code FLASH (Fryxell et al. 2000), the N-body code ph4 (McMillan et al. 2012), the stellar evolution code SeBa (Portegies Zwart & Verbunt 1996), and a treatment of tight multiple systems (usingmultiples Portegies ). In the current study we focus on describing the numerical methods developed to implement stellar feedback of the stars acting on gas, and their consequences. The structure of this study is as follows: In Sect. 2 we describe our particular implementations of radiation, stellar winds and supernovae. Further we discuss our modifications to FLASH for far ultraviolet and cosmic-ray background heating as well as our atomic, molecular and dust cooling approach. In Sect. 3 we describe four simulations conducted with our code, and use them in Sect. 4 to examine the effects of feedback on star formation, cluster structure, and the possibility of triggering star formation through stellar feedback. Finally we close in Sect. 5 with a summary.
The source code for our method, including our new and revised routines for FLASH and the bridge script to couple FLASH and AMUSE, are available at https: //bitbucket.org/torch-sf/torch/. Documentation available is summarized at https://torch-sf.bitbucket.io/. We invite community use of this method and participation in its further development.

FEEDBACK IMPLEMENTATION
In Paper I we described our star formation method. In short, sink particles (Bate et al. 1995;Krumholz et al. 2004;Federrath et al. 2010) form in regions of dense, gravitationally bound gas. As soon as a sink forms, a list of stars is drawn by Poisson sampling from a Kroupa (2001) initial mass function (IMF) with a minimum and maximum mass, using the same method as Sormani et al. (2017). As the sink accumulates enough mass to form the next star on the list, that star is immediately formed with position and velocity chosen based on distributions around the sink values (in the initial models described here we used Gaussians with width given by the local sound speed and the sink radius). The sink then moves to the nearest local density maximum and continues accreting.
The minimum mass is chosen in the models presented to be the hydrogen-burning limit of 0.08M . The maximum mass of a star is correlated with the mass of the final cluster, a result found in observations and parameterized by Weidner et al. (2010). We preserve this correlation by choosing the maximum mass that a star can obtain using their integrated galactic IMF. We calculate the maximum mass of a star assuming a given star formation efficiency sfe for conversion of gas into stars, taken to be unity in our present work, from our initial clouds of mass 10 3 M to 10 5 M . This gives us a maximum mass for our 10 3 M runs of ∼ 30 M and a maximum mass for the 10 5 M runs of ∼ 110 M .

Photoionization
For radiation transport we use the FLASH module FERVENT (Baczynski et al. 2015). This module follows the ray tracing algorithm implemented in ENZO by Wise & Abel (2011). The method creates rays from point sources along directions defined using the HEALPIX (Górski et al. 2005) tiling of a sphere, then traces these rays through the block-structured, adaptively refined grid. The number of rays from a specific source hitting each individual block (usually 8 3 or 16 3 cells) is kept constant by splitting rays when necessary. As each ray intersects a cell the number of photons is reduced by absorption while the gas in the cell is ionized and heated accordingly.
The FERVENT package calculates the ionization fraction due to ultraviolet (UV) radiation (Baczynski et al. 2015) Here x n is the fraction of species n, C cl (T ) is the collisional ionization rate, α B (T ) is the case B recombination coefficient, n H is the number density of neutral hydrogen, n e is the number density of electrons and is the rate of photon ionization, specifically formulated to be photon conservative (Baczynski et al. 2015). Here N γ is the number of photons, V is the volume of a cell, x H + is the hydrogen ionization fraction, shortened hereafter to x, and δt is an ionization time step. Equation (1) was originally solved explicitly. Since we only follow hydrogen ionization, we have modified this method of FERVENT to implicitly solve the ionization evolution, which allows for a solution with a much longer time step. First we rewrite Equation 1 with substitutions for the fractional ionization x everywhere then approximate it with a forward finite difference equation which is quadratic in the ionization x 1 at time t + δt, leading to an algebraic solution for x 1 . The error in this method is given by the next term in the Taylor expansion of the method, which can be used to derive a more accurate estimate of the time step than the original method: where c is a tunable safety parameter that we usually set to 0.8. (Our solutions do not seem sensitive to the exact value of c.) Since ionization depends strongly on the radiation field and the temperature, and rapidly converges once these fields find steady states, integration of the ionization differential equation can capture the proper timescale for each of these events. Basing the time step on the rate of change of ionization allows us to do this, while the implicit solution guarantees stability. We have implemented subcycling on the ray tracing, ionization and heating/cooling within a single MHD time step to allow the gas dynamical time steps to be as large as possible. This results in an overall speedup of FLASH of about an order of magnitude compared to calling the (expensive) gravity and MHD solvers during the short ionization time steps.

Radiation Sources
To calculate stellar radiation fluxes (and winds) based on mass we only consider stars with M >7 M . These stars are evolved using the stellar evolution code SeBa (Portegies Zwart & Verbunt 1996, with updates from Toonen et al. 2016, who included triple evolution), which is integrated into the AMUSE framework. The luminosity N γ and average energy ν γ of ionizing photons from these stars is calculated from their mass and age-dependent surface temperature interpolation of the OSTAR2002 grid from (Lanz & Hubeny 2003) if the star has a surface temperature T * >27.5 × 10 3 K; or else just an estimate from integrating the blackbody emission B(λ) as a function of wavelength λ for the star if T * <27.5 × 10 3 K (e.g. Stahler & Palla 2008). Then the cross section σ for the photons is calculated from ν γ (Osterbrock & Ferland 2006) with σ 0 = 6.304 × 10 −18 cm 2 (Draine 2011a).

Photoelectric Effect from FUV
As a second energy bin for radiation we also include far ultraviolet (FUV) radiation (5.6 eV to 13.6 eV) which is absorbed by dust, ejecting photoelectrons in the process. Especially for lower mass stars (7 ≤ M/M ≤ 13), the power in this radiation bin approaches that of photoionizing radiation. Since the cross section of photoelectric photons is smaller than those that ionize hydrogen, these rays penetrate farther into the gas, heating it farther from the source star than the ionizing radiation.
Although only about one in ten FUV photons ejects an electron from the dust, all impart momentum to the gas that can be important in clearing dense gas around newly formed massive stars. Indeed, at our typical numerical resolution, this radiation pressure is the main process acting to clear gas out of zones with densities n H 10 6 cm −3 surrounding early B and late O type stars with weak stellar winds. If this process is not implemented, the ionizing radiation from the star is trapped in the zone, producing an unphysical ultracompact H ii region.
We limit dust to gas with temperatures less than T sputter = 3 × 10 6 K to ensure that gas does not cool unphysically in regions where dust would have previously been destroyed. This is an estimate based on the temperatures at which dust would sputter within the period that we model (Draine 2011a, eqs. 25.13 and 25.14).
To compute the attenuation of radiation in the FUV with luminosity N γ , we follow the method of the original FERVENT paper (Baczynski et al. 2015), calculating the optical depth τ d = n H ∆rσ d as a function of the path length of the ray ∆r and the dust cross-section σ d = 10 −21 cm −2 H −1 (e.g. Draine 2011b). The attenuated radiation N d = N γ (1 − exp (τ d )), and the flux through a cell is then where the Habing (1968) flux G 0 = 1.6 × 10 −3 erg cm −2 .
We then add the momentum E γ /c from the photons absorbed by each cell to the gas, while we follow Weingartner & Draine (2001) to calculate the heating from the photoelectric electrons ejected into the gas as detailed in Sect. 2.4.
Finally, we also allow for EUV photons to be absorbed by dust. Since the cross section for EUV photons is so much larger for hydrogen (Osterbrock & Ferland 2006, σ H ∼ 6 × 10 −18 cm 2 ) compared to that of dust (Draine 2003, σ d ∼ 1 × 10 −21 cm 2 ), we first compute how many UV photons are absorbed by the gas, then any remaining photons are subject to absorption by dust. As a result, most dust absorption occurs inside H ii regions where the gas is completely ionized, leaving the only the dust optically thick to radiation. Although the dust will eventually be destroyed by sputtering, this occurs on timescales long compared to the expansion time of H ii regions (Arthur et al. 2004;Draine 2011b). Similar to previous studies (Draine 2011b;Kim et al. 2016) we find that dust absorption leads to density gradients within our H ii regions.

Supernovae
We include the possibility of explosions from Type II supernovae from massive stars formed in our molecular clouds, as well as Type Ia supernovae from white dwarfs in the field. Supernovae were previously implemented in FLASH by pure thermal energy injection to study the large-scale driving of turbulence (Joung & Mac Low 2006). However recently several authors have demonstrated that more accurate results can be achieved using injection of both kinetic and thermal energy in a mixture that depends on numerical resolution (Simpson et al. 2015;Kim & Ostriker 2015;Gatto et al. 2015). Simpson et al. (2015) derived an analytic expression for the kinetic fraction f kin based on how well a simulation resolves the pressure-driven snowplow (PDS) as where ∆x is the width of a grid cell, µ is the mean molecular weight, n o is the background number density, E 51 is the supernova energy in units of 10 51 erg, and R PDS and t PDS are the radius and time of the supernova transition into the PDS (Draine 2011a). We have implemented the method of Simpson et al. (2015) for supernova injection into our version of FLASH: cloud-in-cell (CIC) linear interpolation is used to map the energy input into the grid from the supernova onto a 3 3 cube centered at the supernova location. Thermal energy and mass are equally divided among the 27 cells. Kinetic energy is also equally divided and injected in the form of momentum into all but the center cell, where instead the kinetic energy is converted into thermal energy and added to the thermal energy already present. The contents of each zone in the cube are then mapped onto grid zones they overlap (Simpson et al. 2015, Fig. 1). Initial testing shows that even at low resolution, the supernova remnants are nearly spherical and exhibit the proper transitions from Sedov-Taylor to PDS (Draine 2011a), as shown in Figure 1.
An energy plot for the same run is shown in Fig. 2. Here we note that the transitions between the Sedov-Taylor solution (Draine 2011a, when the kinetic energy fraction ∼ 0.25), transition (Haid et al. 2016), and PDS (Draine 2011a) all appear to match the analytic solutions very well. Also the slope of E(t) ∝ t −3/4 during the PDS is well recovered.

Winds
In recent years, the general strategy for injecting winds has been to inject mass and velocity in a region around the star to set the kinetic energy of the wind over the timestep (Pelupessy & Portegies Zwart 2012;Gatto et al. 2017;Rimoldi et al. 2016). This is also the method of Simpson et al. (2015) for supernovae, who use this energy to calculate the momentum for each cell. However, adding momentum to a grid cell that already has mass does not add the same amount of kinetic energy as adding the momentum to an empty cell. Therefore, after adding the momentum of the wind to each cell in the source region, we compute the resulting kinetic energy, and conserve total wind energy by injecting the missing energy as thermal energy into the cell. Stellar wind feedback is implemented using a method of momentum injection of our own design which was inspired by inverting the method of Simpson et al. (2015). The amount of energy deposited by stellar winds onto the grid is given by the mechanical luminosity whereṀ is the stellar mass loss rate, typically 10 −8 M yr −1 to 10 −6 M yr −1 for O and B stars, and v w is the wind terminal velocity, typically 3 × 10 2 km s −1 to 3 × 10 3 km s −1 for the same stars. We consider the update over time ∆t of an individual cell with volume V cell , density ρ old , specific internal energy e intold , and velocity v old . We first calculate the overlap fraction φ between the spherical wind injection region and the cell itself using a 20 3 subgrid of sample points in each cell, following a routine by D. Clarke included in ZEUS-MP (Hayes et al. 2006). This value is normalized to the full volume of the source region (i.e. vol φ = 1). The change in density of the cell The stellar wind input kinetic energy for this cell is The final velocity of the cell can be computed from momentum conservation to be The final change in specific kinetic energy is then so the specific internal energy of the cell needs to be increased to In determining the radius (e.g. the number of cells) across which to inject the winds, both the physical radius of the wind and the ability of the Cartesian grid to resolve the spherical input are important. Here the analytic solution for a stellar wind bubble is our guide. The radius of the wind termination shock (Weaver et al. 1977, Eq. [12]) where t w is the lifetime of the wind and ρ 0 is the background density. Within this region the wind will be free streaming. If this radius is resolved by more than a single cell, we directly inject the energy and momentum calculated as above in the resolved region out to a maximum radius of 6 √ 3 ∆x, a radius at which we find spherical winds to be well resolved. If R 1 < ∆x, we set the radius of the injection region to be ∆x. Note that since the stars are Lagrangian particles not restricted to the cell centers, even wind bubbles smaller than a single cell generally inject not just thermal energy but also momentum and kinetic energy into the grid by straddling cell boundaries.
For each cell within this radius we determine the fractional overlap of the cell with the wind injection region and add that fraction of momentum and thermal energy into the cell, with the momentum and energy evenly distributed throughout the sphere defined by the injection radius. To guarantee that all cells are at maximum resolution in the source region, we add a new criterion that enforces refinement of all blocks that lie within the injection radius of the star.
The mass of hot gas within real stellar wind bubbles is determined by conductive evaporation (Weaver et al. 1977), as well as turbulent mixing, across the contact discontinuity at R C (see Fig. 3) between the hot, rarefied, shocked stellar wind and the dense, radiatively cooled shell. The extra density injected by these mechanisms reduces the temperature in the hot region between R 1 and R C , and thus the sound speed. Capturing this physics exactly is computationally challenging, as conductive evaporation is a diffusive process with Courant timestep ∆t diff ∝ ∆x 2 . However, the Courant timestep ∆t C = C CFL ∆x/max(v, c s ) is dominated by the high sound speed c s in the hot region. Therefore, we introduce the option to mass load the wind to bring its temperature to the correct order of magnitude. We choose a mass-loaded temperature target T ml =5 × 10 6 K in our simulations, set to lie at the low end of the quasi-stable hot gas phase (McKee & Ostriker 1977a). We ensure this temperature by reducing the pre-shock velocity of the wind such that the post-shock temperature (Draine 2011a) We make up the lost energy by adding to the mass of the wind until we recover the proper wind luminosity.
We have found this to be sufficiently hot that the wind bubble in diffuse gas (n ∼ 1 cm −3 ) continues to conserve energy in the shocked gas as expected, while allowing significant gains in the size of the Courant time step. The combination of this radius with the division of kinetic and thermal energy described above leads to bubbles in dense gas primarily injected with thermal energy until the free streaming wind is resolved (and with it the inner boundary of the hot wind bubble), at which point the injection of energy shifts over towards kinetic. This method shows excellent agreement with theoretical predictions for wind bubbles given by Weaver et al. (1977). Comparison between their analytic solutions and a plot from our test runs is shown in Figure 3.
To calculate the mass loss rates for our stars we follow Vink et al. (2000) while for the velocities we use the fitting formula of Kudritzki & Puls (2000). Note that these methods include the bi-stability jump in wind strength in late O and early B stars produced by the higher absorption cross section of Fe iii compared to Fe iv (Vink et al. 2000).

Ultraviolet Radiation Heating
Heating from our stellar sources in the EUV and FUV is calculated by converting the photons to energy fluxes, then applying those fluxes weighted by the probability of an electron of a given energy actually being ejected from an ion or dust grain when it absorbs a photon. For hydrogen ionization this is done by simply differencing the energy of the photon and the ionization potential of hydrogen.
For the background FUV we assume a constant flux of 1.7G 0 and estimate the local visual extinction where N H = 1.87 × 10 21 cm 2 and we take the length scale to be the local Jeans length (Jeans 1902) following the methods of Seifried et al. (2011) and Walch et al. (2015). The fraction of FUV radiation that can heat the gas is then f ext = exp(−3.5A V ).
For the FUV flux we normalize to the Habing flux G 0 following Equation (8) to find G, and then calculate the heating function per unit volume where is a heating efficiency function. We have implemented several different approximations to this function, including detailed fits from Weingartner & Draine (2001) and Wolfire et al. (2003), as well as a simple adjustable parameter following Joung & Mac Low (2006).
The coefficients are taken from their Table 2, using the case with a ratio of visual extinction to reddening R V = 3.1, carbon abundance with respect to hydrogen of b c = 6 × 10 −5 , distribution A, which minimizes the amount of C and Si in grains, and the stellar radiation field of a B0 star, corresponding to a blackbody with 3 × 10 4 K up to 13.6 eV. The flux factor The Wolfire et al. (2003) function is given by with φ PAH = 0.5 following their assumption. Finally, the simplest assumption is to just follow Joung & Mac  Figure 1 of Weaver et al. (1977) and our own test withṀ = 10 −6 M yr −1 and v = 2 × 10 3 km s −1 , at time t = 10 4 yr. The inner white line indicates the analytic solution for the stellar wind termination shock radius R1 and the outer black line the solution for the outer shock radius R2 from Weaver et al. (1977). Note that the shell has cooled and collapsed at this point, with the distance between the contact discontinuity Rc and R2 now resolved by at most two cells. Low (2006) and set a constant value of = 6.5 × 10 −26 erg s −1 . In the example models analyzed in Paper I and here, we use the Weingartner & Draine (2001) approximation (Eq. [21]).

Dust Temperature
Any photons absorbed that do not eject electrons contribute directly to heating the dust. The dust density is assumed to be a constant 1% fraction of the gas density (Draine 2011a). When solving for the dust temperature we use the radiative cooling rate from Goldsmith (2001), assuming the dust is always optically thin at our densities n H < 10 6 cm −3 , while applying photoelectric heating as previously described. To calculate the dust temperature we use Newton's root finding algorithm as in Seifried et al. (2011), which generally converges in less than ten steps.

Cosmic Ray Heating
Cosmic ray heating is applied with an ionization rate of ζ = 10 −17 s −1 and a heating rate of Γ cr /n H = (20 eV)ζ = 3 × 10 −27 erg s −1 as appropriate for the dense regions we are attempting to simulate (Galli & Padovani 2015).

Gas Cooling
For gas cooling we include contributions from atomic and molecular species as well as dust grains. For the atomic contribution we use the piecewise power law in Joung & Mac Low (2006, Fig. 1), itself derived from equilibrium ionization values given by Sutherland & Dopita (1993). For molecular cooling we use tabulated values from Neufeld et al. (1995) which were originally implemented in FLASH by Seifried et al. (2011), while for dust we use the method of Goldsmith (2001) with the cooling equation for dust from Hollenbach & Mc-Kee (1989). Note that dust cooling is also limited to temperatures T < T sputter (see Sect. 2.1.3).

Numerical Solution
To solve the implicit difference equation for the temperature of the gas under all of these heating and cooling sources we have implemented Brent's 1973 method, which we find to be more accurate and stable than the Euler method used by Joung & Mac Low (2006) and Baczynski et al. (2015). In each cell, all heating Γ i ( ) and cooling Λ j ( ) rates are combined to find the rate of change of specific internal energy The cooling rate at the minimum allowed temperature in the simulation (generally 10 K, but could be as low as the CMB background temperature) is calculated and subtracted from the total cooling rate to set a tempera- ture floor. Then, the difference equation is solved for e i+1 by the Brent method. Figure 4 shows the temperature for a cell initially at T = 10 5 K cooling over 10 Myr using the original Euler method of Joung & Mac Low (2006) and our implicit method, as well as a more expensive RK 4(5) method with local extrapolation (Press 2007). For a given time step criterion, the implicit method is generally about twice as accurate as the Eulerian method and ∼ 30% faster than the RK 4(5) method. Given that we call this solver on every iteration of the ray tracing method as we converge to an ionization solution, we have chosen to use the implicit method due to its combination of speed and accuracy.

Tests
In order for our simulations to resemble the real ISM as closely as possible, our heating and cooling solutions should be able to replicate the three-phase ISM (McKee & Ostriker 1977b;Cox 2005), where we have equal pressure solutions for a quasi-static hot medium and warm and cold media. Being able to maintain the different thermal phases of the ISM in pressure equilibrium is important generally. It is particularly important for the initial conditions chosen for our proof-of-concept models, since the background medium for some of our clouds is warm neutral medium, while the clouds themselves always consist of cold neutral medium in pressure equilibrium with the background at the initial time.
To test our heating and cooling methods we therefore created a single cell simulation that iterates over hydrogen number densities with 10 −4 ≤ n H ≤ 10 4 cm −3 , solving for the equilibrium temperature and pressure with =2.0E-17 =-0.021 G=1.69 G 0 P/k T Figure 5. Temperature (blue, solid line) and pressure (red, dashed line) of a single gas cell with 10 −4 ≤ nH ≤ 10 4 cm −3 evolved until equilibrium in the presence of a background UV flux of G = 1.69 G0 and cosmic ray background ionization rate of ζ = 2 × 10 −17 s −1 integrated using the implicit method. The two-phase medium occurs for densities of roughly 10 −1 ≤ nH < 10 3 cm −3 , and a quasi-static third phase at high temperature and low density is then consistent with the pressure.
Milky Way-like background FUV and cosmic ray values. The results are shown in Figure 5, where the threephase medium is shown to be stable in our simulations from P ∼ 4 × 10 3 K cm −3 to 2 × 10 4 K cm −3 , reproducing similar ranges found in Wolfire et al. (2003, Fig . 7) for the solar neighborhood. Note that we can adjust this range by increasing or decreasing our background FUV, as discussed in Hill et al. (2018), which also uses the atomic cooling model our method is based on.

EXAMPLE RUNS
For testing stellar feedback we present four proof-ofconcept runs, three of which include radiation, winds and supernova. However, we terminated these runs for cost reasons before any massive star had exploded as a supernova, so we restrict our discussion to radiation and winds. In Table 1 we show the initial cloud properties, grid resolution and time of initial star formation for these runs.
All four simulations use an initial density field that is spherically symmetric and distributed radially as a Gaussian (Bate et al. 1995;Goodwin et al. 2004) with full width at half maximum of the cloud radius R. The velocity field is initialized with a turbulent Kolmogorov velocity spectrum v(k) = v 0 k −5/3 from wavenumber k = 2 to k = 32 (Wünsch 2015) for the dense gas, where k = 2π/D. We note that choosing the initial conditions does determine a good deal about the subsequent evolution (Goodwin et al. 2004;Girichidis et al. 2011). The surrounding medium is initialized with zero velocity and Table 1. Parameters for each of the four runs described here including cloud mass M , radius R, and central density ρc, in units of ρ = 2.39 × 10 −23 g cm −3 , normalization of the velocity perturbation spectrum v0, the number of refinement levels N ref , cell size ∆x at maximum refinement, and the domain size D. We further give the total number of stars Ns at end of run at time t end when analysis was performed, as well as the time the first star formed t sf . Note that M3 and M3f used different random turbulent patterns initially, explaining their different values of t sf . a Runs ending in "f" include feedback due to radiation and stellar winds.
in pressure equilibrium with the cloud gas. Refinement and derefinement are controlled by the Jeans (1902) criterion as described in Federrath et al. (2010). We now describe individual characteristics of these runs.

M3
In this run, which does not include feedback, two subclusters form that subsequently merge. We highlight the merger event in the relevant plots throughout by including a grey shaded box on each plot covering the time of the merger event. This is our only run among our proof-of-concept models that features a merger.

M3f
In this run several stars form in the main cluster that are massive enough to have ionizing and wind feedback. The growth of an H ii region is initially suppressed by dense gas accretion, creating flickering H ii regions (Peters et al. 2010a,c;De Pree et al. 2014). This continues until the number of massive stars gets large enough for their combined wind and ionization feedback to clear an expanding H ii region around the main cluster, near the end of the run. Around this same time a second cluster forms in the simulation. While the two clusters appear to be falling toward each other, they have yet to merge at the end of the run.

M3f2
This run starts with the surrounding lower density gas in the warm, neutral phase at 4 × 10 3 K, as opposed to the cold phase at roughly 60 K in the other two 10 3 M runs. Therefore, the surrounding gas is less dense than those runs, resulting in slower accretion onto the star forming region. Similarly, feedback is more effective in expelling the gas from the star forming region, since the feedback sees a smaller surface density above it (Grudić et al. 2018). The more effective feedback rapidly shuts down star formation in this run, which therefore only produces 52 stars.

M5f
In our final run we start with an initial sphere of 10 5 M and a radius of 50 pc. The gas outside the sphere is initially in the warm ionized phase at ∼ 8 × 10 3 K. Once the gas collapses and begins to form stars, a large central cluster appears, as well as a secondary, smaller cluster. The central cluster rapidly grows until it stochastically forms a 97 M star, the most massive star formed to date in any simulation we have run. This star rapidly expels all remaining gas from the central cluster, leaving pillars of gas surrounding the star forming region (see Fig. 7 d) that resemble the Eagle Nebula and similar formations (e.g. Hester et al. 1996;McCaughrean & Andersen 2002;McLeod et al. 2015).
A difficulty in performing simulations of large clouds from idealized initial density conditions stems from the initial free-fall time for the gas in the Gaussian sphere, which for this run is only 8.6 Myr. Since turbulence decays within a free fall time t ff (Mac Low et al. 1998), the velocity distribution of the gas became quite smooth by the time star formation commenced in this run at 15.4 Myr. Similar concerns were discussed by Krumholz et al. (2012), who noted that this affects both star formation rate and efficiency. Future models of high mass clouds will need to start with more realistic initial conditions that better model the actual assembly of such structures.

Stellar Group Identification
To identify stars as group members within the simulation we used two methods, HOP (Eisenstein & Hut 1998), which is included in AMUSE, and the Scikit-Learn  HOP determines group membership by the following procedure: 1. Calculate the local density at each particle using its N nn nearest neighbors and the local density gradient at each particle using its N hop nearest neighbors.
2. From each particle, hop to the next particle of the N hop neighbors in the direction of the highest density gradient. Continue until the current particle is the density maximum of the N hop nearest particles.
3. Identify this particle as a group core particle, with this particle's density as the group peak density ρ peak . All particles in the path of hops leading to this particle are added to this group as members. Repeat this process until all particles are in groups.
4. Identify particles that reside on the boundary between two groups by identifying particles where one of its N merge nearest neighbors belongs to a different group. Record the density at these locations as the saddle density, ρ saddle , calculated as the average density between the two boundary particles.
5. Merge any groups where ρ saddle is either greater than an absolute saddle density δ saddle , or whose ratio of saddle to peak densities is less than a given relative saddle density factor threshold f saddle . In mergers the group with the lower peak density ρ peak is merged into the group with the higher peak density.
6. Remove any group whose peak density is lower than the outer threshold density δ outer .
For physical parameters in HOP we use an outer stellar mass density threshold δ outer = 1 M pc −3 , an order of magnitude lower than the average stellar density of an open cluster and an order of magnitude greater the stellar density of the Milky Way (Binney & Tremaine 2011).
For the peak density we use δ peak = 3δ outer as suggested in the original paper. For the saddle density we use a relative saddle density threshold, where the boundary saddle density is compared to the minimum peak density of the two groups, defined by D = δ saddle min (δ peak,1 ; δ peak,2 ) .
If D < f saddle , where f saddle is the saddle density threshold factor, the two groups are merged. For our analysis here we set f saddle = 0.5, slightly more aggressively merging groups than the default value of 0.8. The values (N merge , N hop , N nn ) = (4,16,64), again as suggested in Eisenstein & Hut (1998). DBSCAN on the other hand determines group membership using a simpler procedure: 1. Any particle with at least N min neighbors within a distance ξ is considered a core particle.
2. Any particle that is within distance ξ of at least one core particle, but has fewer than N min neighbors, is considered a boundary particle.
3. All connected core and boundary particles define a group.
4. Any other particles are defined as noise.
For DBSCAN we set the N min = 16 for a core particle to be 16 and we calculate ξ, the maximum neighboring particle separation, to match our physical parameters in HOP. Assuming an average stellar mass of M avg = 0.56 M for the initial mass function (IMF) of Kroupa (2001) and using a stellar density of ρ bg = 1 M pc −3 we compute which we used for the one run we analyzed with DBSCAN here, M5f. Generally we prefer HOP due to its physically motivated thresholds, particularly its ability to compare the relative saddle density between two groups to the minimum peak density of the groups themselves to determine if the two groups should be merged. However it was more practical to use the simpler DBSCAN technique for our largest data set. As with any group-finding numerical technique, both methods struggle to disentangle whether one or two groups exist just before the point of merger. However, we only have one run that experiences a merger of two nearly equal sized groups, while others are either well separated or have mergers where one group is clearly the more massive and dominates the potential. As a check, we examined several times during M5f and verified that we found similar results with HOP and DBSCAN. Both methods provide similar and consistent grouping results when compared on the same data over multiple grouping computations, with an agreement of > 95% on cluster members.

Star Formation
The star formation rate (SFR) as a function of time in our simulations is shown in Figure 8. The data shown as blue points is initially calculated by a second-order central difference of the stellar masses every timestep (as little as 100 yr), which are generally quite noisy. We also show the result of using a Savitzky-Golay (1964) filter convolved with a window size of 51 and fit to a thirdorder polynomial to smooth the data before we take the derivative (black lines). We use this filter on all data hereafter presented with open circles, representing the data taken directly from the simulation, accompanied by line plots, which show the result of the smoothing. For the SFR, we further smooth with a Gaussian filter with 3 kyr variance.
In M3 the SFR generally stays high, only briefly decreasing due to a strong interaction in the main group that placed many of the massive stars on wide orbits and slowed the overall accretion rate of the region. In M3f the data are much noisier due to the flickering H ii regions (see Sect. 3.2), which heat the gas briefly and inject turbulence, but never provide enough outward momentum to the nearby gas to eject it nor ionize enough material to prevent it from cooling again. Run M3f2 shows relatively stable star formation until the first massive stars appear at 5.45 Myr. Their feedback breaks up the filament in which stars are forming, but cannot entirely disperse the dense gas in the region. This allows star formation to continue for another 10 4 yr until an interaction between two massive stars expels one far enough out of the center of the group for a second H ii region to form and expand out of the star-forming region. In run M5f the SFR shows large variations as filaments form and intersect in the first megayear and two separate groups form. At around 17 Myr, a 97 M star forms in the more massive group, emitting radiation and winds that travel throughout the simulation domain, although star formation continues both near and far from the star for another ∼ 3 × 10 5 yr before effectively terminating. The amount of gas available to form stars is presented in Figure 9 where we show both the fraction (by mass) of dense gas (for which n H > 10 4 cm −3 , the limit generally considered for gas to be star forming; Lada et al. 2010) and the fraction of Jeans unstable gas. For all the runs except M5f, the amount of Jeans unstable gas is much smaller, generally by a factor of two or more, than the amount of dense gas. This agrees with results found by Dale et al. (2015a), who concluded that dense gas is a necessary but not sufficient condition for star formation. Also in agrement with Dale et al. (2015a) is our finding that effective stellar feedback, which occurs in runs M3f2 and M5f, actually produces more dense gas, rather than reducing it. Only in run M5f does feedback also seem to increase the amount of Jeans unstable gas, thereby leading to increased star formation near the center of feedback. The difference between the feedback in M3f2 and M5f is that the stellar wind from the extremely massive star in M5f can sweep up a shell dense enough to trap its own H ii region, allowing some triggered star formation in the shell (see Sect. 4.4).

Stellar Group Structural Evolution
We next examine the effectiveness, or lack thereof, of stellar feedback on the evolution of the groups that contain massive stars. Given the expectation that 90 % of all clusters are disrupted (Lada & Lada 2003), presumably by gas expulsion, we might expect any feedback that completely ejects the natal gas to destroy the cluster it formed from by removing the dense gas potential helping to bind the cluster (e.g. Tutukov 1978;Elmegreen 1983;Parmentier et al. 2008;Goodwin 2009;Rahner et al. 2017Rahner et al. , 2019.

Energetics
Figures 10 and 11 show the energy in stars and gas respectively contained within the radius of the main group for each run. These figures show that only runs M3f2 and M5f actually eject their natal gas, when their total gas energy becomes positive. However both stellar groups remain bound with negative total energies after the gas is removed, identifying them as true clusters. This is further confirmed by looking at the virial ratios α = 2T /U of the gas and stars in these groups (Fig. 12), where the group as a whole appears briefly unbound during the time that some of the outer stars escape following gas ejection, but the overall cluster survives and returns to a bound virial ratio in both cases. Indeed, the cluster in run M5f is subvirial at the time of the final snapshot, and the other 3 clusters are close to being virialized, regardless of the current state of the gas in the region defined by the cluster. We do see mass segregation in our runs, as detailed below, which might contribute to their being observed as supervirial, something we will examine in more detail in future work. Our clusters appear likely to survive gas ejection, and mass segregation may assist in their survival, since increasing stellar to gas density ratios increases the likelihood of surviving the gas ejection stage (Kruijssen et al. 2012). Figure 13 shows that the mass in gas dominates the mass in all the groups for most of their evolution, only being driven completely out of the group under the intense feedback of M5f's massive star. Even in M3f2, where the gas actually has positive energy, it has yet to be driven from the group entirely. Indeed, dense gas is growing in the group (Fig. 9 c) as it continues to fall in from the filament and build up along the edge of the H ii region. This infall itself may lead to more star formation, although at the end of the run there has been no similar increase in the amount of Jeans unstable gas, and the total number of stars in the group has been steady for the last 5 × 10 4 yr, as shown in Figure 14. This is similar to our other 10 3 M simulation M3f, where feedback near the main group during the previous 10 5 yr has also stabilized the number of stars. In M5f, feedback eventually leads to the loss of some of the least bound stars in the main group as the gas is ejected, shown by the correlation in the drop in gas mass with the decline in the number of stars in the group.

Radius
In Figure 15 we show the Lagrangian radii for evenly spaced mass bins. In M3 (Fig. 15a), lacking feedback, the group radius drops, aside from a brief bounce when two subclusters merge (grey bar). Gas ejection leading to loss of the least bound stars can be seen in the fast rate of growth of the Lagrangian radii of the central groups for the two runs in which feedback expelled significant amounts of gas: M3f2 (Fig. 15c), and M5 (Fig. 15d). In M5f, the 25% Lagrangian radius of the group only grows by ∼33-50%, but the the outer (100%) and half-mass (50%) radii almost double after the onset of stellar feedback.

Mass Segregation
Next we consider the mass segregation of the groups formed in our simulations. Several methods exist to quantify mass segregation, including looking at the halfmass radii R hm of different mass bins in the group (McMillan et al. 2007;McMillan et al. 2015), the mean or median radius of a subset of massive stars (Bonnell & Davies 1998), and calculating the Gini coefficient of the group (Converse & Stahler 2008;Pelupessy & Portegies Zwart 2012). Allison et al. (2009) pointed out several issues with these methods of computing mass segregation, including how binning can affect the results, the reliance on properly finding the group center, and difficulty in comparison to observations. They presented a new method, based on calculations of N random minimum spanning trees (MSTs) of the group and the MST of the N ms most massive stars contained in the same group as a model independent way of determining the amount of mass segregation. They compared the ratios of the norm of lengths of the random trees, l norm , to the length of the massive star tree, l ms , to obtain the mass segregation ratio where Λ msr > 1 indicates mass segregation in the group. This method is independent of any determination of the group center, always returns the same tree lengths (even if the tree is drawn in a different order), and is simple to implement in both two and three dimensions. Further, since the number of random trees calculated provides a standard deviation of tree length, error for the calculations are straight forward to obtain. We show our calculated values of the three-dimensional value of Λ msr in Figure 16 using N ms = {5, 10, 20} and 50 random samples drawn from each group for the comparison trees. We only start tracking groups once they reach 64 stars in size, with the exception of M3f2 where we start following the group at 24 stars.
Two points can be made with this data. First, all of our runs become mass segregated at early times. This presumably occurs because our groups of N stars, with initially short crossing times t cr = R hm /σ v , have likewise short half-mass relaxation times (Binney & Tremaine 2011) t r = 0.1N t cr / ln N . The wide range in stellar masses then accelerates the dynamical evolution of the cluster to a fraction of the half-mass relaxation time scale t seg ∼ ( m / m hm ) t r , where m and m hm are the mean mass of all and of the high mass stars respectively Portegies Zwart & McMillan (2002). The clumpiness of the stellar distribution helps to preserve this primordial mass segregation throughout the assembly of more massive stellar conglomerates, as was predicted by McMillan et al. (2007) from simulations of small merging subgroups.
Second, feedback seems to be correlated with mass segregation in all of the runs including it. We attribute this to gas expulsion having a stronger effect on low mass, loosely bound stars, causing their orbital radii and kinetic energy to increase more than massive stars and leading naturally to an increase in mass segregation even as the whole group expands. Also, all runs that experience significant mass segregation sustain that segregation over their ten most massive stars ( Fig. 16c and  d), even if the five most massive stars have strong interactions that reduce their ratio Λ msr . This supports the view (e.g. Girichidis et al. 2012b,a) that subgroups will start more mass segregated than dynamics alone can account for if they can survive the ejection of their natal gas, as seen in many observations of young stellar groups (Hillenbrand & Hartmann 1998;de Grijs et al. 2002;Gouliermis et al. 2004;Converse & Stahler 2008).

Triggered Star Formation
A modest level of triggered star formation has been seen both observationally (Thompson et al. 2012;Liu et al. 2017) and numerically (Dale et al. 2012b), although some care must be exercised in observations since the time evolution of the system is not available as it is in simulations, making it difficult to disentangle triggered star formation from formation that would have otherwise occurred naturally due to gravitational collapse (Dale et al. 2015b). Dale et al. (2007) and Dale et al. (2015b) divide triggering into two categories; weak triggering where star formation that would already occur due to normal collapse is accelerated, but without increasing either the overall star formation efficiency or the number of stars created; and strong triggering where collapse is induced in previously stable gas that increases the total star formation efficiency, the number of stars, or both. Apparent triggered star formation occurs in run M5f, which has the strongest feedback. At the time of formation of the 97 M star (hereafter referred to as A * ) in the main group, there are three star-forming sinks present. The first sink (sink # 56) is the one that actually forms A * , and its accretion immediately shuts down. The other two sinks (# 55 and 57) continue accreting gas for another 0.3 Myr from gravitationally unstable regions in the swept up shell driven by the stellar wind from A * (Fig. 17). This appears to be a case of triggering maintaining the global star formation rate temporarily (see Fig. 8 d) even in the presence of strong negative feedback.
In the case of sink 55, star formation was already proceeding at a vigorous rate before A * appeared. There-fore this cannot be considered even weakly triggered star formation, however it seems that the intense feedback was unable to do more than briefly slow the production of stars for about 10 4 yr before the sink returned to producing stars at a rate exceeding 4 × 10 −5 M yr −1 . In the case of sink 57, though, star formation was at less than 1 × 10 −5 M yr −1 when the stellar wind bubble compressed gas in which the sink was embedded. Once this occurred, the rate of star formation grew rapidly. For this to be considered triggered star formation, we should also be able to observe an increase in both the amount of dense gas and Jeans unstable gas in the region surrounding the group occurring as the feedback impacts the region. We show that this increase indeed occurs in Figure 9 (d).
During their ejection, both sinks reached peak speeds of ∼ 30 km s −1 with respect to the group center as they followed the expanding gas driven by feedback. This velocity is noteworthy, since it defines the boundary velocity for massive O and B stars that are considered runaway stars (Gies & Bolton 1986). Generally the production of OB runaways has been considered the result of kicks from a binary partner that goes supernova (Blaauw 1961;Portegies Zwart 2000) or due to dynamical interactions with binaries (Leonard & Duncan 1988;Fujii & Portegies Zwart 2011). However since our sinks (and therefore also the gas they accrete) reach velocities comparable to that of OB runaways, triggered star formation in our simulation shows a third, and not previously considered, method for producing OB runaways.
In this case we produced many lower mass stars, but the total mass produced by the two sinks during this time was over 30 M , therefore the lack of formation of an OB star was simply due to the random selection of our star formation method. The high gas velocity is clearly connected to the feedback of A * , but which physical process contributes the most? The radius and velocity of the D-type front are Spitzer (1978) The velocity of the D front has a maximum value of v d ∼ 15 km s −1 , too slow for our gas, ruling out compression by radiation. Note this also likely rules out radiation driven implosion (Sandford et al. 1982) as a primary trigger. We also considered a champagne flow as a possible method of driving the gas velocities, but as shown in Bodenheimer et al. (1979) and similar to the case of radiation driven implosion, champagne flows only accelerate the lower density gas to high velocities. Even in their case (5), where they allowed a D-type front to move past a small dense cloud, the dense gas was compressed but the dense gas velocities never exceeded 8 km s −1 . This leaves the effect of the winds as the main factor accelerating the gas flow. Normally, the wind bubble would evolve while trapped within the H ii region (Weaver et al. 1977). This means for moderately massive O stars the H ii region dominates the dynamics, since the D-type front strikes the ambient gas first (Mc-Kee et al. 1984). However in the case of winds moving rapidly into a region of dense gas the H ii regions can become trapped within the wind shells (van Buren et al. 1990;Mac Low et al. 1991). To calculate the speed of the shell from the wind, we obtained the luminosity and temperature of A * using SeBa and then calculated the wind luminosity using our stellar wind code as described in Sect. 2.3, finding L w = 2.47 × 10 37 erg. The shell velocity of a stellar wind bubble in a uniform medium is (Weaver et al. 1977) with n o = 10 3 cm −3 to find the wind shell velocity at the peak time of the sink velocities, which is ∼ 3 × 10 4 yr after the formation of A * . This gives us a shell velocity of V 2 = 31 km s −1 , consistent with the maximum sink velocity of 30 km s −1 . Eventually all the star forming filaments that are in close proximity to the main group are disrupted by the feedback of A * . At this point (∼ 17.3 Myr) the overall star formation rate rapidly drops, as all gas in the region becomes warm, low density H ii gas or hot wind shocked gas. Some small amount of star formation still occurs in a smaller secondary group containing sink # 22, but formation here is slow due to the overall smaller fraction of dense gas present in the group, as shown in Figure 18 5. SUMMARY In this paper we have described the implementation of stellar feedback methods in FLASH, whose inclusion within the AMUSE software framework was described in Paper I. We have shown that our implementations reproduce standard benchmarks for ionizing radiation, stellar winds and supernovae. We also include heating from cosmic rays and non-ionizing radiation from both individual stars and the galactic background, and radiative cooling from both gas in collisional equilibrium ionization and dust.
These implementations are part of a larger effort to model the formation and early evolution of star clusters, where we combine magnetohydrodynamics using FLASH, collisional N-body dynamics using ph4, binary and higher-order multiple dynamics using multiples, and stellar evolution using SeBa, with the stellar feedback described here. We have begun to use this framework to study newly formed stellar groups and clusters.
We here report on four proof-of-concept simulations with initial gas masses of either 10 3 M or 10 5 M . From these models we conclude the following: 1. Stellar feedback effectively controls the SFR (Sect. 4.1). The details of cloud structure do matter, however, both for the overall star formation rate and because dense shells swept up by feedback can trigger small amounts of additional star formation.
2. Stellar feedback tends to increase the amount of dense gas present in the star forming region, agreeing with Dale et al. (2015a). Contrary to them, however, we do find a case in which feedback even increases the amount of Jeans unstable gas.
3. Our stellar groups generally form subvirial and end marginally virialized (Sect. 4.2.1). Both groups that ejected their gas (M3f and M5f) went through a period of supervirial expansion but ended subvirial, even while they continue to expand.
4. Feedback results in the ejection of gas from our clusters, but did not disrupt any of the stellar groups created in the runs presented here, although the least bound stars were lost (Sects. 4.2.2 and 4.2.3).
5. Our stellar groups quickly become mass segregated (Sect. 4.3). Feedback-driven gas removal further stratifies the stars according to their current binding energy. After expulsion of their gas, the groups remain mass segregated, consistent with observations.