© The Authors 2024

1 Introduction

There is a consensus that the feedback from stars is key to explaining the physical properties of the observed galaxy population (e.g. Naab & Ostriker 2017). In fact, as seen in numerical simulations, the injection of energy and metals from stars into the surrounding medium regulates the efficiency of star formation (SF; e.g. Hopkins et al. 2011, 2014; Kimm et al. 2017; Marinacci et al. 2019) and prevents the assembly of overly dense stellar bulges (e.g. Agertz & Kravtsov 2015). Feedback also enriches the circumgalactic medium with products of stellar nucleosynthesis (e.g. Muratov et al. 2017; Nelson et al. 2019).

Early models of stellar feedback only accounted for the thermal energy from core-collapse supernovae (CCSNe; e.g. Dekel & Silk 1986; Katz 1992). However, other feedback processes that occur before the SN explosions (e.g. stellar winds and photoionisation) should also make important contributions (e.g. Dib et al. 2011; Hopkins et al. 2011; Brook et al. 2012; Stinson et al. 2012). Evidence for this pre-SN feedback comes from observations that gas is dispersed from the birth places of stars before the first SN explosions take place (e.g. Mathews 1966; Castor et al. 1975; Barnes et al. 2020; Chevance et al. 2022). Theoretical models suggest that preprocessing the immediate stellar environment could also strengthen SN feedback (e.g. Agertz et al. 2013; Geen et al. 2015b; Fierlinger et al. 2016; Hopkins et al. 2018; Kannan et al. 2020; Lucas et al. 2020).

What makes modelling stellar feedback very challenging is the wide range of spatial and temporal scales that regulate the physical processes (see e.g. Geen et al. 2023a). For this reason, different approaches are used to investigate feedback effects on different scales. On subgalactic scales, pre-SN feedback alone can already have a strong effect on the surrounding medium (e.g. Rogers & Pittard 2013; Fierlinger et al. 2016; Gatto et al. 2017), with stellar winds creating cavities of tens of parsecs (pc) in size inside H II regions (Dale et al. 2014). The strong metallicity dependence of the stellar wind strength (Vink 2022) implies that early feedback may be significantly less effective at very low metallicity, but this has not yet been thoroughly explored theoretically or confirmed observationally.

It has been observed that a single stellar system can have a strong influence on its surroundings. For instance, Ramachandran et al. (2019) showed that in the ‘wing’ of the Small Magellanic Cloud (SMC), a single Wolf-Rayet (WR) binary system is dominating the pre-SN feedback over the combined feedback of all OB stars. This highlights the importance of post-main sequence phases when determining stellar feedback.

In addition to improvements in the considered feedback processes, the subgrid injection model has also been improved since its initial implementation, resulting in the existence of different schemes. In some models, the energy is injected completely in the thermal component (e.g. Katz 1992; Stinson et al. 2006; Dalla Vecchia & Schaye 2012; Schaye et al. 2015), while other models inject kinetic energy and momentum (e.g. Dalla Vecchia & Schaye 2008; Dubois & Teyssier 2008; Gatto et al. 2015; Kimm et al. 2015), a mixture of both components (e.g. Davé et al. 2017; Hopkins et al. 2018; Fichtner et al. 2022; Chaikin et al. 2023), or switching between different methods (e.g. Agertz et al. 2020). The choice of the subgrid feedback scheme is often driven by numerical considerations. The strength of stellar feedback is usually either derived from libraries, for example STARBURST99 (Leitherer et al. 1999; Vázquez & Leitherer 2005; Leitherer et al. 2010, 2014) or BPASS (Binary Population and Spectral Synthesis, Eldridge et al. 2017; Stanway & Eldridge 2018), or calibrated by comparing the properties of simulated and observed galaxies.

The direct use of results from population synthesis codes neglects the fact that the feedback is altered by the unresolved surrounding medium with respect to the stellar ejection. Some of the feedback energy is likely dissipated (e.g. Garcia-Segura et al. 1996a,b) at subresolution scales. Cosmological full-box simulations can achieve spatial resolution of the order of 100 pc or higher (e.g. Feldmann et al. 2023), while simulations that zoom onto a selected region can achieve higher resolutions of a few pc or even subpc resolutions (e.g. Kimm et al. 2016; Kim et al. 2019; Agertz et al. 2020; Gutcke et al. 2022; Calura et al. 2022). However, energy dissipation at pc and subpc scale cannot be modelled correctly in these simulations when they are not resolved. In addition, energy dissipation can be due to the variable stellar wind during the lifetime of a star (Vink et al. 2001; Vink 2022), for instance when the fast WR wind catches up to the slow and dense wind of the red supergiant or luminous blue variable phase. However, the strength of the energy dissipation is unclear. Rey-Raposo et al. (2017) and Lancaster et al. (2021) argue that most of the feedback energy is radiated away, while Rosen et al. (2021) find that radiative losses are insignificant.

The interstellar medium (ISM) has a complex, multiphase structure (e.g. Cox & Smith 1974; McKee & Ostriker 1977; Mac Low & McCray 1988; Rosen et al. 1993; Silich et al. 2008), which leads to deviations from spherical symmetry in the evolution of stellar-feedback bubbles, as each bubble can expand further into regions of lower density (e.g. Zamora-Avilés et al. 2019). Finger-like structures (e.g. Garcia-Segura & Franco 1996; Dale et al. 2005) and champagne flows (e.g. Tenorio-Tagle 1979; Franco et al. 1990; Geen et al. 2023b) can also be triggered by radiation-hydrodynamical instabilities. These processes can only be modelled in multidimensional simulations. However, 3D simulations with multiphase structures are numerically expensive, which limits these studies to one or a few models.

In this work, we investigate the energy deposition from the feedback of massive stars at various distances from the stellar sources through high-resolution 1D radiation-hydrodynamics simulations. These simulations allow us to better estimate the radiative losses occurring at small scales and to begin to bridge the gap between stellar and galactic length scales. To this end, we generated a grid of over 2500 simulations of stellar feedback from single and binary stars in a homogeneous ISM. We include the feedback from stellar winds, photoionisation, and CCSNe. Simulations were performed with the grid-based fluid dynamics code PION (Mackey et al. 2021). As input for the PION simulations, we used the stellar models of massive stars described in Fichtner et al. (2022). As well as single stars, the grid includes models of binary systems, as most massive stars are born as binaries (Sana et al. 2012, 2013). We employed models with different initial surface rotation velocity of single stars and binary systems with different orbital periods. We investigated five different subsolar metallicities, with the aim being to replicate the conditions found in high-redshift galaxies. We used these simulations to investigate how much energy reaches different radii at different ambient densities. This is the energy budget that is available to drive feedback in galaxy simulations that do not resolve pc scales. Additionally, we analysed the ratio of kinetic to total energy, which is often a free parameter in subgrid feedback models.

The paper is organised as follows: In Sec. 2, we introduce the grid of stellar models. In Sec. 3, the numerical methods and initial conditions of the PION simulations are described. We present our analysis in Sect. 4, and discuss the implications of our findings in Sec. 5. We then review our main conclusions in Sec. 6.

2 Stellar models

2.1 Stellar grid

We adopt the stellar models from Fichtner et al. (2022), whose main properties we summarise here. The stellar models are computed with the stellar-evolution code Modules for Experiments in Stellar Astrophysics (MESA, v.8845, Paxton et al. 2011, 2013, 2015, 2018, 2019). The grid contains models of single, rotating massive stars as well as models of stars in binary systems. The inclusion of binary systems is of paramount importance as the majority of massive stars are born in binaries (Sana et al. 2012, 2013). Binary stars may experience phases of mass transfer when they fill their Roche lobes. The binary evolution creates evolutionary pathways and spectroscopic appearances which are not achievable for single stars (e.g. De Marco & Izzard 2017). We use stellar models in the mass range log₁₀(M_ini/M_⊙) = 0.9–2.2 in steps of log₁₀(M_ini/M_⊙) = 0.1 for single stars and in the range log₁₀(M_ini/M_⊙) = 1.0–2.2 in steps of log₁₀(M_ini/M_⊙) = 0.2 for the primaries of the binary systems. For the single stars, the grid spans a wide range of initial surface rotation velocities, namely 0 to 600 km s⁻¹ in steps of 100 km s⁻¹. For the binary systems, the initial rotation velocity is set to give tidal locking. Four orbital periods (P = 10^0.6, 10¹, 10², and 10³ days) are considered while, for simplicity, a fixed mass ratio of 0.6 is assumed (see Table 1 for the naming convention used in this paper).

In this work, we investigate the subsolar metallicities Z = 0.0001, 0.0004, 0.0008, 0.001 and 0.004, which are representative of the metallicities observed in high-redshift and/or low-mass galaxies (e.g. Maiolino et al. 2008; Madau & Dickinson 2014). Throughout this work, we state metallicities in absolute units and not relative to the solar metallicity. These metallicities correspond to the lower metallicities used in the binary grid, adding the single stellar models for the metallicity 0.0008 with respect to Fichtner et al. (2022). As in the previous paper, we stop the evolution of the stellar models at core He-exhaustion. In the remainder of the section we highlight the difference with respect to Fichtner et al. (2022).

In contrast to other stellar input properties, the terminal wind velocity υ_∞ is not provided directly by MESA. We calculate υ_∞ as in equation (2) of Fichtner et al. (2022), which was based on Leitherer et al. (1992) and Puls et al. (2008), but we modify the multiplicative prefactor between υ_∞ and the escape velocity. It is changed to follow the prescription of Eldridge et al. (2006), except that a factor of 2.0 is used for WR stars (identified by a temperature above 3 × 10⁴ K and a hydrogen abundance below 0.6).

For binary systems the individual stellar wind and radiation output must be combined into averaged, spherically symmetry quantities that can be input into 1D simulations, achieved by weighting the energy of each wind according to solid angle (see below). We add the mass- and metals-loss rates of the two stellar models and assume the temperature of the hotter star for the combined object. The luminosities are added with a correction factor so that the combined ionising photon luminosity is the sum of the two sources. To calculate the terminal wind velocity, we add the kinetic energy of the wind of the two companions, where we take into account energy dissipation due to wind-wind interactions in the binary system with an improved method (described in the following) compared to Fichtner et al. (2022). We assume that, at large distances from the system, only a fraction of the solid angle allows for the freely expanding wind of each component of the binary system. Based on the models by Eichler & Usov (1993) and Girard & Willson (1987, but see also Pittard & Dawson 2018), the region filled by the weaker stellar wind can be approximated by a cone. Half of the opening angle of that cone (θ) can be estimated numerically from the stellar mass-loss rate $\dot{M}$ and terminal wind velocity υ_∞ via: $$\theta =2.1(1-\frac{{\nu }^{2/5}}{4}){\nu }^{1/3},$$(1)

