An observational test for star formation prescriptions in cosmological hydrodynamical simulations

TL;DR

This study uses observed star formation clustering to evaluate and constrain the density threshold parameter in cosmological simulations, finding that high thresholds better match observed clustering patterns.

Contribution

It introduces a method to use star formation clustering measurements to refine sub-grid star formation models in cosmological simulations.

Findings

Observed clustering favors high star formation density thresholds (>10 cm$^{-3}$).

Lower thresholds (<1 cm$^{-3}$) are inconsistent with observed clustering.

Clustering measurements can effectively constrain star formation prescriptions.

Abstract

State-of-the-art cosmological hydrodynamical simulations of galaxy formation have reached the point at which their outcomes result in galaxies with ever more realism. Still, the employed sub-grid models include several free parameters such as the density threshold, $n$, to localize the star-forming gas. In this work, we investigate the possibilities to utilize the observed clustered nature of star formation (SF) in order to refine SF prescriptions and constrain the density threshold parameter. To this end, we measure the clustering strength, correlation length and power-law index of the two-point correlation function of young ($\tau<50$ Myr) stellar particles and compare our results to observations from the HST Legacy Extragalactic UV Survey (LEGUS). Our simulations reveal a clear trend of larger clustering signal and power-law index and lower correlation length as the SF threshold increases with only mild dependence on galaxy properties such as stellar mass or specific star formation rate. In conclusion, we find that the observed clustering of SF is inconsistent with a low threshold for SF ($n<1$ cm$^{-3}$) and strongly favours a high value for the density threshold of SF ($n>10$ cm$^{-3}$), as for example employed in the NIHAO project.

Figures

Tables

Table 1. Table 1: Clustering parameters for young stellar particles: For each threshold run we show from top to bottom the three fit parameters, scale radius, r s subscript 𝑟 s r_{\rm s} , power-law index γ 𝛾 \gamma and normalization, A 𝐴 A , as well as derived values for the radius r ​ ( ξ = 0 ) 𝑟 𝜉 0 r(\xi=0) where the clustering strength goes to zero and the value of the clustering strength measured at a radius of r = 100 𝑟 100 r=100 pc, 1 + ξ ​ ( r = 100 ​ p ​ c ) 1 𝜉 𝑟 100 p c 1+\xi(r=100\rm{pc}) . From left to right we show the median with uncertainty ranges, 16 th and 84 th percentile and the sample size. Note: The sample size varies because some snapshots are discarded due to too few young stellar particles.

parameter	median	$p^{16th}$	$p^{84th}$	$N_{\mathrm{gal}}$
$n=20.0$
$r_s$ [kpc]	${3.55}_{-0.17}^{+0.71}$	1.67	11.62	189
$\gamma$	${0.074}_{-0.011}^{+0.010}$	-0.051	0.192	189
$A$	${1.57}_{-0.06}^{+0.15}$	0.89	3.24	189
$r\left(\xi =0\right)$ [kpc]	${1.38}_{-0.12}^{+0.10}$	1.13	2.53	189
$1+\xi \left(r=100pc\right)$	${2.19}_{-0.07}^{+0.13}$	1.32	3.68	189
$n=10.0$
$r_s$ [kpc]	${4.77}_{-0.26}^{+0.54}$	2.13	10.13	148
$\gamma$	${0.094}_{-0.007}^{+0.012}$	0.028	0.213	148
$A$	${1.31}_{-0.03}^{+0.07}$	0.97	1.96	148
$r\left(\xi =0\right)$ [kpc]	${1.84}_{-0.09}^{+0.10}$	0.97	2.87	148
$1+\xi \left(r=100pc\right)$	${2.07}_{-0.05}^{+0.13}$	1.53	3.38	148
$n=1.0$
$r_s$ [kpc]	${15.45}_{-0.81}^{+2.24}$	6.84	39.25	168
$\gamma$	${0.032}_{-0.003}^{+0.007}$	0.003	0.106	168
$A$	${1.17}_{-0.02}^{+0.06}$	0.92	1.82	168
$r\left(\xi =0\right)$ [kpc]	${2.69}_{-0.19}^{+0.35}$	0.67	6.36	168
$1+\xi \left(r=100pc\right)$	${1.45}_{-0.02}^{+0.07}$	1.23	2.15	168
$n=0.1$
$r_s$ [kpc]	${10.79}_{-0.34}^{+0.85}$	7.05	20.17	180
$\gamma$	${0.027}_{-0.003}^{+0.004}$	-0.009	0.076	180
$A$	${1.35}_{-0.02}^{+0.06}$	1.10	1.98	180
$r\left(\xi =0\right)$ [kpc]	${3.07}_{-0.16}^{+0.21}$	1.27	5.44	180
$1+\xi \left(r=100pc\right)$	${1.63}_{-0.03}^{+0.05}$	1.32	2.16	180
observations
$\gamma$	$0.21\pm 0.02$
$r\left(\xi =0\right)$ [kpc]	$1.44\pm 0.83$
$1+\xi \left(r=100pc\right)$	$2.13\pm 0.23$

