Energy equipartition between stellar and dark matter particles in cosmological simulations results in spurious growth of galaxy sizes

TL;DR

This paper demonstrates that in cosmological simulations, the mismatch in particle mass between stars and dark matter causes artificial energy transfer, leading to exaggerated galaxy sizes, which impacts the interpretation of simulation results.

Contribution

It reveals how mass segregation due to unequal particle masses causes spurious galaxy size growth in simulations, highlighting a key source of numerical artifact.

Findings

Galaxy sizes increase systematically when dark matter particles are more massive than stellar particles.

Size growth correlates with the ratio of stellar to dark matter particle mass.

Spurious effects are significant when the stellar-to-dark matter particle mass ratio exceeds one.

Abstract

The impact of 2-body scattering on the innermost density profiles of dark matter haloes is well established. We use a suite of cosmological simulations and idealised numerical experiments to show that 2-body scattering is exacerbated in situations where there are two species of unequal mass. This is a consequence of mass segregation and reflects a flow of kinetic energy from the more to less massive particles. This has important implications for the interpretation of galaxy sizes in cosmological hydrodynamic simulations, which nearly always model stars with less massive particles than are used for the dark matter. We compare idealised models as well as simulations from the EAGLE project that differ only in the mass resolution of the dark matter component, but keep sub-grid physics, baryonic mass resolution and gravitational force softening fixed. If the dark matter particle mass exceeds the mass of stellar particles, then galaxy sizes--quantified by their projected half-mass radii, ${\rm R_{50}}$--increase systematically with time until ${\rm R_{50}}$ exceeds a small fraction of the redshift-dependent mean inter-particle separation, $l$ (${\rm R_{50}}\geq 0.05\times l$). Our conclusions should also apply to simulations that adopt different hydrodynamic solvers, subgrid physics or adaptive softening, but in that case may need quantitative revision. Any simulation employing a stellar-to-dark matter particle mass ratio greater than unity will escalate spurious energy transfer from dark matter to baryons on small scales.

Figures

Tables

Table 1. Table 1: Basic numerical parameters used for our cosmological simulations. N 1 subscript N 1 {\rm N_{1}} refers to the number of “dark matter” particles of species 1; N 2 subscript N 2 {\rm N_{2}} to the number of “gas” particles, species 2. The corresponding particles masses are m 1 subscript 𝑚 1 m_{1} and m 2 subscript 𝑚 2 m_{2} , respectively, and their ratio is denoted μ ≡ m 1 / m 2 𝜇 subscript 𝑚 1 subscript 𝑚 2 \mu\equiv m_{1}/m_{2} ; L box subscript L box {\rm L_{box}} is the simulation box size. The run type is also provided: DMO refers to runs for which both species were assumed to be collisionless particles, and “ eagle ” refers to hydrodynamical runs carried out using the Reference model of the eagle project.

	$N_1$	$N_2$	$m_1$	$m_2$	$\mu$	$L_{\mathrm{box}}$	Type
			[${10}^5{M}_{\odot }$]	[${10}^5{M}_{\odot }$]		[Mpc]
	7523	0	1.8	0	–	12.5	DMO
	1883	1883	97.0	18.1	5.36	12.5	DMO
	7$\times$1883	1883	13.9	18.1	0.77	12.5	DMO
	3763	3763	97.0	18.1	5.36	25.0	eagle
	7$\times$3763	3763	13.9	18.1	0.77	25.0	eagle

Equations

$$\frac{\Delta v^2}{v^2}=\frac{8}{{\mathrm{N}}_1}\ln \Lambda \frac{(1+\psi /\mu )}{(1+\psi {)}^2},$$

$$\kappa_{\rm rel}\equiv\frac{t_{\rm rel}}{t_{\rm H}}=\frac{{\rm N_{1}}}{8\ln{\rm N_{1}}}\frac{t_{\rm orb}}{t_{\rm H}}=\frac{\sqrt{200}}{8}\frac{{\rm N_{1}}}{\ln{\rm N_{1}}}\biggr{(}\frac{\overline{\rho}\ }{\rho_{\rm crit}}\biggl{)}^{-1/2},$$

$$t_{\mathrm{seg}}=\frac{t_{\mathrm{rel}}}{\mu }\approx \frac{{\mathrm{N}}_1}{8\ln {\mathrm{N}}_1}\frac{t_{\mathrm{orb}}}{\mu },$$