with $$\nu =\frac{{\dot{M}}_2{\upsilon }_{\mathrm{\infty },2}}{{\dot{M}}_1{\upsilon }_{\mathrm{\infty },1}},$$(2)

the indices 1 and 2 indicate the component with the stronger (given by $\dot{M}{\upsilon }_{\mathrm{\infty }}$) and weaker wind respectively, so that ν ≤ 1. The normalised solid angle of the cone of the weaker model 2 is then $$f_{\mathrm{w}\mathrm{w},2}=0.5(1-\cos \theta ),$$(3)

and the solid angle of the stronger wind model 1 is, consequently: $$f_{\mathrm{w}\mathrm{w},1}=1-f_{\mathrm{w}\mathrm{w},2}.$$(4)

By doing so, we find that the stellar model with the wind of higher momentum covers a larger solid angle, mimicking the increased momentum of its freely expanding outflow at large scales while accounting for energy dissipation in the region of colliding winds.

2.2 Supernova

In order to account for the variety of CCSN types, we employ the compactness parameter to infer which stellar models could explode as a CCSN. This parameter is a measure of the density (and binding energy) of the stellar core at bounce, evaluated at mass coordinate 2.5 M_⊙ defined as 2.5/R_2.5, where R_2.5 is the radius at bounce that encloses 2.5 M_⊙ in units of 10³ km. Models with a value of the compactness parameter below 0.45 can likely explode, while models ≥ 0.45 undergo direct black hole formation (O’Connor & Ott 2011) (although this simple prescription has been challenged recently; see Wang et al. 2022). To evaluate the compactness parameter, the stellar models are evolved beyond core He-exhaustion, except for the most massive models which are not expected to have a CCSN.

For the SN energies, we use a linear fit to the numerical results in Schneider et al. (2021), where the explosion energies appear correlated to the compactness parameter. For a simplified estimate of the ejected mass, we assume that the innermost 2 M_⊙ of the stellar models are locked in the remnant, while the remaining mass is ejected. The ejected mass of hydrogen, helium and metals is calculated from the respective masses of these components outside the innermost 2 M_⊙.

The compactness parameter is highly non-monotonic with stellar mass, resulting in the distribution of explosion energies shown in Fig. 1. All single stellar models with masses around 10 M_⊙ explode, while some of the models around 20 M_⊙ do not. Which models explode varies with metallicity and initial rotation velocity and there is no simple relationship between explosion (or not) and the initial stellar parameters. Consistent with part of the literature, for an initial mass of 25 M_⊙ all models explode. For binary systems, the lowest initial primary mass of a system where a CCSN takes place (either primary and secondary) is around 16 M_⊙. In systems with an initial primary mass of 16 M_⊙ and 25 M_⊙, all primaries and secondaries undergo a CCSN. For higher masses, CCSNe become rare for the primaries; however, we find (in broad agreement with Schneider et al. 2021) that, in binary systems, higher mass stars also explode (up to 60 M_⊙) compared with single stars. Nearly all secondaries explode where the initial primary mass is 40 M_⊙; above that mass CCSNe only happen for a few models.

The SN energies span an order of magnitude, ranging from 2.5 × 10⁵⁰ erg to 2.5 × 10⁵¹ erg, where more massive models overall tend to have higher energies. All models with a mass of 8 M_⊙ explode with nearly the same energy of around 2.5 × 10⁵⁰ erg. No clear trends with initial rotation velocity for single stars or orbital period for binary systems are visible.

By employing the compactness parameter to infer SN energies, the SN energy ejected by a stellar population varies with metallicity (because the compactness at core collapse depends in a non-linear way on many initial stellar properties, including metallicity). As an example, a population of single non-rotating stars ejects at Z = 0.004 an energy of 4.1 × 10⁴⁹ erg ${\mathrm{M}}_{\odot }^{-1}$ from CCSNe, while the same population ejects at Z = 0.001 only 2.4 × 10⁴⁹ erg ${\mathrm{M}}_{\odot }^{-1}$ (weighted by a Salpeter initial mass function (IMF) normalised to the total mass in stars above 8 M_⊙; see Sect. 4.2 for a description of the adopted IMF weighting). If instead one would employ the simple assumption that just all the single models (up to the highest mass of 160 M_⊙) undergo a CCSN with the canonical energy of 10⁵¹ erg, the population would provide 5.0 × 10⁴⁹ erg ${\mathrm{M}}_{\odot }^{-1}$, while for the whole grid (i.e. assuming every considered model undergoes a CCSN with 10⁵¹ erg and a binary mass fraction of 70 per cent), the population would provide 4.7 × 10⁴⁹ erg ${\mathrm{M}}_{\odot }^{-1}$, independent of metallicity.

Ejected masses increase with initial stellar mass and range from a few M_⊙ to around 25 M_⊙ and 40 M_⊙ for single and binary models, respectively. The few cases of extreme mass ejection of binary SN might originate from adopting the compactness parameter as the explosion criterion. Overall, single models with a lower rotation velocity have higher SN ejecta masses due to the lower mass-loss rates over the course of their lifetime. In closer binary systems, less mass is ejected by the primary during the CCSN and more by the secondary compared to wider systems, due to the stronger mass transfers. However, the extent of these differences varies with metallicity. The metal fractions of the ejected masses also increase with initial mass. For single stars, the ejecta of fast rotating stars has higher metal fractions, up to where nearly all the ejected mass is found in metals for the most massive, H-depleted stars at the lowest metallicity. For binary systems, the CCSNe of primary stellar models have on average a higher metal fraction in their ejecta than the CCSNe of their secondaries. In binary systems where both the primary and secondary undergo a CCSN, the metal fraction of the primary is higher than that of the secondary, due to the higher mass-loss rates by both Case A mass-transfer¹ as well as stellar wind. Low mass secondaries have the lowest metal fractions. In the most extreme cases of stellar models undergoing SN explosions, up to nearly 18 M_⊙ (single) and 30 M_⊙ (binary) of metals would be ejected and then enrich the surrounding ISM. The progenitors of these cases are stripped of their outer layers already during their lifetime.

Overall, SN energy, mass and metal ejecta vary with metallicity. However, no strong trends with metallicity are visible (except a decrease of the mass ejecta with increasing metallicity for single stars). Binary systems give rise to CCSNe from progenitors with higher initial mass than single stars, leading to more energetic events. Moreover, binary star systems can potentially host a sequence of two successive SN events.

Fig. 1 Explosion energy of each stellar model if it has a CCSN, where the colour of the squares indicates the SN energy. White squares indicate models that do not have a CCSN. The left column shows the single models, while the two right columns show the primary and secondary models of binary systems, respectively. Each row shows one metallicity. The x-axis shows the different names according to Table 1. The initial stellar masses are shown on the y-axis; we note that the binary grid has a coarser sampling in initial mass than the single grid. The highest initial masses are not shown as they have no CCSN.

3 Numerical methods

3.1 One-dimensional simulations

We perform 1D, spherically symmetric simulations with the radiation-hydrodynamics code PION, using the source term for evolving stellar winds described in Mackey et al. (2021). The stellar models described in Sect. 2 provide the stellar mass, luminosity, effective temperature, mass-loss rate, rotation velocity, radius and element abundances as a function of time. These quantities are used to calculate (i) the ionising photon luminosity and mean energy per ionisation for the photoionisation calculation, and (ii) wind density, pressure and terminal velocity imposed in a boundary region around the origin. The stars impact the surrounding ISM during their lifetime via kinetic winds injected into the innermost two cells, and via photoionisation energy injected throughout the H II region, where most is injected at the Strömgren radius due to the high opacity. The wind terminal velocity is calculated as described in Sect. 2.1. The ionising photon luminosity, Q₀, is calculated by matching the temperature-dependent Q₀ data tabulated in Diaz-Miller et al. (1998) using a power-law fit below 33 000 K, scaled to the appropriate stellar radius. At higher temperatures a blackbody provides a good estimate of Q₀. The mean energy per ionisation is calculated assuming a blackbody spectral shape above 13.6 eV. The non-equilibrium photoionisation scheme was tested in Mackey (2012) and benchmarked against analytic solutions and other codes in Bisbas et al. (2015).

The radiative heating and cooling scheme was described in Mackey et al. (2013) and consists of: radiative heating by photoionisation, metal-line cooling assuming collisional ionisation equilibrium (CIE) with scaled solar abundances, Bremstrahlung in ionised gas, fine-structure line cooling of C⁺, O in warm neutral gas, together with line cooling from CNO ions in photoionised gas, similar to the model of Henney et al. (2009).

For this project we modified the 2D/3D nested grid algorithms (Mackey et al. 2021) to also work for 1D spherical coordinates, focused on the origin. The 1D nested grid is a simple reduction of the multidimensional algorithm with the addition of appropriate geometric source terms. We follow the description in Mackey et al. (2021) of prolongation and restriction (Tóth & Roe 2002) at the coarse/fine level boundaries, and ensure consistency between levels using the flux correction algorithm of Berger & Colella (1989).

The coarsest level 0 covers the full domain of size ℒ with radial coordinate r ∈ [0, ℒ]; level 1 covers r ∈ [0, ℒ/2], level j covers r ∈ [0, ℒ/2^j], down to the finest level. The coordinate singularity at r = 0 is treated with reflecting boundary conditions (the surface areas of this boundary is anyway zero so the flux through r = 0 is always zero), and in the pre-SN phase the stellar-wind boundary condition is imposed on the first two cells adjacent to r = 0. The outer boundary of each level j is obtained by interpolating data from the appropriate cells of level j − 1 and these are updated every timestep. After each step of level j, averaged states are sent to level j − 1 to overwrite cell data in the subdomain r ∈ [0, ℒ/2^j] of level j − 1. An integration scheme is used that is accurate to second order in time and space (Falle et al. 1998), using a TVD linear interpolation and minmod slope limiter.

3.2 Initial conditions and simulation parameters

The simulation suite is performed in a computational domain of ≈ [0, 66] pc, which ensures that the feedback bubble does not leave the domain during the star’s lifetime. We use six levels of static mesh-refinement with 256 grid cells per level and a factor of 2 refinement with each level as described above, resulting in a cell size of 0.0081 pc near the origin to 0.26 pc on the coarsest level at r > ℒ/2. In cases where the free-streaming region of the wind at 10 per cent of the (primary) stellar lifetime is not resolved with at least 5 grid cells (happening only for the lowest-mass stellar models), we add additional refinement layers in the centre until the criterion is met.