Equations

$$\xi (r)=\frac{DD(r)}{RR(r)}-1$$

$$1+\xi (r)=A{\left(\frac{r}{r_{\mathrm{s}}}\right)}^{-\gamma }\exp \left(-\frac{r}{r_{\mathrm{s}}}\right)$$

Full text

††pagerange: An observational test for star formation prescriptions in cosmological hydrodynamical simulations–References††pubyear: 2019

1 Introduction

It is widely assumed that stars primarily form in clusters across the stellar disks of galaxies. Because current cosmological simulations barely resolve the sites and processes of star formation (SF) itself, most state-of-the-art simulations implement phenomenological recipes for SF on scales of giant molecular clouds (GMCs). Although little is known about the small scale details of the SF process, such as the collapse and onset of GMC fragmentation or the impact of turbulent and magnetic pressure on the gas clouds, the larger scale galactic effects are well studied and understood (Kennicutt & Evans, 2012). E.g. the typical masses of GMCs, the sites of SF, are $\sim 10^{2}-10^{4}{\rm M}_{\odot}$ with typical densities of $\sim 10^{2}-10^{5}$ cm*-3*.

In order to identify analogues of molecular clouds in state-of-the-art hydrodynamical simulations several recipes have been introduced (see e.g. Hopkins et al., 2013). Due to the huge dynamical ranges involved in the formation of galaxies (from kilo-parsec scales down to sub-parsec scales) several key physical ingredients need to be modelled as sub-grid models with free parameters. One of these sub-grid models deals with the formation of stars out of the interstellar gas where SF is usually localized via a density threshold $n$. This threshold varies over more than 5 orders of magnitude across different simulations from $n=0.01-10^{3}$ cm*-3*. Only very few works focused on the impact of different localization criteria for SF on the global galaxy properties (Hopkins et al., 2013; Dutton et al., 2019a) and even fewer works compared the detailed outcome of the adopted SF criteria to the observed properties of star forming regions in galaxies.

Although current cosmological simulations lack the resolution to resolve the details of SF inside the molecular clouds, the large scale effects and properties of these clouds are well resolved in the highest resolution simulations today, e.g. APOSTLE (Sawala et al., 2016), AURIGA (Grand et al., 2017), FIRE (Hopkins et al., 2018) or NIHAO (Buck et al., 2019) and especially the observed spatial clustering behaviour of young stellar complexes (Grasha et al., 2015, 2017) should be well resembled by current models of SF.

Here we set out to introduce and describe a conclusive test for SF models in cosmological simulations by means of the clustering signal of young star clusters (in both simulations and observations). We use the results on clustered SF to put constraints on (up to now) free parameters of current cosmological simulations of galaxy formation. This paper is structured as follows: in §2 we describe the simulation suite, in §3 we lay out the methodology used in this work, in §4 we present our constraints on the SF parameter and finally in §5 we discuss our results and present our conclusions.