Full text

11footnotetext: E-mail: aaron.ludlow@icrar.org

1 Introduction

Cosmological simulations of collisionless dark matter (DM) make reliable predictions for the innermost structure of DM haloes. Such simulations incur relatively modest computational cost and have been repeated at ever increasing resolution, exposing the limits of their reliability (see, e.g., Stadel et al., 2009; Navarro et al., 2010). Controlling for other numerical parameters–such as time-stepping, integration accuracy and gravitational softening–their main impediment is 2-body relaxation, which sets a lower-limit to the spatial resolution of any N-body simulation (Power et al. 2003; hereafter P03; Ludlow et al. 2018; hereafter LSB18). This limitation is well understood and readily accounted for, leading to widespread agreement on the innermost structure of DM haloes.

Such simulations provide the rudimentary infrastructure for modelling galaxy formation, offering a tangible connection to observational astrophysics. Current approaches to this problem follow semi-analytic or halo occupation methods–here the physics of galaxy formation is divorced from the evolution of DM–or simultaneously model the co-evolution of DM and baryonic fluids. In both approaches, sub-resolution models for galaxy formation require calibration against observables before sensible predictions for galaxy populations can be made. This may overshadow the complex non-linear coupling between numerical and subgrid parameters, and may mask subtle numerical effects.

One possible issue–which we highlight in this letter–is the importance of 2-body relaxation for the stellar component of simulated galaxies. Stars are treated as collisionless particles in cosmological simulations and, like DM, their dynamics must be subject to 2-body scattering. Galaxies formed in cosmological simulations, while calibrated to resemble observed systems, may therefore evolve in a way that is subject to numerical artefact.

In Section 2 we discuss the importance of 2-body scattering in N-body simulations, emphasising differences between uniform resolution runs and those involving mixtures of DM and stars of unequal mass, which is the conventional approach. We present simple numerical experiments that illustrate the effects. In Section 3.1 we describe the cosmological simulations used to test the impact of 2-body scattering on the evolution of stellar systems; their results are presented in Sections 3.2 and 3.3. We provide some closing remarks in Section 4.

2 2-body relaxation in an idealised galaxy-halo model

Cosmological simulations involve mixtures gas, stars and DM particles typically of unequal mass. When collisions cannot be ignored, their evolution is subject to 2-body scattering and, when masses are unequal, to energy equipartition (e.g. Spitzer & Hart, 1971). The net energy exchange between species due to these processes can be described by a diffusion equation, with coefficients that depend on their initial phase-space distributions, and the ratio of particle masses.

Following Binney & Tremaine (2008), we consider the collisional relaxation time of such a system, neglecting the gas component. We define the particle mass ratio, $\mu\equiv m_{1}/m_{2}\geq 1$, and the fraction of mass in $m_{2}$: $\psi\equiv{\rm M_{2}/M_{1}}={\rm N_{2}}\,m_{2}/{\rm N_{1}}\,m_{1}$, where ${\rm N}_{i}$ are the number of particles of species $i$. A test particle that traverses a system of size $R$ will experience $\delta n=\delta n_{1}+\delta n_{2}\approx 2\,\pi\,(\Sigma_{1}+\Sigma_{2})\,b\,db$ collisions with impact parameters in the range $(b,b+db)$, where $\Sigma_{1}={\rm N_{1}}/\pi\,R^{2}$ and $\Sigma_{2}=\psi\,\mu\,{\rm N_{1}}/\pi\,R^{2}$ are the surface densities of species 1 and 2, respectively. From the impulse approximation, any single encounter results in a small velocity perturbation ($\delta v\ll v$) perpendicular to the particle’s direction of motion; its trajectory is unaltered. Regardless of its mass, velocity perturbations are of order $|\delta v_{i}|\approx 2\,G\,m_{i}/(b\,v)$ for encounters with particles of mass $m_{i}$. Such encounters add incoherently and their cumulative effect will be given by integrating $\delta v_{1}^{2}\delta n_{1}+\delta v_{2}^{2}\delta n_{2}$ over some range of impact parameters, $b_{\rm min}$ to $b_{\rm max}$. The relative square velocity change after traversing the system is given by