We use a uniform ISM with four different densities: 2.2, 4.4, 11, and 22 × 10⁻²² g cm⁻³ (corresponding to a hydrogen number density of around 100, 200, 500, and 1000 cm⁻³). These densities are selected to approximate the range of SF densities in zoom-in simulations of galaxies. While in general the ISM is actually turbulent and structured, this cannot be captured in spherically symmetric calculations and a uniform ISM is as realistic as any other feasible configuration. Multidimensional simulations would be preferable, but at the cost of orders of magnitude more computation and storage. Even if these are possible, 1D simulations provide a valuable point of reference to investigate the importance of multidimensional effects such as dynamical instability and turbulent mixing.

The temperature of the ambient medium is initialised at 50 K, which is within a factor of two of the equilibrium temperature for the different setups. A temperature floor of 10 K is set. Initial ambient metallicities are set to the initial surface values of the stellar models. Tracers for H, He, C, N, O, Z, and H+ are passively advected, where in Z the remaining metals are stored. Non-equilibrium ionisation of H⁺ is calculated including pho-toionisation, collisional ionisation and radiative recombination, using a backward-difference implicit method. Ionising photon radiation is calculated by raytracing outwards from the star at the origin to track the attenuation. Outputs are generated every 10¹¹ s and the simulations are run until the end of the stellar lifetime. In total, the PION simulation grid consists of 1960 runs of single stars and 560 runs of binary systems.

To simulate the subsequent SN, we use the PION output at the end of the stellar lifetime. We overwrite the innermost 10 cells, using a 2-component profile for the density and velocity following Whalen et al. (2008), where the two parameters of the model are found iteratively. For the SN of the primary in a binary system, the SN injection is shifted by two cells to avoid conflicts with the wind injection of the secondary. The hydrogen and helium fractions are injected into their corresponding tracer fields, while the metals of the SN are traced with a separate field. Outputs are generated in multiples of the time step, where the early evolution is sampled with a higher frequency. The SN simulations are evolved until 5 Myr, or longer until the shell of material swept up by the SN reaches 20 pc. For the binary grid, the simulation is restarted from the last output of the primary lifetime to simulate the remaining secondary lifetime, adding, if required, a primary SN. In case the secondary has a CCSN, its evolution is simulated as in the single star case. In cases with a primary SN and no secondary SN, we ensure that the simulation is run for at least 5 Myr after the primary SN.

4 Results

4.1 Evolution of selected models

First, we look at selected models in more detail: These models are single stellar models with Z = 0.004 and without stellar rotation, in a uniform medium of 100 cm⁻³.

4.1.1 Evolution of a 32 M_⊙ star

We show the density, temperature and velocity profiles of one PION simulation to examine the general evolution in the circumstellar medium around a massive star. In the top panel of Fig. 2, five regions can be seen, which are indicated by the encircled numbers above the top panel. These are from left to right: 1: free-expansion of the stellar wind, 2: hot shocked wind material, 3: H II region, 4: cold swept-up medium, and 5: unperturbed surrounding medium. In the middle panel, which shows the last output during the star’s life, shells of wind ejecta can be seen as locally overdense and cooler regions (in Fig. 2 at around 10 pc). At the end of the stellar lifetime, the feedback bubble is around 28 pc. Therefore, the pre-SN feedback displaced already around 3 × 10⁵ M_⊙ of circumstellar material. The bottom panel shows the profile at the end of the SN simulation which was evolved for 5 Myr, containing a hot, low-density medium surrounded by a shell of cold swept-up gas. The overdense shells and other structures found at the end of the star’s life are swept by the SN ejecta. In the 5 Myr after the SNe, the bubble expands to around 40 pc, displacing around 9 × 10⁵ M0. In Appendix A, we show profiles for the evolution of a binary system with two SNe for comparison.

We shall now discuss the energy flow through surfaces of spheres at different radii, and equate the energy that is available for feedback in cell sizes above a few pc. These radii are denoted r_s and we choose r_s ∈ [5, 10, 20] pc for this work because they cover the full range of behaviours (i.e. barely any dissipation to strong dissipation) at the gas densities considered. The kinetic and thermal energy densities, u_kin and u_th, at a radius of r_s and time t since a star’s birth are $$u_{\text{kin }}(r_{\mathrm{s}},t)=0.5\rho (r_{\mathrm{s}},t)\upsilon {(r_{\mathrm{s}},t)}^2,$$(5) $$u_{\mathrm{t}\mathrm{h}}(r_{\mathrm{s}},t)=\frac{p(r_{\mathrm{s}},t)}{(\gamma -1)},$$(6)

with ρ the (mass) density, υ the 1D-velocity, p the pressure, and $\gamma =\frac{5}{3}$. The energy flux through the spherical surface, e_i(r_s, t) is then $$e_{\mathrm{i}}(r_{\mathrm{s}},t)=u_{\mathrm{i}}(r_{\mathrm{s}},t)\upsilon (r_{\mathrm{s}},t)4\pi r_{\mathrm{s}}^2,$$(7)

and the cumulative energy flow through this sphere is the discretised time integral of the flux: $$E_{\mathrm{i},\mathrm{c}}(r_{\mathrm{s}},t)={\displaystyle \sum _{t_{\mathrm{j}}\le t}}{e}_{\mathrm{i}}(r_{\mathrm{s}},t_{\mathrm{j}})\mathrm{\Delta }t,$$(8)

with i indicating the kinetic or thermal energy density and Δt the time interval with respect to the last output. We only consider radii of r_s = 20 pc and below, as we find (see e.g. Sect. 4.2) that nearly all energy is dissipated already for 20 pc in the cases of higher densities. A caveat in the flow calculations is introduced by the photoionisation energy. It is injected mostly at the Strömgren radius; if that radius is larger than r_s, the photoionisation energy is not accounted for.

Figure 3 shows the cumulative energy flow (total, and split into the kinetic and thermal component) through r_s = 5, 10 and 20 pc of this model. In the pre-SN flow, before the end of the stellar lifetime at around 6.1 Myr, we can see that the kinetic energy only flows through each sphere at one instance in time when the cold swept-up mass crosses that radius, while the thermal energy flow is more continuous. However, also the thermal energy shows a noticeable rise at two times. The first occurs when the swept-up mass crosses r_s, (where also the kinetic energy increases strongly) and the second rise occurs when the first shocked wind material reaches r_s. Higher radii are reached by the swept-up mass at a later time, for example, it takes around 1 Myr until the first feedback energy reaches 10 pc. The kinetic energy is higher at higher radii at the end of the stellar lifetime (i.e. before around 6.1 Myr), while the total energy shows a more complicated behaviour with radius due to the interplay of energy dissipation and the injection of photoionisation energy at radii larger than r_s. The fraction of kinetic energy to the total energy at the end of the stellar life increases from around 4 per cent at 5 pc to around 50 per cent at 20 pc. The shocked wind material never reaches 20 pc, so we see there mainly the contribution of the photoionisation energy to the total energy.

The subsequent SN explosion of the stellar model at around 6.1 Myr leads to a sharp rise in the cumulative energy flow, dominating the total energy budget at all radii. Due to the low density inside the feedback bubble at the end of the stellar lifetime, the energy from the SN explosion can rapidly reach 20 pc without conspicuous radiative energy losses. At 5 pc, nearly all the SN energy is kinetic, while at 20 pc part of the energy has been converted to thermal energy through shock heating.

Although SN energy is the predominant factor in the overall energy budget, the pre-SN feedback occurring during the star’s lifespan is crucial. We present a comparison involving the same SN, but situated in a uniform medium unaffected by the star’s pre-SN feedback, in Fig. 3. This is depicted with symbols positioned at the conclusion of the fiducial simulation. For all radii, the simulation without pre-SN feedback has a lower total energy. The differences increase with radius. The case without pre-SN feedback has dissipated the energy by one and two orders of magnitude for 10 and 20 pc, respectively. On the contrary, even at 20 pc, our simulation with pre-SN feedback has retained most of its energy. This results in a difference of two orders of magnitude in the energy at 20 pc. Further comparisons with the case without pre-SN feedback are shown in the Sect. 4.5.

Fig. 2 Density, temperature, and velocity profiles of a PION simulation with a single star model of Z = 0.004, υini = 0 km s−1, and log10(M/M⊙) = 1.5 in a uniform medium of n = 100 cm−3 at 5 Myr (top), at the end of the stellar lifetime (middle), and 5 Myr after the subsequent SN (bottom). Dashed lines indicate negative velocities. The circled numbers at the top indicate the five different regions found in the density profile due to pre-SN feedback.

Fig. 3 Cumulative energy flow through rs = 5, 10, and 20 pc for a single star model of Z = 0.004, υini = 0 km s−1, and log10(M/M⊙) = 1.5 in a uniform medium of n = 100 cm−3. The solid line shows the total energy, while the dashed and dotted lines show the kinetic and thermal energy component, respectively. A SN explosion occurs around 6.1 Myr. The star symbols (labelled ‘uniform’) show a comparison of the total energy with respect to a SN run in a uniform medium without pre-SN feedback beforehand.

Fig. 4 Density profiles of our simulations with single star models of Z = 0.004, υini = 0 km s−1, in a uniform medium of n = 100 cm−3 at the end of the star’s lifetime. For clarity, we show only the profiles of very second model of this sequence of stellar models. The labels indicate the logarithm of the initial masses in M⊙.

4.1.2 Density profiles of one sequence of stellar models

Figure 4 shows the density profiles of the half of the stellar models at Z = 0.004, υ_ini = 0 km s⁻¹, in a uniform medium of n = 100 cm⁻³ at the end of the stellar lifetime. It is remarkable that the size of the feedback bubble (i.e. the location of the peak of the swept-up mass at around 27 pc) is similar across nearly the whole mass range. To illustrate the origin of the similar size, we compare the size of the feedback bubble (measured by the outer edge of the swept-up shell) to the relation for the radial evolution of an H II region (Spitzer 1978) in Fig. 5. For estimating the Spitzer radius, we use the properties (luminosity, temperature, lifetime) of the simulated stars and the temperature of the H II region, but we neglect any time evolution of these. The H II region radius calculated from the Spitzer solution is also nearly independent from the stellar mass for the highest mass models. The small dependence that is visible is similar in both and the radii agree reasonable well considering the approximations that were made. The reason for the nearly mass independence is that the lower luminosities of ionising photons for smaller stellar masses is compensated by longer lifetimes (in that stellar mass range). We note that the H II regions of the lower mass models have reached pressure equilibrium with their surroundings, and so the neutral shell of swept-up material becomes geometically thick, invalidating the assumptions in the Spitzer solution (see also Williams et al. 2018). For higher mass stars this is not the case and a thin shell of swept-up gas bounds the H II region (see Fig. 4).