2 Cosmological Simulations

For this work we use 17 sets of simulations with stellar masses from 10^{7}\mathrel{\hbox{\hbox to0.0pt{\hbox{\lower 4.0pt\hbox{\sim}}\hss}\hbox{<}}}M_{\rm star}\mathrel{\hbox{\hbox to0.0pt{\hbox{\lower 4.0pt\hbox{\sim}}\hss}\hbox{<}}}10^{11} ${\rm M}_{\odot}$ . For each simulation we analyze the most massive galaxy at 13 snapshots equally spaced in time between redshift $z=0$ and $z=0.5$.

These simulations are a sub-sample of the NIHAO hydrodynamical cosmological zoom-in suite (Wang et al., 2015) run with the smooth-particle-hydrodynamics (SPH) code Gasoline2 (Wadsley et al., 2017). The galaxy formation model is described in detail in the papers of Wang et al. (2015) and Dutton et al. (2019a). The most important parameter to our study is the SF threshold, $n$, which we describe below.

The fiducial NIHAO simulations adopt as SF threshold of $n=10$ cm${}^{-3}\sim 50m_{\rm gas}/\epsilon_{\rm gas}^{3}$ where 50 denotes the number of neighbours in the SPH smoothing kernel, $m_{\rm gas}$ the initial gas particle mass, and $\epsilon_{\rm gas}$ the gravitational force softening Stinson et al. (2006, 2013). In NIHAO $\epsilon_{\rm gas}\propto m_{\rm gas}^{1/3}$, so that $n$ is independent of gas particle mass. Each fiducial simulation is complemented by three additional SF threshold runs: $n=0.1$, $n=1.0$ and $n=20.0$. The former two are further described and used in Dutton et al. (2019a). The threshold of $n=0.1$ is similar to the values adopted by the EAGLE/APOSTLE (Schaye et al., 2015; Sawala et al., 2016) and ILLUSTRIS/AURIGA (Vogelsberger et al., 2014; Grand et al., 2017) projects111Note: All values of $n$ used here are well below typical densities of GMCs ($n\sim 100$ cm*-3*) and their SF cores ($n\sim 10^{4}$ cm*-3*).. While the SF threshold is not the only parameter than can be varied, in this study we focus on this particular parameter given its strong impact on galaxy morphology and halo response (e.g. Dutton et al., 2019a; Benitez-Llambay et al., 2018; Bose et al., 2018). We like to highlight that the $n=20.0$, $n=10.0$ and $n=1.0$ simulations are run with exactly the same parameters except for the star formation threshold (see Dutton et al., 2019a, for more details) and still follow the same abundance matching relation across cosmic time although with very different morphology for the young stars. The $n=0.1$ runs were recalibrated by reducing the efficiency of the early stellar feedback from $e=0.13$ to $0.04$ in order to make the galaxies fit the abundance matching relation again.

Haloes in the zoom-in regions are identified using the halo finder AHF2 (Knollmann & Knebe, 2009; Gill et al., 2004). The virial mass, $M_{200}$, denotes the mass of all particles within a sphere of virial radius, $R_{200}$, containing $\Delta$ = 200 times the cosmic critical matter density, $\rho_{\rm crit}$. In general, galaxy properties are measured from all hydro particles within $0.2R_{200}$. Dark matter particle masses and force softenings are chosen to resolve the mass profile at \mathrel{\hbox{\hbox to0.0pt{\hbox{\lower 4.0pt\hbox{\sim}}\hss}\hbox{<}}}1\% of $R_{200}$ with dark matter force softenings of $\epsilon_{\rm dark}=207-931$ pc, and hydro softenings of $\epsilon_{\rm gas}=88-400$ pc. Gas particle masses range between $3.5\times 10^{3}{\rm M}_{\odot}-3.2\times 10^{6}{\rm M}_{\odot}$222The results presented in this work have been compared against the simulations with 8 times higher mass resolution presented in Buck et al. (2019) and we found good convergence of the clustering behaviour of young stars.. All simulations in the NIHAO project, including the ones used here, employ a pressure floor to keep the Jeans mass of the gas resolved and suppress artificial fragmentation (see also appendix A1 of Smith et al., 2018). Our implementation follows Agertz et al. (2009) which is equivalent to the criteria proposed in Richings & Schaye (2016) and fulfils the Truelove et al. (1997) criterion at all times. Thus, the Jeans mass in our simulations is resolved with $\sim 4$ SPH kernel masses.