[TABLE]

where we have assumed a typical velocity $v^{2}\approx G\,{\rm N_{1}}\,m_{1}\,(1+\psi)/R$, and $\Lambda\equiv b_{\rm max}/b_{\rm min}$ is the Coulomb logarithm.

For cosmological simulations eq. 1 can be simplified if we identify species 1 with DM and species 2 with stars; $\psi$ is then the stellar-to-DM halo mass ratio, typically \lower 2.15277pt\hbox{;\buildrel<\over{\sim};}0.05. Assuming equal numbers of baryon and DM particles $\mu=(\Omega_{\rm M}-\Omega_{\rm bar})/\Omega_{\rm bar}$, where $\Omega_{\rm M}$ and $\Omega_{\rm bar}$ are the cosmic densities of matter and baryons, respectively. In this case $\mu\geq 1$ and $\psi\ll 1$, and the ratio of bracketed terms in eq. 1 is close to unity and may be ignored. If we further assume $b_{\rm max}=R$ and set $b_{\rm min}=b_{90}=G\,(m_{1}+m_{2})/v^{2}$ as the impact parameter yielding 90∘ deflections, then $\Lambda={\rm N_{1}}\,(1+\psi)/(1+\mu^{-1})\approx{\rm N_{1}}$ and eq. 1 reduces to $\Delta v^{2}/v^{2}\approx 8\ln{\rm N_{1}}/{\rm N_{1}}$. The relaxation of both species is driven by encounters with massive particles.

The number of orbits a particle must complete so that $\Delta v^{2}/v^{2}\approx 1$ defines the relaxation time, $t_{\rm rel}=t_{\rm orb}/(\Delta v^{2}/v^{2})$. In units of the Hubble time (roughly the orbital time at the radius***We define $r_{200}$ as the size of a sphere centred on the particle with the minimum potential energy that encloses a mean density of $200\times\rho_{\rm crit}$, where $\rho_{\rm crit}(z)=3\,H(z)^{2}/(8\,\pi\,G)$ is the critical density. The corresponding mass is $M_{200}=(800/3)\,\pi\,r_{200}^{3}\rho_{\rm crit}$ and circular velocity, $V_{200}=\sqrt{G\,M_{200}/r_{200}}$., $r_{200}$), $t_{\rm H}\approx 2\,\pi\,r_{200}/{\rm V_{200}}$, this can be expressed

[TABLE]

where $N\equiv N(r)$ is the enclosed particle number, $\overline{\rho}(r)$ the enclosed density, and $t_{\rm orb}=2\,\pi\,r/V$ is the local orbital time†††A more common definition of the relaxation time is based on the number of crossings a particle must execute such that $\Delta v^{2}/v^{2}\approx 1$, which differs from our definition by a factor $t_{\rm orb}/t_{\rm cross}=\pi$. We adopt the orbital time to define $t_{\rm rel}$ for consistency with P03: in this case, $t_{\rm orb}/t_{\rm H}\approx(r/V)/(r_{200}/V_{200})$ results in eq. 2. (P03). When other numerical parameters are chosen wisely, $t_{\rm rel}$ sets a minimum resolved spatial scale within which collisions cannot be ignored. The solution to eq. 2 thus defines a “convergence radius”, $r_{\rm conv}$, which marks the location at which $\kappa_{\rm rel}\sim 1$ (see, e.g., P03; LSB18).

The value of $\kappa_{\rm rel}$ corresponding to a certain level of convergence must be obtained empirically by comparing simulations of differing mass resolution. P03 found that, for DM-only simulations, the circular velocity profile, $V_{\rm c}(r)$, of an individual Milky Way-mass halo converges to $\approx 10$ per cent at the radius where $\kappa\approx 0.6$; similar convergence in the average $V_{\rm c}(r)$ profiles requires a less conservative value, $\kappa\approx 0.18$, regardless of halo mass (LSB18). A convenient approximation is given by $r_{\rm conv}=0.174\,\kappa_{\rm rel}^{2/3}\,l$, where $l={\rm L_{box}}/N_{\rm part}^{1/3}$ is the mean inter-particle spacing in physical units, and $\kappa_{\rm rel}=0.18$ (LSB18).