Fig. 5 Comparison of the radius of the outer edge of the shell of swept-up material to an estimate from the Spitzer radius. The radius is calculated from the single star models of Z = 0.004, υini = 0 km s−1, in a uniform medium of n = 100 cm−3 at the end of the star’s lifetime.

4.1.3 Evolution of a 100 M_⊙ star

The energy flows for our single-star models with a metallicity of Z = 0.004 and no stellar rotation are depicted in Fig. 6. The model with 100 M_⊙ shows that the swept-up mass reaches the different r_s earlier compared to the 32 M_⊙ model. The progression of energy flow in both models is similar, yet with higher values for the 100 M_⊙, up to about 3 million years, which marks the onset of the WR phase. Following this, there is a pronounced increase in cumulative energy, with shorter intervals compared to the initial gradient. At a radius of r_s = 5 pc, the energy composition before the WR phase is predominantly thermal as the termination shock lies within this radius. After the WR phase, the kinetic energy becomes dominant as the termination shock moves beyond this radius. Conversely, at r_s = 20 pc, the trend is opposite: kinetic energy is prominent before the WR phase, but post-phase, thermal energy takes precedence. This indicates that the kinetic energy from the WR wind is transformed into thermal energy as it travels outward. The shocked wind material extends to 20 pc, with the WR wind pushing a shell of early wind material beyond even 20 pc. Notably, in this model, the total energy diminishes as the radius increases, and it does not undergo a SN explosion.

4.1.4 Evolution of a 10 M_⊙ star

As opposed to the 32 M_⊙ model, Fig. 6 shows that the cumulative thermal energy flow driven by a 10 M_⊙ star dominates over the kinetic energy flow at all radii during most of the star’s lifetime. Only when the swept-up mass reaches r_s is the kinetic energy important. This time is later than in the 32 M_⊙ case. The total energy is increasing with radius in the pre-SN phase and exceeds the injection of the stellar wind for 10 and 20 pc as the energy budget is completely dominated by photoionisation. The photoionisation energy couples mostly at the Strömgren radius, and therefore the energy is not yet injected at the smaller radii. The subsequent SN dominates the total energy flow by more than one order of magnitude and is decreasing with distance. As the stellar model created a smaller wind cavity compared to the 32 M_⊙ model, the SN energy remains mostly in the kinetic component only out to ~5 pc. At 10 pc, most of SN energy has been converted to thermal energy. Further out at 20 pc, the total energy decreases by over one order of magnitude due to thermal dissipation, only a small kinetic energy flow remains.

Fig. 6 Cumulative kinetic (top) and thermal (bottom) feedback energy for different stellar initial mass as a function of time for the single star models with Z = 0.004 and no initial rotation velocity is shown with coloured bars for the whole time of the PION simulation. The different columns show different initial ambient densities, while the lower three rows vary rs. The star symbol indicates the time of the SN, while black regions visualise outflowing energies below 1046 erg or inflows. The upper row shows IMF-weighted results (normalised to the total mass in stars above 8 M⊙). Additionally, the sum of the energy inputs from SN, stellar winds and an estimate of the photoionisation energy is shown (labelled ‘Star’), and the energy input excluding the photionisation energy (labelled ‘no ph’).

Fig. 7 IMF-weighted (normalised to the total mass in stars above 8 M⊙) cumulative energy flow for the single star models with metallicity Z = 0.004 and no initial rotation with a background density of 500 cm−3. Additionally, the sum of the energy inputs from SN, stellar winds and an estimate of the photoionisation energy is shown (labelled ‘Star’), and the energy input excluding the photionisation energy (labelled ‘no ph’).

4.2 Overview of the cumulative energy flow

In addition to the flows of the different stellar models, Fig. 6 shows IMF-weighted values using a Salpeter IMF for our models between 10^0.9 M_⊙ for single stars, or 10 M_⊙ for the primary star mass in binary systems, to 10^2.2 M_⊙. We note that we normalise for a population of only massive stars with masses between 8 M_⊙ and 10^2.2 M_⊙ to avoid the more uncertain low mass end of the IMF. For binary systems, we apply the IMF to the primary star. When including the SN runs in the IMF, we neglect the evolution after the start of inflows characterised by a decrease in the energy.

We compare the IMF-weighted cumulative flows with the time-integrated input energy from the stellar wind, SN and photoionisation energy. The photoionisation energy is estimated by integrating the Planck spectrum from 13.6 eV to 100 eV, while subtracting 13.6 eV from each photon, which is the maximum power available to drive feedback. This assumes that most ionising photons ionise a H atom and heat the gas with the remaining photon energy, and neglects radiative losses.

4.2.1 Radial distance evolution of the cumulative energy flow

First we look at single star models with Z = 0.004 and no initial rotation velocity. We start with the total IMF-weighted cumulative flows of the density 500 cm⁻³, visible in Fig. 7. The energy reaching r_s decreases with distance, resulting in orders of magnitude lower energy at 20 pc. In addition, the energy crosses more distant values of r_s later. This means that the first time the energy arrives at a large distance is a few million years after the birth of the stellar population. As for the separate models, the SN energy also dominates the IMF-weighted energy budget over the energy crossing r_s during the stellar lifetime.

Compared to the total energy input of winds, photoionisation, and SN, the energy flows from 1D simulations are significantly lower (see the red line in Fig. 7). However, the energy is similar to the injection of only wind and SN energy in both total energy and its time evolution, when evaluated at r_s = 5 pc. This highlights the low coupling efficiency of the photoioninzation energy. For the smallest radius energy losses are negligible (especially at lower densities). At early times (below around 2 Myr), the energy flows are also lower with respect to only wind and SN injection. This highlights that main sequence stellar winds are not efficient, while later phases and SNe are.

4.2.2 Overview of the cumulative energy flow at fixed metallicity

We compare the different background densities and split the energy into the kinetic and thermal components (see Fig. 6). Before the onset of the first SN, corresponding to the early times in Fig. 6, the late evolutionary phases of the most massive stars are the biggest contributors to the kinetic energy. This can be seen from the sharp increase in the IMF-weighted kinetic energy at the time when they enter the post-main-sequence phase. Overall, the cumulative IMF-weighted kinetic energy flow before SNe can be well approximated by a single injection due to the steep increase caused by the post-main sequence phase. Similar to bigger values of r_s, higher background densities reduce the kinetic energy reaching r_s as both higher densities and larger distances increase the mass within r_s. The kinetic energy (including SNe) at 5 pc is nearly independent of density, as the pre-SN feedback created a lower density feedback bubble beforehand, and therefore the central densities are independent of the initial densities. Where the IMF-weighted kinetic energy decreases strongly with increasing measurement radius, r_s, it indicates that the SN becomes radiative in between these radii. As expected, this radius decreases with increasing density.

Similar to the IMF-weighted kinetic energy, the IMF-weighted thermal energy before the first SNe is dominated by the feedback from the most massive models. For high background densities or large radii, the late phases of the most massive models become more important for the total thermal energy because part of the wind energy is converted from kinetic to thermal through shock heating. At low densities and large radii (see lower left and IMF-weighted subpanels in Fig. 6), the total energy is dominated by thermal energy from the late stages of the massive stars. The thermal energy flowing through 5 pc does not vary strongly with density.

The IMF-weighted thermal to kinetic energy ratio changes with measurement radius. The intermediate mass stars (25 M_⊙ to 40 M_⊙) are the strongest contributors to the total energy as they are the ‘best of both worlds’. They have SNe like stars around 10 M_⊙, but the SN energy is higher and their winds are stronger. While the most massive stars around 100 M_⊙ have even stronger winds, they do not have SNe. This, in combination with the IMF, makes them important for the total energy feedback from a stellar population.

Figure 8 is the same as Fig. 6 except that we show results for the grid of binary systems with the shortest period. Trends with density and radius are similar for the binary systems as for the single models. Also the comparison of the input energy with the measured energy flow is similar. Some results, however, are specific for the binary systems. In the most massive binary systems, mainly the primary star contributes to the kinetic and thermal energy flow. Some models experience a phase of strong mass loss during mass-transfer events, which can lead to the formation of circumstellar shells. During the time when only the secondary star lives, there is barely any energy increase for most models (without SN). Compared to the single models, for binaries there is also the special case of runs with two subsequent SN when both primary and secondary undergo a CCSN. In these cases, the secondary SN evolves in a cavity formed by the primary SN, which reduces energy dissipation. This is especially important for the case of the highest density, where the SN of the secondary star of the 10^1.4 M_⊙ primary retains nearly 10⁵¹ erg even until 20 pc, while the first SN in the system lost most of its energy. However, this is not the case for the secondary SN of the system with a 10^1.2 M_⊙ primary due to the lower feedback.

Overall, when considering only the evolution during the stellar lifetime, we can say that the hot bubble is resolved with 5 pc. However, with 10 pc it is not very well resolved, and with 20 pc not at all. On large scales these results favour a kinetic energy injection method, because the hot bubble is not resolved. This is in agreement with recipes of feedback that do not inject thermal energy as it is being radiated away, but instead kinetic energy (e.g. Dalla Vecchia & Schaye 2008). The plots indicate that the main sequence of the stellar model is not a main energy contributor. The thermal energy is dominated by photoionisation, the kinetic energy by the post-main-sequence phase during the star’s lifetime. Far away from the star, a delay of the arrival of significant energy of more than 1 Myr occurs. The kinetic to total energy ratio changes with time, density, radius and stellar track.

Fig. 8 Same as Fig. 6, but for binary systems with Z = 0.004 and log10(P [days]) = 0.6. The star symbol indicates the time of a SN (primary or secondary), while the diamond symbol indicates the end of the primary lifetime without a SN.

Fig. 9 IMF-weighted (normalised to the total mass in stars above 8 M⊙) cumulative energy flow for the single star models with metallicity Z = 0.0001 and no initial rotation. From top to bottom, the rows in each panel show kinetic, thermal, and total energy flow, while the columns show different initial ambient densities. Additionally, the sum of the energy inputs from SN, stellar winds and an estimate of the photoionisation energy is shown (labelled ‘Star’), and the energy input excluding the photionisation energy (labelled ‘no ph’).