NIHAO galaxies provide a perfect test sample to study the effects of the SF parameters on the structure of galaxies since all simulations used here (runs for $n=0.1,1.0,10.0$) follow the abundance matching relations since $z=4$ (see Dutton et al., 2019a). They form the right amount of stars (Wang et al., 2015) and they follow the cold gas versus stellar mass relation at $z=0$ (Stinson et al., 2015; Buck et al., 2017). Furthermore, they are in very good agreement with the Local Group stellar mass function (Buck et al., 2019) as well as the velocity function (Macciò et al., 2016; Dutton et al., 2019b) and form MW analogues matching the the detailed kinematics of the MW bulge (Buck et al., 2018).

3 Method: Two-Point Correlation function of young star clusters

We use the two-point correlation function (TPCF), $\xi(r)$, to measure the magnitude of clustering of young stellar particles as a function of separation $r$ in a suite of cosmological hydrodynamical simulations. The two-dimensional correlation function $1+\xi(r)$ is defined as the probability above Poisson noise of finding two stellar particles with a separation $r$ as $dP=N^{2}[1+\xi(r)]d\mathbf{x_{1}}d\mathbf{x_{2}}$, where N is the surface density of young stellar particles within two infinitesimal spatial elements $d\mathbf{x_{1}}$ and $d\mathbf{x_{2}}$, separated by distance $r$ (see e.g. Peebles, 1980, for a similar definition in spherical coordinates). For a truly random Poisson distribution, the TPCF will be flat across all separations, such that $1+\xi(r)=1$. On the other hand, a clustered stellar distribution will have $1+\xi(r)>1$ at small values of $r$ and will be declining with increasing $r$ towards that of a flat, non-clustered distribution. The correlation function of a fractal (self-similar) distribution is described with a single power-law as $1+\xi(r)=(r/r_{\rm s})^{-\gamma}$, where $r_{\rm s}$ is the correlation length of the clustering and $\gamma$ describes the hierarchical ordering (Calzetti et al., 1989; Larson, 1995).

In a two-dimensional self-similar distribution, the total number of clusters $N$ increases with radius $r$ as $N\propto r^{D_{2}}$, where $D_{2}$ is the two-dimensional fractal dimension (Mandelbrot, 1982). The number of clusters for every radial aperture will increase as $N\propto r^{-\gamma}\times r^{2}\propto r^{-\gamma+2}$. Thus, the power-law slope $\gamma$ of the correlation function $1+\xi(r)$ determines the two-dimensional fractal dimension as $D_{2}=-\gamma+2$. A flat, non-clustered distribution of $\gamma=0$ will result in a fractal geometric dimension of $D_{2}=2$; a steep slope will indicate a clustered distribution with a fractal dimension less than 2. Indeed, the distribution of SF and interstellar gas is shown to exhibit a projected, two-dimensional fractal dimension of $D_{2}\sim 1.2-1.6$ over a large range of environments (e.g. Sánchez et al., 2005; Elmegreen et al., 2006; Sánchez & Alfaro, 2008).

We implement the calculation of the TPCF using the standard estimator formula:

[TABLE]

Here $DD$ denotes the number of data-data pairs and $RR$ the number of random-random pairs with separation between $r$ and $r+\delta r$. From this we determine the slope and the correlation length of the TPCF by fitting a power-law with an exponential cut-off to account for the finite extension of the stellar disc in the form of:

[TABLE]

$A$ is a normalization parameter, $r_{\rm s}$ the scale radius and $\gamma$ the power law index.

Recent measurements of the clustering of young stellar complexes are obtained from the HST Legacy Extragalactic UV Survey (Calzetti et al., 2015; Adamo et al., 2017; Messa et al., 2018, LEGUS). Grasha et al. (2017) measure the TPCF of young stellar complexes at scales from $5$ pc to $10$ kpc and fit it with a double power law measuring the correlation length, clustering strength and the power law index333Note: We define $\gamma$ without the minus sign which leads to a sign flip between our results and the ones of Grasha et al. (2017). below and above the break radius of $r\sim 100$ pc.

So in what follows we use their average power-law fit (from their table 4) at large radii in the form of $1+\xi(r)=5.6(\pm 0.3)r^{-0.21\pm 0.02}$ to compare with NIHAO. We use the results from Grasha et al. (2017) for young star clusters with ages of \mathrel{\hbox{\hbox to0.0pt{\hbox{\lower 4.0pt\hbox{\sim}}\hss}\hbox{<}}}40 Myr and sizes of roughly \mathrel{\hbox{\hbox to0.0pt{\hbox{\lower 4.0pt\hbox{\sim}}\hss}\hbox{<}}}300 pc (their classes 1, 2, 3) which is well in agreement with the spatial extent of stellar tracer particles in our cosmological simulations. From the fits of Grasha et al. (2017) and the fits to our simulation data we derive the clustering strength at a scale of $r=100$ pc, $1+\xi(r=100\rm{pc})$, and the correlation length as the radius where the correlation function $\xi$ goes to zero, $r(\xi=0)$ and compare the two results. A complete list of fitted and derived clustering parameters for the simulations and observations is given in table 1.

4 Results: Constraints on the star formation threshold

In figure 1 we show three versions of the same galaxy simulated with different density threshold values increasing from left to right. Open white circles show the positions of young stars ($\tau<50$ Myr) and color-coding shows the surface density of neutral hydrogen gas computed following the results from Rahmati et al. (2013) as implemented by Gutcke et al. (2017). Their total number as well as their number in the inner 5 kpc is indicated in each panel and is roughly equal for all three runs. We see that the high threshold run forms stars that are more clustered compared to the lower threshold runs. Especially the equal number of young stars within the central region shows that this not due simply more SF in the very center of the galaxy. This figure confirms the naive expectation that a higher threshold might lead to a more clustered SF compared to a low threshold. In the following section we analyse the effects of the SF threshold on the spatial distribution of young stars by fitting equation 2 to the TPCF of young stellar particles derived from the simulations. An example of this power law fit to the TPCF derived from one of our galaxies is shown in figure 2 with the fit parameters and their uncertainties shown in the legend. Typical fitting uncertainties on the derived parameters ($r_{\rm s}$ and $\gamma$) are on the order of \mathrel{\hbox{\hbox to0.0pt{\hbox{\lower 4.0pt\hbox{\sim}}\hss}\hbox{<}}}20\% for $\gamma$ and \mathrel{\hbox{\hbox to0.0pt{\hbox{\lower 4.0pt\hbox{\sim}}\hss}\hbox{<}}}5\% for $r_{\rm s}$. A stronger correlation of young stars in the high threshold run can be appreciated from the shorter correlation length, $r_{\rm s}$, and the overall stronger correlation signal compared to the other two threshold runs.

In figure 3 we compare the distribution of the clustering strength $1+\xi(r=100\rm{pc})$, the power-law index $\gamma$ and the correlation length $r(\xi=0)$ for our four sets of simulations. The $n=0.1$ runs are shown with green histograms, the $n=1.0$ runs in blue, the $n=10.0$ ones with red histograms and the $n=20.0$ runs with orange histograms. Vertical coloured lines indicate the median values and gray vertical bands show the observations. There is a clear trend of larger clustering strength and power-law index and lower correlation length with increasing SF threshold. While for all simulations the power-law index $\gamma$ is lower than observed, the clustering strength and the correlation length of the high ($n=10$) SF threshold runs are in agreement with the observed results. This shows that indeed SF is more clustered in the high threshold simulations.