When $\mu\neq 1$, two-body collisions also lead to a segregation of the two components: massive particles will, on average, lose energy to less massive ones, causing them to congregate in halo centres while heating the low-mass component. This “mass segregation” signals the onset of energy equipartition.

The simple 2-component toy model of Spitzer (1969) suggests that the segregation timescale, $t_{\rm seg}$, is shorter than $t_{\rm rel}$ by a factor roughly equal to the ratio of the particle masses:

[TABLE]

Homogeneous mixtures of particles of different mass will therefore segregate at radii $r_{\rm seg}\geq r_{\rm conv}$ provided $\mu\geq 1$. A simple estimate of $r_{\rm seg}$ therefore follows from eq. 2 (or from $r_{\rm conv}=0.174\,\kappa_{\rm rel}^{2/3}\,l$) if $\kappa_{\rm rel}$ is replaced by $\kappa_{\rm seg}=\mu\,\kappa_{\rm rel}$. Whether equipartition can be reached, however, depends on the ratios of particle mass, $\mu$, and of the total mass of each component, ${\rm M_{1}/M_{2}}$, and $r_{\rm seg}$ should therefore be viewed as an upper limit. (Simple analytic estimates and numerical results suggest that full equipartition may not be possible if {\rm M_{1}}\lower 2.15277pt\hbox{;\buildrel>\over{\sim};}{\rm M_{2}}\,\mu^{-3/2}, which is almost always the case in DM-dominated galaxies.)

As $\mu\rightarrow 1$, the importance of mass segregation diminishes. Nevertheless, different species may still structurally evolve through 2-body scattering and, as we show below, this evolution is sensitive to the initial segregation of each component (The sensitivity to initial segregation becomes clear when comparing the idealised tests presented in Figure 1 to the collisionless cosmological runs in Figure 2, which start with both species equally mixed.).

Figure 1 shows results from numerical experiments designed to illustrate these effects. We consider here idealised equilibrium systems composed of a galaxy embedded within a DM halo. Both are modelled as spherical, collisionless Hernquist 1990 spheres with galaxy-to-halo mass ratio $M_{\rm gal}/M_{\rm h}=0.027$ (close to the “peak” galaxy formation efficiency of Behroozi et al. 2013) and ratio of scale radii $r_{\rm half}/a_{\rm h}=0.25$ ($r_{\rm half}$ and $a_{\rm h}$ are the galaxy half-mass radius and halo scale radius, respectively). Initial conditions, constructed using GalIC (Yurin & Springel, 2014), differ only in the stellar-to-DM particle mass ratio, $\mu$. GalIC assigns particle phase-space coordinates iteratively in order to optimally match an underlying analytic distribution function, which is itself a stable solution to the collisionless Boltzmann equation and therefore in collisionless equilibrium‡‡‡We have verified that GalICs yields stable particle configurations for our experiments by carrying out tests for which $\mu=1$ (to eliminate mass segregation) and ${\rm N_{1}}=10^{6}$ (to suppress 2-body scattering). Our tests confirm that the mass profiles of species 1 and 2 remain stable at r\lower 2.15277pt\hbox{;\buildrel>\over{\sim};}r_{\rm conv} for a Hubble time.. We adopt ${\rm N_{1}}=5\times 10^{4}$ (for DM), and consider $\mu=1$, 2, 5 and 25. All runs used the same softening length, $\epsilon/l_{h}=0.1$ ($l_{h}=[3/4\,\pi\,N_{1}]^{1/3}$ is the Wigner-Seitz radius), and were evolved using gadget-2 (Springel, 2005) for $t\approx 13.3$ Gyr. (We have verified that our simulation results are robust to small changes in timestepping and softening length.) Because these systems are initially in collisionless equilibrium, any evolution away from the initial state must be driven by 2-body scattering.

Different panels correspond to different $\mu$, as indicated. Solid blue curves show $V_{c,1}(r)$ for the DM, and solid red curves show $V_{c,2}(r)$ for stars; tints and shades encode the time evolution, which increases linearly from $t=0$ (light) to $t\approx 13.3\,{\rm Gyr}$ (dark). Dashed lines of corresponding colour show the initial profiles used to construct the galaxy/halo models. For each curve (except the initial profiles) thick lines extend down to the convergence radius expected from eq. 2 (for $\kappa_{\rm rel}=0.6$); thin lines extend to the radius enclosing 100 DM particles.