4.2.3 Metallicity comparison of the cumulative energy flow

At lowest metallicity (see Fig. 9), the SN dominates the kinetic energy flow by orders of magnitude over the pre-SN feedback (stellar winds and photoionisation). The thermal energy is similar for all densities for 5 and 10 pc, indicating that the SN stays adiabatic until 10 pc. The wind bubble is smaller than 5 pc for 1000 cm⁻³ as can be seen from the IMF-weighted thermal energy after SN. The decrease of the total energy between 10 and 20 pc for the highest density is due to energy losses as the SN becomes radiative.

Next, we compare the thermal and kinetic energy for the metallicities Z = 0.0001 with the already discussed case of Z = 0.004. Relative to Z = 0.004, the kinetic energy is similar in Z = 0.0001 for 5 pc and 100 cm⁻³. The reason is that the SN energies are similar for the metallicities and the wind bubbles all extend out to 5 pc making the initial conditions less important. At the combination of 5 pc and 1000 cm⁻³ and 10 pc and intermediate densities, the kinetic energy is lower at lower metallicity as SN energy is thermalised already at a smaller radius. In addition, the lack of a WR phase is reducing the kinetic energy input. The kinetic energy evolution with radius of the density 100 cm⁻³ has the peculiar behaviour that it is higher at r_s = 20 pc than at 10 pc, while the thermal energy for this density shows the opposite behaviour. One can see how the thermal to kinetic energy ratio changes with metallicity, for example for 500 cm⁻³ and 10 pc the kinetic energy fraction is lower at lower metallicity. The total energy for Z = 0.0001 shows barely any energy loss between 10 and 20 pc for the density 100 cm⁻³ as the region inbetween is completely dominated by the H II region at the end of the stellar lifetime. In contrast, the total energy for the ISM density 500 cm⁻³ decrease between 10 pc and 20 pc. In that density, the peak of the swept-up mass is just around 20 pc at the end of the stellar lifetime, and the lowest mass models even exhibit densities of the order of the ambient density then. The higher densities lead to stronger energy dissipation.

As the SN energy varies with the radius for the different densities, but the injection is the same for the combinations of radius and density, the pre-SN wind is important as it sets the cavity size. There is less radiative cooling at lower metallicity, so the H II region is hotter which leads to stronger feedback. Therefore there is still feedback at 20 pc for 500 cm⁻³, as also SNe expand to a bigger radius in a fixed time at lower metallicity. At 5 pc, the total energy for nearly all cases is the same, only the thermal to kinetic ratio is different at the largest densities. It is surprising that the total energies look so similar with metallicity at 5 pc; however, the 5 pc results only show how much radiative cooling one has. If one resolves 5 pc (so cell size is around 1 pc) one can maybe directly inject the stellar feedback (wind, photoionisation, SN) without using an intermediate step such as the 1D simulations presented here.

Fig. 10 Cumulative total energy flow for all models through rs = 10 pc, weighted by an IMF (normalised to the total mass in stars above 8 M⊙). The different subpanels show different combinations of metallicity and ambient density. The grey horizontal line separates single from binary models.

4.2.4 Comparing the cumulative energy flow of all models

A comparison of the IMF-weighted total energy flows for all single and binary models can be found in Fig. 10, with the example of the energy flux measured at 10 pc from the star. The sharp energy increases visible are due to SN, and the earlier increases in same models are due to post-main-sequence phase. Additionally, a time delay until the first energy flow is visible, especially for the lower densities. A striking result is the remarkable similarity over the whole grid of models in their total energy feedback. This especially means that the feedback at low metallicity is comparable to the one at high metallicity due to SNe. There is no strong dependence on metallicity. Differences can be seen at early time due to the different wind strength at the different metallicities. However, a trend with rotation velocity is visible, with fast-rotating models that tend to have higher feedback (as they lose more mass, as we already saw beforehand). Variations in the feedback for the binary system seem to be mainly driven by the SNe. Also in the overview, the feedback is decreasing, as expected, with the radius at which it is measured.

Regarding the fraction of cumulative kinetic energy to the cumulative total energy (see Fig. 11), lower radii have higher fractions as the kinetic energy gets thermalised outwards. For 5 and 10 pc lower densities have higher fractions as less mass can get shock heated. Additionally, at 10 pc a trend with the stellar rotation velocity can be seen (except for the highest metallicity). As there are more stellar models with a WR phase at higher rotation velocity, more kinetic energy is injected during the stellar lifetime. This leads to a bigger cavity, which enables more of the kinetic SN energy to be retained. At 20 pc, the overall energy dissipation removes trends, especially at the higher densities. At early times, the kinetic fraction tends to be lower, while at later times SN and post-main-sequence winds lead mainly to higher kinetic energies, increasing this fraction.

4.3 Feedback efficiency

The fraction of input energy (in the form of stellar winds, SN and photoionisation) flowing through r_s between 5 and 20 pc, the feedback efficiency, is mostly around a few per cent, and remarkably similar for all different tracks (see Fig. 12). This reflects the low coupling efficiency of photoionisation. The two highest densities at 20 pc have nearly all lower efficiencies, as radiative losses are strong. We note that the binary and fast rotating single models have an efficiency that is around one order of magnitude higher than the slower rotating single models at 20 pc and highest density. The metallicity Z = 0.004 without stellar rotation has a higher efficiency as it has more SN, while the one with 300 km s⁻¹ has the least number of SN. A caveat in the efficiency estimation is that the photoionisation energy can be deposited at a higher radius than 5 pc, which is therefore not captured in the kinetic and thermal energy flow measurements. However, it is included in the input energy estimation for all radii. We saw already that the photoionisation has such a low efficiency and the total flow at 5 pc is following closely the input energy of winds and SN. Therefore, it is worthwhile to look at the efficiency without considering the photoionisation energy in the input (see Fig. 13; however, the PION simulation flows all include it). Considering only stellar winds and SNe is also similar to typical subresolution models of feedback in galaxy formation simulations without radiation transfer. This efficiency therefore gives an indication of the different total energy budget available for driving feedback when injecting stellar winds and SN directly or using results from circumstellar medium simulations like these.

At 100 cm⁻³ and 5 pc, the feedback energy taken directly from MESA are similar to those obtained with PION so (nearly) all the energy is retained. This indicates that if one can resolve 5 pc and can include photoionisation, one can directly use the results from the stellar models. However, at the same density, but a radius of 20 pc, only around half of the original input energy is still available. At higher density, the energy ratio drops further. At 1000 cm⁻³ and 20 pc, especially the slower rotating single models contain orders of magnitude less energy, and binary models up to 30 per cent of the input energy. In these cases of orders of magnitude less energy, a subresolution model might just inject the remaining SN momentum. While it would be possible, this efficiency never reaches values above one as the energies are dominated by SNe. Additionally, the photoionisation (that is neglected in this plot) is mostly not injected at the lowest radii, while at larger radii energy dissipation is already important, limiting the achievable efficiency.

Fig. 11 Ratio of cumulative kinetic energy flow through three rs values to the cumulative total energy flow at the final time. The values are IMF-weighted. Each subpanel shows all combinations of tracks and metallicity, while the different subpanels show different combinations of ambient density and radius.

Fig. 12 Ratio of cumulative total energy flow through three rs values to the cumulative energy input through winds, SN, and photoionisation at the final time. The values are IMF-weighted. Each subpanel shows all combinations of tracks and metallicity, while the different subpanels show different combinations of ambient density and radius.

Fig. 13 Same as Fig. 12, but the energy input only considers winds and SN, but no photoionisation. We note the different colour scale.

4.4 Feedback energy of stellar populations

We can further combine the results into stellar populations. For this, we use the metallicity-dependent initial rotation velocity distribution for single stars and the orbital period distribution for binary stars of Fichtner et al. (2022). The first was based on observations of the rotation velocities of young stars in the Galaxy, Large Magellanic Cloud and SMC, while the orbital period distribution was based on the results of Kobulnicky et al. (2014) and Almeida et al. (2017). Additionally, we assume a binary mass fraction of 70 per cent. The remaining parameters are then metallicity, ambient density and r_s dependent.

We show the ranges of the total energy and the kinetic fraction in Table 2 for the three r_s values. Additionally, Fig. 14 shows the cumulative population energy as a function of metallicity for the different densities and r_s values. At r_s = 5 and 10 pc, all populations have remarkable similar total energies, and only at 20 pc a clear density dependence can be seen. For the highest densities at 20 pc, the energy budget is completely dominated by the binary stars due to their order of magnitude higher retained energy compared to slow rotating single stars. Opposite to what might be the expected behaviour, the energy is lowest for the highest metallicity at this radius and density. The reason is that cooling is stronger at higher metallicity. For the contribution of the single stars, there is the additional effect that lower metallicities contain higher fractions of faster rotating stars. The dip visible in Fig. 14 for the highest density and radius is originating from the binary models. They have a higher energy at Z = 0.001 than Z = 0.0008 due to differences in the number of SNe and feedback bubble size.

The typical kinetic energy fraction decreases with radius. At r_s = 5 pc the kinetic energy is for all cases dominating over the thermal energy, with the lower densities having the higher kinetic energy fractions. Additionally, a metallicity dependence due to cooling can be seen for r_s = 10 pc, where stronger cooling leads to higher kinetic energy fractions.

4.5 Comparison to a case without pre-SN feedback

In Sect. 4.1.1, we see a steep increase of the total energy at the time of the SNe, which gives rise to the question of the importance of the pre-SN feedback phase. To investigate this, we simulate the SNe of single stars with the highest metallicity and no rotation velocity in an unperturbed medium. Fig. 15 shows that less energy reaches the different radii without pre-SN feedback. For most cases the differences are more than one order of magnitude, as the low density feedback bubble helps to retain energy. For the lower two densities at 5 pc, the total energy without SN is 10 to 30 per cent of the one in the fiducial case. In the other two cases with a higher fraction (highest two densities at 20 pc) also the standard model has radiated away most of the energy already. Overall, it is visible that the simulations without pre-SN feedback become radiative at a lower radius than the ones with pre-SN feedback. The dips in the ratio are due to the SN energy reaching r_s earlier in the case with a cavity through pre-SN feedback than without pre-SN feedback. The order of magnitude difference in the energy feedback when including or not pre-SN feedback clearly highlight the importance of the pre-SN feedback.