4.1 Clustered star formation

We show the results for the correlation length, $r(\xi=0)$, and the clustering strength, $1+\xi(r=100\rm{pc})$, derived from the fits of Grasha et al. (2017) for a set of 6 LEGUS galaxies with gray stars or gray shaded bands in the figures.

In figure 4 we compare the correlation length (upper panels) and the clustering strength (lower panels) as a function of stellar mass (left panels) and specific SFR (sSFR, right panels) to observational results. We show with gray stars the results for individual galaxies where present (upper panels) and with gray bands the averaged observational results for the six LEGUS galaxies as presented in Grasha et al. (2017, lower panels). For our simulations we further show a linear regression to test for any correlation and display the correlation coefficients in the legend. In general, we recover the findings from figure 3 that the two lower threshold runs exhibit both larger correlation length and smaller clustering strength. We find a weak correlation of $r(\xi=0)$ with stellar mass and a weak anti-correlation of $1+\xi(r=100\rm{pc})$ with sSFR. We find that $r(\xi=0)$ derived from the simulations is slightly larger compared to the observed values of $\sim 1$ kpc but consistent with the observed values (gray stars) for the largest stellar masses and the highest threshold runs (red). On the other hand, for the $n=10.0$ runs $1+\xi(r=100\rm{pc})$ is in agreement with the observed range of $\sim 1.90-2.36$ (gray bands). The $n=0.1$ (green) and $n=1.0$ (blue) runs are inconsistent with the observed values for both $r(\xi=0)$ and $1+\xi(r=100\rm{pc})$ although we note there is a large scatter.

These findings are further supported by figure 5 which plots the correlation length $r(\xi=0)$ vs. the clustering strength, $1+\xi(r=100\rm{pc})$. For the simulations we show median values over all galaxies and redshifts. The error bars show the uncertainty on the median defined as: $\sigma_{\rm high}=1.22(p_{84}-p_{50})/\sqrt{N_{\rm gal}}$ and $\sigma_{\rm low}=1.22(p_{50}-p_{16})/\sqrt{N_{\rm gal}}$ respectively. Here $p_{84}$, $p_{16}$ and $p_{50}$ denote the $84^{\rm th}$, $16^{\rm th}$ percentile and the median. $N_{\rm gal}$ denotes the number of galaxies. The observed parameter range is shown with a gray box. The $n=0.1$ and $n=1.0$ runs are located in the upper left corner of this plot while the $n=10.0$ and $n=20.0$ runs show smaller correlation length and larger clustering strength, respectively. This indicates that only the high threshold runs ($n=10.0$ and $n=20.0$) result in a clustering of young stellar particles which is consistent with the observed clustering of young star clusters. This diagram shows a trend of increasing clustering of SF with increasing density threshold $n$ and thus advocates for the usage of a high density threshold in simulations of galaxy formation; much higher than what is used in the AURIGA or APOSTLE simulations.

5 Discussion and Conclusion

State-of-the-art cosmological hydrodynamical simulations of galaxy formation have reached the point at which their outcomes result in galaxies with ever more realistic properties (e.g. Buck et al., 2019; Grand et al., 2017). Still, the employed models are rather simplistic and involve several free parameters. One of these parameters is the density threshold $n$ used to localize the SF gas in the simulation. This parameter varies over several orders of magnitude among different simulation groups.