DM profiles are reasonably stable for r\lower 2.15277pt\hbox{;\buildrel>\over{\sim};}r_{\rm conv}, as is the case for stars if $\mu=1$. Note, however, that as $\mu$ increases, the curves deviate systematically from their initial profile at radii \lower 2.15277pt\hbox{;\buildrel>\over{\sim};}r_{\rm conv}; this is particularly true for stars. The arrows mark $r_{\rm seg}$ calculated from eq. 2 after replacing $t_{\rm rel}$ by $t_{\rm seg}=t_{\rm rel}/\mu$ . For \mu\lower 2.15277pt\hbox{;\buildrel<\over{\sim};}5, these arrows track more closely the radii at which $V_{c}(r)$ profiles first show noticeable differences from their initial values. Note also that the segregation of the stars and DM is much more prominent when $\mu$ is large: DM haloes develop denser centres while the stellar component gradually expands. Importantly, even for $\mu=1$ there is considerable evolution in ${\rm V_{c,2}}(r)$ for r\lower 2.15277pt\hbox{;\buildrel<\over{\sim};}r_{\rm conv}. This is because 2-body collisions will tend to homogenise populations that are initially segregated. 2-body scattering can be thought of as a diffusion process with different coefficients describing the first- and second-order processes. For particular initial configurations–which depend on $\mu$, $M_{1}/M_{2}$ and the spatial segregation of each species–energy equipartition may give rise to mass segregation. For the particular case of $\mu=1$ the diffusion coefficients are equal, and 2-body scattering will lead to a mixing of the two components. In that case, a centrally compact stellar component will tend to become more diffuse with time, as seen in the upper left panel of Figure 1, even if it was constructed to be in collisionless equilibrium initially (note that scattering-driven diffusion is largely confined to within the convergence radius, as expected for $\mu=1$).

3 Cosmological Simulations

3.1 Simulation set-up

DM haloes and their associated galaxies form hierarchically through accretion and mergers and are, at best, *quasi-*equilibrium structures. It is therefore worthwhile to test the importance of mass segregation and 2-body scattering in cosmological simulations that include two particle species. The remainder of the paper will focus on such simulations.

All cosmological runs used parameters consistent with the Planck Collaboration et al. (2014) data release: $h\equiv H_{0}/(100\,{\rm km\,s^{-1}\,Mpc^{-1}})=0.6777$ is the Hubble parameter; $\sigma_{8}=0.8288$ the ($z=0$) rms density fluctuation in 8 $h^{-1}{\rm Mpc}$ spheres; and $\Omega_{\rm M}=1-\Omega_{\Lambda}=0.307$ and $\Omega_{\rm bar}=\Omega_{\rm M}-\Omega_{\rm DM}=0.04825$, are the energy density parameters in units of $\rho_{\rm crit}$.

We use three cosmological simulations that echo those used by Binney & Knebe (2002) to investigate 2-body scattering in cosmological DM-only simulations. The first evolves the DM with ${\rm N_{part}=752^{3}}$ equal-mass particles ($m_{\rm DM}=1.8\times 10^{5}\,{\rm M_{\odot}}$). The second uses two particle species of equal abundance, ${\rm N_{1}=N_{2}=188^{3}}$, but with a mass ratio $\mu=\Omega_{\rm DM}/\Omega_{\rm bar}\approx 5.36$; this run is analogous to most cosmological hydrodynamical simulations (DM and baryons are sampled with equal particle numbers) but differs by modelling both species as collisionless fluids. The masses of DM (species 1) and “gas” particles (species 2) are, respectively, $m_{\rm DM}=97.0\times 10^{5}\,{\rm M_{\odot}}$ and $m_{\rm gas}=18.1\times 10^{5}\,{\rm M_{\odot}}$. The final run also adopts a two collisionless components, but with unequal particle numbers: ${\rm N_{1}}/7={\rm N_{2}}=188^{3}$ and hence $\mu=(1/7)\,\Omega_{\rm DM}/\Omega_{\rm bar}\approx 0.77$ (or $m_{\rm DM}=13.9\times 10^{5}\,{\rm M_{\odot}}$). Initial conditions were created from those described above by splitting DM particles into seven equal-mass particles and arranging them on a cubic lattice in a manner that preserves a force-free unperturbed particle load. All runs used a linear box size of $L_{\rm box}=12.5$ Mpc (comoving) and identical phases and amplitudes for mutually resolved modes. They differ only in the number of particles of species 1, and hence in $\mu$.