Fig. 14 Cumulative total energy for stellar populations as a function of metallicity for different density and rs. The values are IMF-weighted (normalised to the total mass in stars above 8 M⊙), and a rotation velocity distribution for single and orbital period distribution for binary stars, and a binary fraction is assumed.

5 Discussion

5.1 Comparison to previous work

Previous studies have already examined the interaction between feedback from a single source and the ambient medium. For instance, Fierlinger et al. (2016) investigated the efficiency of the stellar wind and SN for a 60 M_⊙ stellar model with solar metallicity in an ambient density of 100 cm. Without pre-SN feedback, they find that the retained kinetic energy, measured when the velocity of the highest density cell has decreased to the sound speed of the ambient medium, is about 0.1 per cent of the input 10⁵¹ erg. Similarly, in our test without pre-SN wind at the same ambient density, albeit at a lower metallicity of Z = 0.004, we find retained kinetic energies between 0.04 per cent and 0.4 per cent (increasing with increasing SN energy) at 20 pc. When including stellar winds, the retained energy increases by a factor of around 3 at 5 pc, 19 at 10 pc and 27 at 20 pc for our highest mass model, which has a SN in the test case. Fierlinger et al. (2016) found an increase of the retained kinetic energy by a factor of six, which agrees with the range found in this work. However, the comparison of our results must be done with caution, as we also include photoionisation in our models.

Freyer et al. (2003) modelled the evolution of a uniform medium with number density 20 cm⁻³ around a 60 M_⊙ star at solar metallicity, affected by the feedback of wind and photoionisation. Due to a different treatment of the photoionisation energy, we recalculate their efficiency, the ratio of energy in their domain to the cumulative input, according to the method presented for the photoionisation energy in Sect. 4.2 by using their average stellar temperature and luminosity. With this method, the efficiency in Freyer et al. (2003) is around 1.4 × 10⁻² at the end of the stellar lifetime. The efficiency of the cumulative energy flow of our most similar model, log₁₀(M/M_⊙) = 1.8 without stellar rotation, is around 9.8 × 10⁻³ at a distance of 5 pc at the end of the stellar lifetime. This is lower but within a factor of two of Freyer et al. (2003). However, our simulation has a lower metallicity of Z = 0.004 and a higher ambient density of 100 cm⁻³. Fig. 12 displays a weak trend of higher efficiency for higher metallicities and lower densities, albeit for whole stellar populations with SNe. Considering this, we find no disagreement between the efficiency of our model and the one of Freyer et al. (2003).

5.2 Uncertainties and caveats

Despite the computational effort, uncertainties remain in the models used in this work. Among these are the stellar modelling (for a discussion see Fichtner et al. 2022), where for instance the rate of mass-loss is still under debate (e.g. Björklund et al. 2023; Sabhahit et al. 2022), the coupling of the stellar models into PION, the SN rate and energy, the composition of stellar populations and the 1D simulations with PION. Spherically symmetric simulations cannot capture turbulence or its dissipative effects through mixing hot and cool gas phases, which is especially effective at the wind-ISM interface (Mackey et al. 2015; Lancaster et al. 2021). Approximate methods have been applied to the expansion of superbubbles driven by SN feedback (El-Badry et al. 2019), but these should always be calibrated to 3D simulations and observations, and in addition the metallicity/magnetic-field dependence of such schemes is unclear.

Geen et al. (2015b) found that 3D simulations with uniform initial conditions result in a circumstellar medium with almost spherical symmetry and they do not expect large differences compared with 1D simulations with similar initial conditions, but a comparison of feedback efficiency for 3D vs. 1D was not carried out. We cannot compare our results meaningfully with Geen et al. (2015b) because of the different stellar mass and metallicity. The Geen et al. (2015b) simulations produced spherically symmetric structures because (i) the initial conditions were uniform with uniform external pressure and (ii) the 15 M_⊙ star in their calculations only has main sequence and red supergiant phases, during which strong dynamical instabilities do not develop for expansion into a uniform ISM. Arthur et al. (2011) showed that substruture is diminished for H II regions expanding around B stars compared to more massive O stars, because they emit a softer radiation field, producing a weaker temperature jump between ionised and neutral gas. While this and other diffusive processes may destroy substructure in circumstellar nebulae and favour a more spherically symmetric evolution (McKee et al. 1984; Arthur et al. 2011), there are also physical processes and environmental factors that drive the solution away from spherical symmetry. Simulations of a more massive star than used in Geen et al. (2015b) with wind-wind interactions from a WR phase would have resulted in the development of strong dynamical instabilities breaking spherical symmetry, as in Freyer et al. (2006). Further deviations from spherical symmetry can be caused, for example, by a structured ISM, the movement of the star providing feedback (e.g. Mac Low et al. 1991; Mackey et al. 2015), turbulent clouds (e.g. Medina et al. 2014), self-gravitating clouds (e.g. Geen et al. 2015a), stellar rotation (Chita et al. 2008) or the rotational and gravitational field of the cloud (e.g. Dalgleish et al. 2018).

In particular, the multiphase or stratified, non-uniform density structure of the ISM can lead to varied expansion rates and efficiencies of feedback bubbles, which could allow the feedback bubble to expand further during the stellar lifetime, resulting in larger radii than shown in Fig. 5. These effects of ISM structures and dynamical instabilities are expected to be stronger for higher mass stars with more intense ionising radiation fields (Arthur et al. 2011). This pre-SN bubble size is important for the radius at which the SN feedback efficiency drops. Fractal clouds with strong density stratification have been shown to increase the effectiveness of SN feedback (e.g. Walch & Naab 2015; Haid et al. 2018), and analytic models and simulations show that such structures also influence pre-SN circumstellar nebulae (McKee et al. 1984; Rogers & Pittard 2013; Haid et al. 2019). The growth of the H II region can be further influenced, for instance through champagne flows, leading to changes in the detailed evolution (e.g. of Fig. 5). H II region expansion into power-law density profiles can differ from expansion into a uniform medium through the development of ionisation-front instabilities (Garcia-Segura & Franco 1996). Inclusion of stellar winds further complicates the hydrodynamics of circumstellar material in 3D simulations (e.g. Geen et al. 2023b).

It is not possible to capture these inherently multidimensional effects with our 1D simulations because of the imposed symmetry. Future work making quantitative comparison between 1D and 3D solutions for more realistic environments would be valuable to assess the accuracy of our approach to subresolution feedback modelling. In this aspect, the presented results may be seen as lower limits of the possible feedback for r_s values higher than the pre-SN bubble radius. We note, however, that our 1D approach represents an effort to improve on current subgrid implementations of stellar feedback in simulations of galaxy evolution. These implementations are effectively ‘0D’ in the sense that stellar winds, radiation and SNe are used directly from stellar evolution calculations.

Furthermore, each of our 1D simulations includes only one source (a single star or a binary system), whereas star formation is often clustered. Clustering could alter the large-scale impact of pre-SN and SN feedback. Multiple stars in close proximity could drive a wind bubble and H II together, creating a larger cavity (Mac Low & McCray 1988; Rahner et al. 2017). This could allow more SNe to provide efficient feedback and the feedback energy to be retained until larger radii. Superbubbles created by the correlated feedback of stellar populations can escape out of the galactic disc (Mac Low & McCray 1988). This can launch galactic fountains and outflows, through which energy, mass and metals can be transported into the circumgalactic medium and beyond. We plan to use the results of our 1D calculations in a subgrid model for stellar feedback to investigate galactic winds in large-scale simulations of galaxy formation.

We do not consider SN Type Ia in this work. They originate from lower mass stars and explode with a large time delay after the stellar birth. Due to this, the explosion does not occur close to the birth place. Therefore, their evolution is not strongly influenced by the early feedback of massive stars and can be modelled separately.

Fig. 15 Ratio of cumulative total energy flow for a model without pre-SN feedback to the standard model with pre-SN feedback for single star models with Z = 0.004 and no initial rotation velocity. The values are both IMF-weighted (normalised to the total mass in stars above 8 M⊙) and the SNe are placed at the same time as in the standard model.

5.3 Consequences of our results for stellar feedback

In nearly all cases for the considered parameters, one star or binary system is able to create a cavity of at least 10 pc size. Similarly, Chevance et al. (2022) found that a populations of stars can disperse giant molecular clouds already in the first few Myr after the first massive stars form. The importance of the post-main-sequence phase is in agreement with the findings in the superbubble N 206 in the SMC. In this low-density region, one WO binary is completely dominating the pre-SN feedback budget over the 320 OB stars (Ramachandran et al. 2019).

Rotating and binary systems are thought to be important for galaxy formation due to their higher feedback. While the feedback energy from stellar winds of a population of rotating and binary stars of our grid was found to be around one order of magnitude higher than from single non-rotating stars at the lowest metallicity of Z = 0.0001, the global properties of a galaxy simulated with a feedback subgrid model based on these different populations were remarkable similar (Fichtner et al. 2022). Consistent with that, we find, in our 1D simulations of the circumstellar medium, similar energies (within a factor of two) for a population of single non-rotating stars, rotating single stars, binary stars and a population combining single and binary stars when evaluated for the lowest density at r_s = 5 pc. This suggests that the importance of rotating and binary stars might be overestimated. It is enough if the pre-SN feedback is able to create a cavity, and then the SN energy input is the main contributor to feedback energy reaching larger scales. Contrary to expectations, main-sequence winds of massive stars do not dominate the pre-SN feedback; the main contributors are post-main-sequence winds and photoionisation.

Populations of binaries, but also some populations of single stars with rotation velocities of 500 or 600 km s⁻¹, retain around one order of magnitude more feedback energy than single non-rotating models at r_s = 20 pc for the highest density (see e.g. Fig. 12). In these populations, the feedback of some stellar models is able to evacuate the star’s surrounding to at least 20 pc before the (secondary) SN occurs. For the population of binaries, the models with two SNe with a primary mass above 10^1.4 M_⊙ are of particular importance to create such large cavities. Therefore, the different strength of pre-(secondary)-SN feedback is mainly important through the size of the cavity. However, this is mainly relevant when measuring the energy at a r_s similar to the cavity size. Far inside the cavity the similar SN energies lead to similar cumulative energies flows, while far outside the energy dissipation by cooling dominates.