In this work we set out to investigate the prospects of using the observed clustered nature of SF in order to refine the SF prescriptions in simulations and constrain the density threshold parameter. We find that the observed strong clustering of SF only agrees with simulations that employ a high threshold for star formation (n\lower 2.15277pt\hbox{;\buildrel>\over{\sim};}10) like e.g. the NIHAO (Wang et al., 2015) and FIRE (Hopkins et al., 2018) simulations, the super-bubble feedback runs of Keller et al. (2014) or simulations using an H2 based SF prescription which naturally leads to a high density for star forming gas (e.g. Christensen et al., 2012; Munshi et al., 2019). Such high SF thresholds runs have been used to simulate realistic bulgeless galaxies (Governato et al., 2010) and dwarf galaxy populations of Milky Way mass galaxies (Zolotov et al., 2012; Brooks & Zolotov, 2014). On the other hand, observations disfavour simulations with a low density threshold like the AURIGA/ILLUSTRIS (Grand et al., 2017; Vogelsberger et al., 2014; Pillepich et al., 2018), APOSTLE/EAGLE (Sawala et al., 2016; Schaye et al., 2015), ROMULUS (Tremmel et al., 2017) or VELA (Zolotov et al., 2015) simulations. A lower value of the star formation threshold, $n$, results in a significant decrease in simulation runtime but comes at the cost of a less realistic star formation model. Furthermore, $n$ strongly impacts the predicted dark matter halo structure. As recently shown by Dutton et al. (2019a) a high value of $n$ naturally resolves the Too-Big-Too-Fail and the cusp-core tension of $\Lambda\textrm{CDM}$. On the other hand, this implies that a low value of $n$ does not lead to accurate predictions for the structure of $\Lambda\textrm{CDM}$ haloes.

In this study we focussed on varying the star formation threshold keeping other parameters like e.g. the star formation efficiency and the number of SPH kernel masses in the pressure floor fixed. Thus, it remains to be seen whether the clustering of the young stars for low star formation thresholds can be increased by choosing different numerical values for these parameters. However, blindly varying parameters does not teach us something about the underlying physics of galaxy formation. We thus see the presented method as an observational test of future star formation prescriptions in cosmological simulations. Such models have to recover not only the star formation histories and the total stellar masses (integrated star formation histories) but also the observed spatial distribution of star formation itself.

Finally, while it is a priori evident that a simple threshold value must be a too simplistic description of the SF process, our current computational resources prohibit more complex first principle prescriptions, at least for fully cosmological simulations. Improved models for isolated galaxy simulations without cosmological context exist (Hu et al., 2017; Semenov et al., 2017; Emerick et al., 2018) and updated models for cosmological simulations are explored (Applebaum et al., 2018). Therefore, we believe that the presented correlations between spatial clustering of young stellar complexes and the density threshold can be used to gauge the value of $n$ in future simulations and remove one free parameter from the models. Furthermore, given the current successes in galaxy formation models and the wealth of observational data obtained from large scale surveys the onus is now to compare model predictions and observations beyond simple integral properties such as total stellar mass or galaxy size and use the full range of galactic morphologies to test numerical and theoretical models. Only in this way we will improve our understanding of galaxy formation.

Acknowledgments

We thank the referee for an insightful report which helped us improve the manuscript. TB acknowledges support from the Sonderforschungsbereich SFB 881 “The Milky Way System” (subproject A2) of the German Research Foundation (DFG) and by the European Research Council under ERC-CoG grant CRAGSMAN-646955. The authors gratefully acknowledge the Gauss Centre for Supercomputing e.V. (www.gauss-centre.eu) for funding this project by providing computing time on the GCS Supercomputer SuperMUC at Leibniz Supercomputing Centre (www.lrz.de). This research was carried out on the High Performance Computing resources at New York University Abu Dhabi; Simulations have been performed on the ISAAC cluster of the Max-Planck-Institut für Astronomie and the HYDRA and DRACO clusters at the Rechenzentrum in Garching. This research made use of the pynbody Pontzen et al. (2013), tangos Pontzen & Tremmel (2018), matplotlib (Hunter, 2007), SciPy (Jones et al., 01) and NumPy, IPython and Jupyter (Walt et al., 2011; Pérez & Granger, 2007; Kluyver et al., 2016) python packages.