Resolving Star Cluster Formation in Galaxy Simulations with Cosmic Ray Feedback

TL;DR

This study addresses how cosmic-ray feedback influences star-cluster formation in galaxies by leveraging high-resolution tallbox AREPO simulations with dynamically coupled cosmic rays and a multiphase ISM. The authors compare three CR transport implementations within the Crisp gas physics framework, identifying star clusters with a 4D Friends-of-Friends method and analyzing their mass functions, environmental dependence, and SN clustering. They find that CRs reduce the star-formation rate and shift cluster demographics within the range allowed by observations; while cluster radii and velocity dispersions show modest changes, the virial parameters trend lower in CR-enabled runs, indicating more bound clusters due to less turbulent support. Overall, the CRISP approach demonstrates predictive power in linking CR transport physics to star-cluster demographics, while the idealized setup highlights areas for future realism, such as radiative feedback and global galactic dynamics.

Abstract

Star clusters host the massive stars responsible for feedback in star-forming galaxies. Stellar feedback shapes the interstellar medium (ISM), affecting the formation of future star clusters. To self-consistently capture the interplay between feedback and star formation, a model must resolve the parsec-scale star formation sites and the multiphase ISM. Additionally, the dynamical impact of cosmic rays (CRs) on star formation rates (SFRs) must also be considered. We present the first simulations of the formation of an ensemble of star clusters with dynamically-coupled CRs, near-individual star particles, and a feedback-regulated ISM. We analyze tallbox simulations performed using the CRISP model in the moving-mesh code AREPO. We apply varied implementations of CR transport under the theory of self-confinement. We find that CRs simultaneously reduce the SFR, the power law slope of the cluster mass function, and the cluster formation efficiency. Each simulation is compatible with observations, and CR feedback tends to move results along observed star cluster relations. We see only modest changes in cluster radius and velocity dispersions, but significant differences in the virial parameters. Ultimately, the primary impact of CRs is to reduce SFRs. Lower SFRs imply fewer supernovae, and consequently a lower turbulent energy budget for gas. Star clusters formed in a CR-regulated ISM have lower velocity dispersions, and are therefore more bound under self-gravity. The effective clustering of SNe is unchanged by CRs. Through this work, we demonstrate the predictive power of the CRISP feedback model, despite this idealized setup.

Figures (9)

• Figure 1: Star particles in the CR-NL-IN simulation at $t\approx250\;\mathrm{Myr}$, colored by the age of the star particle. Bluer particles are younger, redder particles are older. Three filled contours of molecular hydrogen surface density are displayed in gray, with star-forming gas cells marked with a purple $\times$. We see a variety of clusters of star particles with similar ages, suggesting these stars formed together and remained bound. The youngest star clusters are primarily localized near molecular clouds and star-forming regions. In contrast, the oldest stars and star clusters appear throughout the full domain.
• Figure 2: Star cluster mass function using the young cluster sample described in Section \ref{['sec:cluster_identification']}, where the star clusters are identified with an age of $\sim10\;\mathrm{Myr}$. The mass functions are fit with power law using a Poisson likelihood, as described in Section \ref{['sec:mass_function_sub']}. The best fits are displayed as a dashed line, and the power law slope for each fit is reported in the legend. We neglect star clusters below a mass of $10^{3}\;\mathrm{M}_\odot$ in the fit, as the mass function appears to change slope at masses less than this threshold. We also include a reference line with power law slope $-2$. Each run produces a cluster mass function with a power law slope between $-2$ and $-2.5$, consistent with observationally-determined slopes in environments with $\Sigma_\mathrm{SFR}$ values similar to this work adamo_probing_2015johnson_panchromatic_2017adamo_star_2020messa_young_2018. Our results are reported in Table \ref{['tab:data_table_appendix']}.
• Figure 3: The best fit and error for the power law slope of the mass function within an SFRD bin for the three cases, from this work and relevant observational works. The errors here are determined by bootstrapping. The power law slopes inferred from the simulations show a positive correlation with the SFRD and are broadly consistent with the observational values and trends. Observed values are from M31 johnson_panchromatic_2017, M51 messa_young_2018, M83 adamo_probing_2015, the HiPEEC sample adamo_star_2020-1, and the LEGUS sample adamo_star_2020. Our results are reported in Table \ref{['tab:data_table_appendix']}.
• Figure 4: CFE ($\Gamma$) vs. SFRD range ($\Sigma_\mathrm{SFR}$) for the three cases and three SFRD subsamples. The CFE is defined as the percent of stellar mass that is formed in massive bound clusters ($M_\mathrm{bound}>1000\;\mathrm{M}_\odot$, squares; $M_\mathrm{bound}>5000\;\mathrm{M}_\odot$, circles) relative to the total stellar mass formed in the SFRD bin. The bound cluster criteria is described in Section \ref{['sec:cfe_subsub']}. The values of $\Gamma$ from our simulations with the lower mass threshold of $5000\;\mathrm{M}_\odot$ are an excellent match to the theoretical prediction from kruijssen_fraction_2012 and broadly agree with the data from observations. Observed values, in order of appearance in the legend, are from adamo_probing_2011, adamo_probing_2015, annibali_cluster_2011, bastian_star_2008, fensch_massive_2019, goddard_fraction_2010, ginsburg_high_2018, johnson_panchromatic_2016, larsen_luminosity_2002, lim_star_2015, messa_young_2018, pasquali_infrared_2011, ryon_snapshot_2014, rafelski_star_2005, and silva-villa_star_2011. Fits are from goddard_fraction_2010, kruijssen_fraction_2012, and chandar_fraction_2017. Labels with a $\dagger$ symbol in the legend have SFRD values determined from their respective references and $\Gamma$ values from goddard_fraction_2010. Our results are reported in Table \ref{['tab:data_table_appendix']}.
• Figure 5: Upper panel: Distribution of ISM volume as a function of density ($n_\mathrm{H}$). The sample includes gas within $200\;\mathrm{pc}$ of the midplane over the time range $112\;\mathrm{Myr}$ to $250\;\mathrm{Myr}$. Lower panel: Distribution of SN injection environment densities, over the same time range. The lightly shaded regions are the fraction that occur less than $3\;\mathrm{Myr}$ after $t_{60}$ (the approximate beginning of a cluster's formation). For both distributions, the fraction occurring below densities $10^{-2}\;\mathrm{cm}^{-3}$ are included as a percentage. Despite the significant differences in ISM gas distributions, approximately half of SNe occur in diffuse gas environments, irrespective of the CR model. Additionally, the first few SNe tend to occur in very dense gas, tracing the conditions of star particle formation.