The feedback energy return is nearly independent of metallicity. At large radii even an opposite trend than expected can be seen: the energy feedback is higher at lower metallicity due to the weaker metal cooling. This trend might hold also true for lower metallicity. However, the question arises if in a metalfree, population III setting a similar behaviour is possible. To approximate such a situation, we rerun one simulation of a non-rotating 32 M_⊙ stellar model from the lowest metallicity without any cooling from metals and no stellar wind (Krtička & Kubát 2006, 2009). We find that the energy at r_s = 5 pc is lower with 1.2 × 10⁵¹ erg compared to the fiducial case with 1.7 × 10⁵¹ erg, likely due to the missing stellar winds. This is still higher than the 7.7 × 10⁵⁰ erg of the SN-only test case (which has, however, a higher metallicity and therefore another, but higher, SN energy). In contrast, at r_s = 20 pc, the feedback energy of the population III-test case is higher than the fiducial case, with 7.0 × 10⁵⁰ erg to 5.2 × 10⁵⁰ erg, while the SN-only test case has radiated away nearly all the energy. The feedback from photoionisation and SN in a metal-free ISM is able to create a larger feedback bubble than in the fiducial case, leading to the higher energy at 20 pc. The lower energy injection due to the lack of stellar winds is (over-)compensated by a lower energy dissipation due to less cooling. Therefore, even for population III stars the pre-SN feedback might be able to create a feedback bubble that is able to support a SN energy return to scales of tens of pc. This might reduce SF in the surroundings of the first stars (Mackey et al. 2003) and expel the gas from the low-mass host minihalo (Klessen & Glover 2023).

Our results show that the details of radiative and mechanical pre-SN feedback are not crucial to determining the overall feedback energy returned to the ISM, but the primary dependence is on which stars explode and with how much energy. We use the compactness parameter to determine explodability and explosion energy following O’Connor & Ott (2011) and Schneider et al. (2021), but recently Wang et al. (2022) found that their 2D SN explosion simulations disagree quite strongly with this prescription in some aspects. In particular most of their stars explode successfully regardless of the compactness parameter and, for the more massive stars, a compact core is actually required to successfully initiate an explosion. Burrows et al. (2024) have started a similar analysis of a suite of 3D CCSN simulations, finding similar results. Given the rapid development of this field and current lack of consensus in the literature, this must be seen as one of the largest uncertainties in determining feedback from stellar populations.

As well as these recent developments in SN modelling, the topic of winds from massive stars and the upper mass limit has received renewed attention (Vink 2022). In particular, as very massive stars approach the Eddington limit, their mass-loss rates should increase more strongly with mass (Sabhahit et al. 2022) than standard scaling laws for O and B stars. Higgins et al. (2023) showed that if the upper mass limit is on the order of 500 M_⊙ then the mass and metal yields from winds of very massive stars (>100 M_⊙) may contribute significantly to the total, even when integrated over an IMF (albeit at high metallicity). We have not included stars with mass M > 160 M_⊙ in this work, nor stellar models representative of binary mergers, but it will be interesting to see what effect such stars would have on the cumulative energy feedback rates presented here.

5.4 Implications for a subgrid feedback model

The use of 1D simulations of the interaction with the surrounding medium leads to changes compared to taking values directly from stellar models. The first difference is a lower feedback energy due to energy dissipation on scales that are not resolved. For example, in cosmological simulations this might lead to an overestimate of the effect of stellar feedback on a galactic scale. In addition, there is a time delay between the ejection of the stellar feedback and its arrival at distances further away from the star. Furthermore, the ambient medium alters the kinetic-to-total energy ratio, which decreases with increasing distance from the source.

As discussed earlier, we find nearly no metallicity dependence of the feedback energy. This has implications for high-redshift or low-mass galaxies which have a low metallicity, as it increases the impact of feedback on these systems. The inclusion of photoionisation and stellar winds, especially postmain-sequence winds, is important for feedback schemes, as they change the timing and efficiency of energy transfer to the ISM. The energy released by SNe varies among different populations, such as populations with lower metallicity. This is particularly important as the SN energy dominates the energy budget in the circumstellar medium.

We compare our results to the widely used Feedback In Realistic Environments (FIRE) project (Hopkins et al. 2023). To derive their feedback energy, we integrate their equations (1), (2), (4), and (5) in their Sect. 3.3 until 1 Gyr. We then renormalise the IMF stated in their Sect. 3.3 to include only massive stars above 8 M_⊙, to match the population definition used in this work. Our feedback energy is 30 to 100 times lower for the highest density and largest distances, due to energy dissipation in the ambient medium (see Fig. 16). Even at low density and distance, where we do not see strong energy dissipation, our energy feedback is a few times lower. This could originate from the differences in the treatment of CCSNe. In Hopkins et al. (2023), all massive stars (up to 120 M_⊙) are assumed to undergo a CCSN, while in our work only some models undergo a CCSN depending on the compactness parameter of the stellar model. The stars around 8–10 M_⊙, which are the most numerous, have the lowest-energy explosions when using our SN prescription, whereas Hopkins et al. (2023) assume that the explosion energy is independent of mass.

Subgrid models of SN feedback commonly inject energy only in one component, for example thermal energy, or keep a fixed ratio between different components. However, we unveil a much more complicated picture. The fraction of energy in the kinetic component is a complex function of time, evolutionary phase, metallicity, density and distance, which should be accounted for to create more realistic subgrid models of stellar feedback.

The density and resolution in galaxy formation simulations are connected to each other, in the sense that higher resolutions enable higher densities. For instance, in large-scale simulations like IllustrisTNG-100 gravitational softening lengths of the order of 100 pc with star formation thresholds of 0.1 cm⁻³ are reached (Pillepich et al. 2018b,a), while recently zoom-in simulations of Gutcke et al. (2022) use a gravitational softening of 0.5 pc and 4 pc for gas and stars, respectively and a star formation threshold of 1000 cm⁻³. For the available feedback budget, we find that a higher density and a larger measurement radius, r_s, have a similar effect: they both lead to more energy dissipation. Therefore, the impact of larger distances (i.e. lower resolution) in galaxy simulations could be reduced by the corresponding lower achieved densities and might help to maintain a higher feedback energy budget despite the lower resolution.

Implementation of our feedback prescription into detailed galaxy-formation simulations are required to assess the impact of the gas-altered feedback on a whole galaxy, which may depend on the resolution adopted in the galaxy formation simulation. For resolutions of pc-scale the effect will likely be smaller than for lower resolutions, because the energy dissipation and time delay are small. We plan to implement a subresolution model with the feedback from our 1D simulations to test against standard feedback implementations.

Fig. 16 Feedback energy from a stellar population as a function of metallicity for different density and radii is shown with blue cross symbols. The red dots indicate the energy feedback from CCSNe and stellar winds in FIRE-3 (Hopkins et al. 2023), while the orange dots include additionally photoionisation feedback that is assumed to have a coupling efficiency of 3 per cent. All values are IMF-weighted (normalised to the total mass in stars above 8 M⊙).

6 Summary

Simulations of galaxy formation rely on the inclusion of stellar feedback through subgrid models, where the feedback energy is often calculated directly from the stellar ejection. However, the circumstellar medium can alter the feedback – with respect to that injected by stars – on scales below the resolution of the simulations. In this work, we investigate the resolution gap between stellar and galaxy simulations by quantifying the energy that reaches pc to tens of pc scales.

We compute a large grid of simulations of the circumstellar medium with the code PION, including the feedback from stellar winds, photoionisation, and SNe of massive stars. We adopt an extended version of the detailed grid of stellar evolutionary models presented in Fichtner et al. (2022), which are computed with MESA. The grid includes different initial rotation velocities for single stars and orbital periods for binary systems. We use four different ambient densities and five different metallicities. For the SNe, we perform computations to identify which models undergo a CCSN according to the compactness parameter, a commonly used criterion for explodability. Additionally, we assume a relation between SN energy and the compactness parameter.

In the grid of PION simulations, we evaluate the thermal and kinetic energy flowing through spheres of three exemplary radii: 5, 10, and 20 pc. Assuming a standard Salpeter IMF for stars of 8 to 160 M⊙, our findings can be summarised as follows:

1. The use of the criterion for explodability, as opposed to the standard approach of injecting 10⁵¹ erg per massive star, leads to a metallicity-dependent SN energy return. The kinetic energy input of stellar winds decreases sharply with metallicity, especially for non- or slowly rotating single stars. Despite the metallicity dependence, nearly all the stars considered are capable of producing a cavity of at least 10 pc in radius, even at low metallicity;

2. The location of the shell of ISM material swept up by photoionisation feedback at the end of the stellar lifetime is independent of the stellar initial mass (except for the lowest mass models) in a population where all other parameters are the same;

3. The energy budget at all radii is dominated by the SN energy. Despite this, pre-SN feedback strongly increases the efficiency of SN feedback (mostly by more than a factor of 10) when comparing the fiducial simulations to test runs without pre-SN feedback. The amount of retained SN energy is strongly influenced by the size of the pre-SN feedback cavity. These findings highlight the importance of pre-SN feedback and the necessity of their inclusion. Although we treat each star in isolation in the 1D simulations, the effect of stellar clustering will increase the efficiency of energy transfer to the ISM because of the larger cavities formed;

4. The time when the feedback energy front crosses 20 pc is significantly later than the time it is ejected from the star. The length of this delay is significant with respect to the cloud-formation timescale. The direct injection of feedback at the time of the stellar ejection on a scale of tens of pc could therefore halt SF too early;

5. Stellar masses of around 30 M_⊙ make the greatest contribution to the total energy, as they are massive enough for strong winds, but still explode as CCSN, unlike more massive stars. Additionally, such stars are not as rare as the most massive stars;

6. Overall, the main sequence stellar winds are of limited importance for the energy budget. Before SN, the post-main sequence winds dominate the kinetic energy, and photoionisation dominates the thermal energy;

7. The strongest dependence of the energy feedback return is on the radius at which it is measured, with higher radii retaining lower energies. Additionally, a decrease in energy with increasing density is found, as expected. The feedback available for a subgrid model also therefore depends on the properties of the galaxy simulation (resolution, maximum resolved gas density), suggesting that feedback models should be developed that take a larger number of dependencies of the energy into account;

8. The total feedback energies are similar for single stars with different rotation velocities and binary systems with different orbital periods, although fast-rotating single stars have the highest feedback energies overall. Especially noticeable is the near independence of the feedback energy on metallicity, even though stellar winds are strongly metallicity dependent. This is explained by the similar SN energies across the grid and by the fact that almost all stellar models do produce some stellar wind bubble and H II region cavities that allow SNe to propagate to larger radii without losing significant energy. Binary systems with two SNe form bigger cavities, leading to higher retained energies even at 20 pc;