The particle mixture models were repeated for the second set of simulations (in a larger volume; ${\rm L_{box}}=25$ Mpc) but with species 2 treated as a gaseous fluid ($N_{\rm gas}=N_{2}=376^{3}$). These runs employ cooling, star formation and feedback from stars and active galactic nuclei in accord with the Reference model of the eagle program (see Schaye et al., 2015). They differ only in the number of DM particles: one has $N_{\rm DM}=N_{\rm gas}=376^{3}$ ($\mu\approx 5.36$) and the other $N_{\rm DM}=7\times N_{\rm gas}=7\times 376^{3}$ ($\mu=0.77$). At $z=0$, star particles have an average mass of roughly 65 percent of the primordial gas particle mass (stars lose mass to gas particles through stellar winds as they evolve). As a result, $\mu_{\star}\equiv m_{\rm DM}/\langle m_{\star}\rangle\approx 1.18$ in our runs, but we neglect this small departure from unity.

All runs used the same softening length for both species, which is a fixed fraction of the mean baryonic inter-particle separation§§§For $N_{\rm p}=752^{3}$, we use the DM inter-particle spacing.: $\epsilon/l=0.04$ (comoving) for $z>2.8$, and $\epsilon/l=0.01$ (physical) thereafter. Haloes were identified using subfind (Springel et al., 2001), which returns the coordinate of the particle with the minimum potential energy, $\mathbf{x}_{\rm MB}$, as well as $r_{200}$, ${\rm M_{200}}$ and ${\rm V_{200}}$.

As in Figure 1, we focus our analysis on the circular velocity profiles of each mass component, and use subscripts to denote the relevant species. (For example, $V_{c,1}(r)$ refers to the circular velocity profile of particles of mass $m_{1}$.) Hereafter, for clarity, we drop explicit reference to DM or baryonic particles, even in the hydrodynamic runs, but instead identify DM with species 1 and stars with species 2. We do not consider the mass profiles of gas particles in our hydrodynamic simulations.

3.2 Cosmological simulations with unequal-mass collisionless particles

Figure 2 shows the median ($z=0$) circular velocity profiles of haloes in four separate mass bins in our collisionless cosmological runs. Grey curves correspond to the uniform resolution (${\rm N_{part}}=752^{3}$) simulation, which can be used to assess convergence in the lower-resolution runs. Blue curves correspond to the run with $N_{1}=N_{2}=188^{3}$ ($\mu=5.36$); orange to the one with $N_{1}/7=N_{2}=188^{3}$ ($\mu=0.77$). Thick lines extend down to $r_{\rm conv}=0.055\,l$ (LSB18), where $l$ is the mean inter-particle spacing¶¶¶If we were to use instead $l={\rm L_{box}/N_{tot}}$ , where ${\rm N_{tot}=N_{1}+N_{2}}$, $r_{\rm conv}$ would be smaller by a factor of $2^{1/3}\approx 1.26$. This will not affect the interpretation of our results, so we opt for the more conservative estimate of $r_{\rm conv}$. of particles of mass $m_{1}$; thin lines to the $r_{\rm conv}$ expected from eq. 2 with $\kappa_{\rm rel}=0.6$. To aid the comparison, all curves have been normalised to $V_{0}=\sqrt{G\,M_{200,i}/r_{200}}$, where $M_{200,i}$ is the mass of species $i$ enclosed by $r_{200}$.

This figure prompts a few comments. First, notice that, for simulations involving particle mixtures, the $V_{c,1}(r)$ profiles agree reasonably well with those of equal-mass haloes in the uniform resolution run (the solid coloured curves align closely with the grey curves). The largest differences in $V_{c}(r)$ are \lower 2.15277pt\hbox{;\buildrel<\over{\sim};}10 per cent for $r>r_{\rm conv}$, as expected. Particles of mass $m_{2}$, however, behave differently depending on $\mu$. For $\mu\approx 0.77$ (dashed orange lines), the circular velocity profiles of species 1 and 2 are similar: both deviate by \lower 2.15277pt\hbox{;\buildrel<\over{\sim};}10 per cent from the high-resolution run for all $r>r_{\rm conv}$ and all halo masses considered. This is expected: since $\mu$ is close to 1, and both species are initially well-mixed, it follows that $t_{\rm seg}\approx t_{\rm rel}$ and $r_{\rm seg}\approx r_{\rm conv}$, and both species should remain approximately homogeneous at r\lower 2.15277pt\hbox{;\buildrel>\over{\sim};}r_{\rm conv} at all times. For $\mu\approx 5.36$, however, this is not the case: for r\lower 2.15277pt\hbox{;\buildrel<\over{\sim};}r_{\rm seg}, $V_{\rm c,2}(r)$ is considerably lower than what is expected for a purely collisionless system, consistent with mass segregation driven by 2-body scattering (this confirms the results of Binney & Knebe 2002). The radius below which this suppression becomes significant coincides roughly with $r_{\rm seg}$ (downward pointing arrows), approximated by $r_{\rm seg}=0.174\,(\mu\,\kappa_{\rm rel})^{2/3}\,l\approx 3.1\,r_{\rm conv}$ (assuming $\kappa_{\rm rel}=0.188$; LSB18).

3.3 Impact of 2-body scattering on galaxy sizes

Many cosmological simulations spawn one star particle per gas particle, which typically have comparable masses but are $\approx\Omega_{\rm bar}/(\Omega_{\rm M}-\Omega_{\rm bar})$ times less massive than the DM particles. Other simulations attempt to increase the resolution of the stellar component by generating multiple star particles per gas particle which are considerably less massive (e.g. Dubois et al., 2014; Revaz & Jablonka, 2018). Galaxies formed in both types of simulations may be subject to equipartition effects, which may have important implications for the interpretation of galaxy sizes, among other properties.

What impact does equipartition have on galaxy sizes in cosmological hydrodynamical simulations? Figure 3 summarises the results of our tests. Each panel shows the median projected half-stellar mass radii, ${\rm R_{50}}$, as a function of galaxy stellar mass (masses are defined using bound stellar particles within a 100 physical kpc aperture centred on $\mathbf{x}_{\rm MB}$) at four different redshifts: $z=0$, 0.5, 1 and 2. We use blue curves for $\mu=5.36$ and orange curves for $\mu=0.77$. The vertical dashed lines correspond to 100 primordial gas particles, dotted lines to 2000. These runs use identical baryonic mass resolution, force softening (arrows indicate $2.8\times\epsilon$) and sub-grid physics models; they differ only in DM particle mass.

Galaxy sizes show clear differences between these runs, both in their mass and redshift dependence. Consider first $z=0$ (upper-left panel). For $\mu=5.36$, the median size-mass relation flattens abruptly for stellar masses {\rm M}_{\star}\lower 2.15277pt\hbox{;\buildrel<\over{\sim};}2000\,m_{\rm gas} (dotted vertical line) below which ${\rm R}_{50}\approx 2.8\,{\rm kpc}$, regardless of ${\rm M_{\star}}$. For $\mu\approx 0.77$ this is not observed: sizes continue to decrease monotonically with decreasing ${\rm M}_{\star}$ to the lowest mass-scale considered ($\approx 10$ stellar particles). Similar results are seen at $z=0.5$ for $\mu=5.36$, although in this case ${\rm R}_{50}$ levels-off at lower mass (${\rm M}_{\star}\approx 10^{8.7}\,{\rm M_{\odot}}$), and correspondingly smaller size (${\rm R}_{50}\approx 2\,{\rm kpc}$). For $\mu=0.77$ galaxy sizes evolve very little from $z=0.5$ to $z=0$ (thin lines, repeated in all panels, show the $z=0$ size-mass relations for comparison).

Note as well that, for the different $\mu$ values, sizes begin to converge at higher redshift: by $z=2$, for example, they are virtually indistinguishable for galaxies resolved with more than $\approx 100$ particles. Intriguingly, convergence is attained at all $z$ provided sizes exceed the physical convergence radius of haloes in the $\mu=5.36$ run (shown here as $r_{\rm conv}=0.055\,l$ and highlighted using a red horizontal line; see LSB18).