9. The efficiency of the feedback (the retained energy as a fraction of the input) is around a few per cent, except for the higher densities at 20 pc, where nearly all the energy is radiated away. The overall low efficiency reflects the low coupling efficiency of the photoionisation energy input. When taking the ratio without including the photoionsation energy in the input energy, this energy is nearly completely retained at 5 pc and around 100 cm⁻³ (and probably below);

10. The fraction of energy in the kinetic component is not a constant – as assumed in standard subgrid models – but depends on the time, ambient density, metallicity, radius, and the considered stellar track. When assuming a binary mass fraction of 70 per cent (with a rotation velocity and an orbital period distribution), the fraction of kinetic to total energy ranges between 0.66 and 0.95 close to the star, while further away the thermal energy becomes more important (kinetic fractions of 0.22 to 0.62) due to shock heating.

Despite its lower energy, pre-SN feedback is essential even at low metallicities (e.g. like those found in high-redshift galaxies) as it preprocesses the star’s surroundings before the succeeding SN. We do not find a strong decrease in the feedback energy with decreasing metallicity, because all calculated models produce sufficient pre-SN feedback that SNe explode into a low-density cavity and therefore remain in the adiabatic phase for longer. Our work shows that it is not the amount of pre-SN feedback energy that is of major importance, but rather the feedback mechanism implemented, and notably the existence of pre-SN feedback, binary systems with two SNe, and the timing of the feedback. The ambient medium influences, for example, the energy budget, feedback delay, and the kinetic energy fraction already at scales that are subresolution in galaxy formation simulations. These effects should be accounted for in order to create more realistic numerical models of stellar feedback.

The scripts to generate the PION simulations as well as the flow data of the shown simulations at the three discussed radii are publicly available.

Acknowledgements

YAF, ERD and CP gratefully acknowledge the Collaborative Research Center 1601 (SFB 1601 sub-projects C5 and C6) funded by the Deutsche Forschungsgemeinschaft (DFG, German Research Foundation) – 500700252. This work was also carried out within the Collaborative Research Centre 956, sub-project C04, funded by the Deutsche Forschungsgemeinschaft (DFG) – project ID 184018867. YAF is part of the International Max Planck research school in Astronomy and Astrophysics and guest at the Max-Planck-Institute for Radio Astronomy. JM acknowledges support from a Royal Society-Science Foundation Ireland University Research Fellowship (20/RS-URF-R/3712). JM is grateful to Prof. N. Langer for the use of a guest office at the AIfA during the time when this research was carried out. CP is grateful to SISSA, the University of Trieste, and IFPU, where part of this work was carried out, for hospitality and support. In addition to the software cited in the main text, this research made use of: the SILO library for data I/O (https://software.llnl.gov/Silo/), the SUNDIALS ODE solver library (Hindmarsh et al. 2005; Gardner et al. 2022), PYPION (Green & Mackey 2021), NUMPY (Harris et al. 2020), MATPLOTLIB (Hunter 2007), SCIPY (Virtanen et al. 2020), ASTROPY (Astropy Collaboration 2022) and SEABORN (Waskom 2021).

Data availability

The data underlying this article are available in Zenodo at https://doi.org/10.5281/zenodo.13285149.

Appendix A Evolution of a binary system

We show the profile around a binary system that has two SNe (from the primary and secondary star) during the evolution in Fig. A.1, in comparison to the single star case in Fig. 2. The mass-loss of the secondary star prevents the primary SN bubble from reaching as low densities as in the single case or after the secondary SNe (middle panel). Furthermore, the central region does not have the lowest density but rather the highest densities inside the bubble. As for the single star case, the secondary SN removes all structures inside the bubble (lower panel). The feedback bubble extends out to nearly 50 pc 5 Myr after the secondary SN, while in the case of the corresponding single star only 40 pc are affected 5 Myr after the SN. This implies that the volume of the cavity by the binary SNe is twice as big as in the single case.

Fig. A.1 Density, temperature and velocity profile of a PION simulation with a binary system of Z = 0.004, log10(P [days]) = 1.0 and log10(Mp [M⊙]) = 1.4, in a uniform medium of n = 100 cm−3 at the end of the primary lifetime (top), at the end of the secondary lifetime (middle) and 5 Myr after the secondary SN (bottom). Dashed lines indicate negative velocities.

References

All Tables

All Figures

Fig. 1 Explosion energy of each stellar model if it has a CCSN, where the colour of the squares indicates the SN energy. White squares indicate models that do not have a CCSN. The left column shows the single models, while the two right columns show the primary and secondary models of binary systems, respectively. Each row shows one metallicity. The x-axis shows the different names according to Table 1. The initial stellar masses are shown on the y-axis; we note that the binary grid has a coarser sampling in initial mass than the single grid. The highest initial masses are not shown as they have no CCSN.

In the text

Fig. 2 Density, temperature, and velocity profiles of a PION simulation with a single star model of Z = 0.004, υini = 0 km s−1, and log10(M/M⊙) = 1.5 in a uniform medium of n = 100 cm−3 at 5 Myr (top), at the end of the stellar lifetime (middle), and 5 Myr after the subsequent SN (bottom). Dashed lines indicate negative velocities. The circled numbers at the top indicate the five different regions found in the density profile due to pre-SN feedback.

In the text

Fig. 3 Cumulative energy flow through rs = 5, 10, and 20 pc for a single star model of Z = 0.004, υini = 0 km s−1, and log10(M/M⊙) = 1.5 in a uniform medium of n = 100 cm−3. The solid line shows the total energy, while the dashed and dotted lines show the kinetic and thermal energy component, respectively. A SN explosion occurs around 6.1 Myr. The star symbols (labelled ‘uniform’) show a comparison of the total energy with respect to a SN run in a uniform medium without pre-SN feedback beforehand.

In the text

	Fig. 4 Density profiles of our simulations with single star models of Z = 0.004, υini = 0 km s−1, in a uniform medium of n = 100 cm−3 at the end of the star’s lifetime. For clarity, we show only the profiles of very second model of this sequence of stellar models. The labels indicate the logarithm of the initial masses in M⊙.
In the text

	Fig. 5 Comparison of the radius of the outer edge of the shell of swept-up material to an estimate from the Spitzer radius. The radius is calculated from the single star models of Z = 0.004, υini = 0 km s−1, in a uniform medium of n = 100 cm−3 at the end of the star’s lifetime.
In the text

Fig. 6 Cumulative kinetic (top) and thermal (bottom) feedback energy for different stellar initial mass as a function of time for the single star models with Z = 0.004 and no initial rotation velocity is shown with coloured bars for the whole time of the PION simulation. The different columns show different initial ambient densities, while the lower three rows vary rs. The star symbol indicates the time of the SN, while black regions visualise outflowing energies below 1046 erg or inflows. The upper row shows IMF-weighted results (normalised to the total mass in stars above 8 M⊙). Additionally, the sum of the energy inputs from SN, stellar winds and an estimate of the photoionisation energy is shown (labelled ‘Star’), and the energy input excluding the photionisation energy (labelled ‘no ph’).

In the text

Fig. 7 IMF-weighted (normalised to the total mass in stars above 8 M⊙) cumulative energy flow for the single star models with metallicity Z = 0.004 and no initial rotation with a background density of 500 cm−3. Additionally, the sum of the energy inputs from SN, stellar winds and an estimate of the photoionisation energy is shown (labelled ‘Star’), and the energy input excluding the photionisation energy (labelled ‘no ph’).

In the text

	Fig. 8 Same as Fig. 6, but for binary systems with Z = 0.004 and log10(P [days]) = 0.6. The star symbol indicates the time of a SN (primary or secondary), while the diamond symbol indicates the end of the primary lifetime without a SN.
In the text

Fig. 9 IMF-weighted (normalised to the total mass in stars above 8 M⊙) cumulative energy flow for the single star models with metallicity Z = 0.0001 and no initial rotation. From top to bottom, the rows in each panel show kinetic, thermal, and total energy flow, while the columns show different initial ambient densities. Additionally, the sum of the energy inputs from SN, stellar winds and an estimate of the photoionisation energy is shown (labelled ‘Star’), and the energy input excluding the photionisation energy (labelled ‘no ph’).

In the text

	Fig. 10 Cumulative total energy flow for all models through rs = 10 pc, weighted by an IMF (normalised to the total mass in stars above 8 M⊙). The different subpanels show different combinations of metallicity and ambient density. The grey horizontal line separates single from binary models.
In the text

	Fig. 11 Ratio of cumulative kinetic energy flow through three rs values to the cumulative total energy flow at the final time. The values are IMF-weighted. Each subpanel shows all combinations of tracks and metallicity, while the different subpanels show different combinations of ambient density and radius.
In the text

	Fig. 12 Ratio of cumulative total energy flow through three rs values to the cumulative energy input through winds, SN, and photoionisation at the final time. The values are IMF-weighted. Each subpanel shows all combinations of tracks and metallicity, while the different subpanels show different combinations of ambient density and radius.
In the text

	Fig. 13 Same as Fig. 12, but the energy input only considers winds and SN, but no photoionisation. We note the different colour scale.
In the text

	Fig. 14 Cumulative total energy for stellar populations as a function of metallicity for different density and rs. The values are IMF-weighted (normalised to the total mass in stars above 8 M⊙), and a rotation velocity distribution for single and orbital period distribution for binary stars, and a binary fraction is assumed.
In the text

	Fig. 15 Ratio of cumulative total energy flow for a model without pre-SN feedback to the standard model with pre-SN feedback for single star models with Z = 0.004 and no initial rotation velocity. The values are both IMF-weighted (normalised to the total mass in stars above 8 M⊙) and the SNe are placed at the same time as in the standard model.
In the text

Fig. 16 Feedback energy from a stellar population as a function of metallicity for different density and radii is shown with blue cross symbols. The red dots indicate the energy feedback from CCSNe and stellar winds in FIRE-3 (Hopkins et al. 2023), while the orange dots include additionally photoionisation feedback that is assumed to have a coupling efficiency of 3 per cent. All values are IMF-weighted (normalised to the total mass in stars above 8 M⊙).

In the text

	Fig. A.1 Density, temperature and velocity profile of a PION simulation with a binary system of Z = 0.004, log10(P [days]) = 1.0 and log10(Mp [M⊙]) = 1.4, in a uniform medium of n = 100 cm−3 at the end of the primary lifetime (top), at the end of the secondary lifetime (middle) and 5 Myr after the secondary SN (bottom). Dashed lines indicate negative velocities.
In the text