Although using $\mu\approx 1$ will minimise the spurious transport of energy between particle species, we emphasise that by itself it does not guarantee that the simulations are immune to numerical effects. Convergence tests that simultaneously increase both the DM and baryonic resolution, and use $\mu\approx 1$, are required to test in which regime the results are robust.

4 Summary and Discussion

Previous studies of galaxy sizes in cosmological simulations report trends similar to those in Figure 3 for $\mu=5.36$. In eagle, Furlong et al. (2017) note that galaxy sizes increase systematically with increasing ${\rm M_{\star}}$ and with decreasing redshift. They identified a sample of passive galaxies between $z=1.5$ and 2 that remain quiescent centrals until $z=0$: all grow in size between their identification redshift and $z=0$. They report that compact centrals at $z=2$ grow secularly by “stellar migration” to the present day.

Campbell et al. (2017) present convergence tests of projected half-mass radii in the Apostle simulations (Sawala et al., 2016, $\mu\approx 5.36$). Comparing low-, intermediate- and high-resolution runs they show that ${\rm R_{50}}$ flattens at a characteristic scale comparable to the spline softening length. Our results indicate that sizes are subject to spurious growth below scales comparable to the convergence radius ($\approx 0.055\,l$) which are close to those quoted by Campbell et al. (2017). We can distinguish between softening and 2-body scattering as the culprit for this resolution dependence using the redshift evolution of ${\rm R_{50}}$. If softening is the cause, then ${\rm R_{50}}\approx\epsilon$ should set the minimum size at all redshifts, whereas 2-body scattering would give rise to a slow growth of ${\rm R_{50}}$ for poorly-resolved galaxies. Our results support the latter explanation (Figure 3).

Similar results were recently reported for galaxy sizes in the Illustris TNG50 simulation (Pillepich et al., 2019). In TNG50, for which $\mu\approx 5.3$, the sizes of low-mass galaxies flatten at systematically larger physical scales, and at higher stellar masses, as mass-resolution decreases. These results are not consistent with softening setting a minimum physical size to low-mass galaxies. In TNG100, Genel et al. (2018) also report a flattening of sizes for low-mass galaxies, an effect that becomes more pronounced for quiescent systems. They also note that, during quenched phases, galaxy sizes increase systematically with time, particularly among poorer-resolved low-mass systems, despite little growth in stellar mass over the same period. The secular growth of sizes of non-star forming galaxies (e.g. dwarfs or ellipticals) is an expected consequence of (spurious) energy equipartition between stellar and DM particles.

Our explanation is that 2-body scattering leads to a slow diffusion of stellar particles out of the dense central regions of galaxies. This is consistent with the simulations of Revaz & Jablonka (2018, $\mu=21.9$), in which quenched dwarf galaxies grow systematically in size with decreasing $z$, despite their passive evolution. Indeed, Revaz & Jablonka (2018) hypothesise that this result is due to 2-body scattering.

We note that hydrodynamical simulations involve gas particles that may also be subject to mass segregation if their masses differ from those of the DM, but in this case there is a compounding effect: collisions with DM particles also tend to heat gas particles as the kinetic energy associated with velocity perturbations thermalises (see Steinmetz & White, 1997, for details). Disentangling these effects (mass segregation and gravitational heating) will be challenging, and requires further study.

Finally, we note that assessing the impact of equipartition on galaxy sizes does not require time-consuming, high-resolution simulations of large volumes. Since the effect appears limited to haloes/galaxies of relatively low-particle number it can be gauged by comparing runs in relatively small-volumes that reach the target stellar mass resolution but vary $\mu$. 2-body scattering may also affect velocity dispersion and anisotropy profiles, angular momentum distributions and gas fractions. These issues will be addressed in a follow-up paper.

Acknowledgements

We thank Adrian Jenkins and Chris Power for useful conversations, and our referee, Alexander Knebe, for a useful report. ADL is supported by the Australian Research Council (project Nr. FT160100250). This work used the DiRAC@Durham facility managed by the Institute for Computational Cosmology on behalf of the STFC DiRAC HPC Facility (www.dirac.ac.uk). The equipment was funded by BEIS capital funding via STFC capital grants ST/K00042X/1, ST/P002293/1, ST/R002371/1 and ST/S002502/1, Durham University and STFC operations grant ST/R000832/1.