The kinematics of star clusters undergoing gas expulsion in Newtonian and Milgromian dynamics

TL;DR

This study compares star cluster kinematics during gas expulsion under Newtonian and Milgromian dynamics, revealing distinct velocity dispersion profiles and implications for galaxy formation theories.

Contribution

It provides the first detailed comparison of gas expulsion effects on star cluster kinematics in Newtonian and Milgromian (MOND) frameworks.

Findings

Gas expulsion explains NGC 2419's velocity dispersion in Newtonian dynamics.

MOND shows less steep central velocity dispersion slopes than Newtonian.

Low star formation efficiencies can produce ultra-diffuse objects.

Abstract

We study the kinematics of stars in clusters undergoing gas expulsion in standard Newtonian dynamics and also in Milgromian dynamics (MOND). Gas expulsion can explain the observed line-of-sight (LoS) velocity dispersion profile of NGC 2419 in Newtonian dynamics. For a given star formation efficiency (SFE), the shapes of the velocity dispersion profiles, which are normalised by the velocity dispersion at the projected half-mass radius, are almost indistinguishable for different SFE models in Newtonian dynamics. The velocity dispersion of a star cluster in the outer halo of a galaxy can indeed have a strong radial anisotropy in Newtonian dynamics after gas expulsion. MOND displays several different properties from Newtonian dynamics. In particular, the slope of the central velocity dispersion profile is less steep in MOND for the same SFE. Moreover, for a given SFE, more massive embedded cluster models result in more rapidly declining central velocity dispersion profiles for the final star clusters, while less massive embedded cluster models lead to flatter velocity dispersion profiles for the final products. The onset of the radial-orbit instability in post-gas-expulsion MOND models is discussed. SFEs as low as a few percent, typical of molecular clouds, lead to surviving ultra-diffuse objects. Gas expulsion alone is unlikely the physical mechanism for the observed velocity dispersion profile of NGC 2419 in MOND.

Tables

Table 1. Table 1: Parameters of ICs for the embedded (pre-gas expulsion) clusters with a Plummer density profile: the 1 st subscript 1 st 1_{\rm st} through 3 rd subscript 3 rd 3_{\rm rd} columns are the ID of the ICs, the total mass, M b subscript 𝑀 b M_{\rm b} , and the Plummer radius, r P subscript 𝑟 P r_{\rm P} . The 4 th subscript 4 th 4_{\rm th} through 6 th subscript 6 th 6_{\rm th} columns are the projected half-mass radius of the ICs ( R h subscript 𝑅 h R_{\rm h} ) and of the final star clusters after undergoing gas expulsion ( R h , f subscript 𝑅 h f R_{\rm h,f} ). The projected central density, Σ 0 subscript Σ 0 \Sigma_{\rm 0} for the ICs and Σ 0 , f subscript Σ 0 f \Sigma_{\rm 0,f} for the final star clusters, are shown in the 7 th subscript 7 th 7_{\rm th} through 9 th subscript 9 th 9_{\rm th} columns, respectively. The 10 th subscript 10 th 10_{\rm th} through 12 th subscript 12 th 12_{\rm th} columns represent the line-of-sight velocity dispersion at projected half-mass radii. The LoS direction is along x 𝑥 x -axis for the normalisation constants in the 7 th − 12 th subscript 7 th subscript 12 th 7_{\rm th}-12_{\rm th} columns. The fraction of the mass of bound stars in the final star clusters to the mass of the initial embedded clusters are listed, respectively, in 13 th subscript 13 th 13_{\rm th} and 14 th subscript 14 th 14_{\rm th} columns for different SFE.

ICs	$M_b$	$r_P$	$R_h$	$R_{h,f}$	$R_{h,f}$	${\Sigma }_0$	${\Sigma }_{0,f}$	${\Sigma }_{0,f}$	${\sigma }_0$	${\sigma }_{0,f}$	${\sigma }_{0,f}$	$f_b$	$f_b$
	$\left(M_{\odot }\right)$	$\left(\mathrm{pc}\right)$	$\left(\mathrm{pc}\right)$	$\left(\mathrm{pc}\right)$	$\left(\mathrm{pc}\right)$	$\left(M_{\odot }{\mathrm{pc}}^{-2}\right)$	$\left(M_{\odot }{\mathrm{pc}}^{-2}\right)$	$\left(M_{\odot }{\mathrm{pc}}^{-2}\right)$	$\left(\mathrm{km}{s}^{-1}\right)$	$\left(\mathrm{km}{s}^{-1}\right)$	$\left(\mathrm{km}{s}^{-1}\right)$
$ϵ$				0.5	0.4		0.5	0.4		0.5	0.4	0.5	0.4
N1	${10}^7$	5.0	5.0	15.0	28.9	$8.21\times {10}^4$	$6.88\times {10}^3$	$1.99\times {10}^3$	27.82	9.18	3.75	0.73	0.34
M1	${10}^7$	5.0	5.0	16.2	34.7	$7.97\times {10}^4$	$6.91\times {10}^3$	$1.96\times {10}^3$	28.08	11.29	7.69	0.96	0.87
M2	${10}^6$	5.0	5.0	10.0	13.3	$8.21\times {10}^3$	$1.18\times {10}^3$	$6.35\times {10}^2$	9.28	5.55	4.80	1.00	0.98
M3	${10}^5$	5.0	5.0	7.1	8.0	$8.67\times {10}^2$	$2.16\times {10}^2$	$1.25\times {10}^2$	3.66	2.77	2.58	0.99	0.98
M4	${10}^5$	10.0	10.0	12.6	13.7	$2.13\times {10}^2$	$6.34\times {10}^1$	$3.85\times {10}^1$	3.23	2.63	2.48	0.97	0.94

Table 2. Table 2: The observed values of the velocity dispersion profile of NGC 2419 (Ibata et al., 2011b ) .

Sample A
R (pc)	8.29	22.05	32.40	47.83	63.18	122.54
${\sigma }_{\mathrm{LoS}}\left(\mathrm{km}{s}^{-1}\right)$	${6.573}_{-0.575}^{+1.413}$	${3.816}_{-0.467}^{+1.032}$	${2.645}_{-0.378}^{+0.826}$	${2.662}_{-0.454}^{+0.731}$	${0.997}_{-0.549}^{+0.576}$	${0.756}_{-0.318}^{+0.436}$
${\sigma }_{\mathrm{LoS}}/{\sigma }_{R_h=19pc}$	${1.55}_{-0.14}^{+0.33}$	${0.90}_{-0.11}^{+0.24}$	${0.62}_{-0.09}^{+0.19}$	${0.63}_{-0.11}^{+0.17}$	${0.23}_{-0.13}^{+0.14}$	${0.18}_{-0.07}^{+0.10}$
Sample A+B
R (pc)	8.51	22.25	32.76	48.63	63.50	126.37
${\sigma }_{\mathrm{LoS}}\left(\mathrm{km}{s}^{-1}\right)$	${7.061}_{-0.698}^{+1.420}$	${3.778}_{-0.698}^{+1.002}$	${3.896}_{-0.504}^{+1.079}$	${3.109}_{-0.486}^{+0.779}$	${1.994}_{-0.651}^{+0.664}$	${1.295}_{-0.381}^{+0.683}$
${\sigma }_{\mathrm{LoS}}/{\sigma }_{R_h=19pc}$	${1.46}_{-0.14}^{+0.29}$	${0.78}_{-0.09}^{+0.21}$	${0.81}_{-0.10}^{+0.22}$	${0.64}_{-0.10}^{+0.16}$	${0.41}_{-0.13}^{+0.14}$	${0.27}_{-0.08}^{+0.14}$

Table 3. Table 3: Parameters of ICs for the embedded (pre-gas expulsion) clusters with a Plummer density profile to reproduce the surface density and LoS velocity profiles of NGC 2419 in Newtonian dynamics: the 1 st subscript 1 st 1_{\rm st} through 4 th subscript 4 th 4_{\rm th} columns are the ID of the ICs, the total mass, M b subscript 𝑀 b M_{\rm b} , the Plummer radius, r P subscript 𝑟 P r_{\rm P} and the half-mass radii, R h subscript 𝑅 h R_{\rm h} , of the ICs. The 5 th subscript 5 th 5_{\rm th} through 7 th subscript 7 th 7_{\rm th} columns are the star formation efficiency, ϵ italic-ϵ \epsilon , the half-mass radii of the final star clusters after undergoing gas expulsion ( R h , f subscript 𝑅 h f R_{\rm h,f} ) and the total mass of the final star clusters within the tidal radius, M f subscript 𝑀 f M_{\rm f} . The mass-to-light ratio (V-band) of the models are listed in the 8 th subscript 8 th 8_{\rm th} column.

ICs	$M_b$	$r_P$	$R_h$	$ϵ$	$R_{h,f}$	$M_f$	${\rm Y}$
	$\left(M_{\odot }\right)$	$\left(\mathrm{pc}\right)$	$\left(\mathrm{pc}\right)$		$\left(\mathrm{pc}\right)$	$\left({10}^5{M}_{\odot }\right)$
N2	$5\times {10}^6$	3.25	3.25	0.4	18.0	7.48	1.70
N3	$7\times {10}^6$	3.25	3.25	0.4	17.5	10.29	2.34
N4	$3\times {10}^6$	6.50	6.50	0.5	18.5	10.83	2.46

Equations

$${\rho }_{\mathrm{b}}(r)=\left(\frac{3M_{\mathrm{b}}}{4\pi r_{\mathrm{P}}^3}\right){\left(1+\frac{r^2}{r_{\mathrm{P}}^2}\right)}^{-5/2},$$

$$r_{\mathrm{vir}}\equiv \sqrt{G{M}_{\mathrm{b},\mathrm{f}}{a}_0}/g_{\mathrm{ext}},$$

$$I_{\mathrm{xx}}(\mathcal{R})=\frac{{\sum }_i{m}_i(y_i^2+z_i^2)/{\mathcal{R}}_i^2}{{\sum }_i{m}_i}.$$

$$I_{\mathrm{xy}}(\mathcal{R})=\frac{{\sum }_i-m_i{x}_i{y}_i/{\mathcal{R}}_i^2}{{\sum }_i{m}_i}.$$

$$a^{\prime }(\mathcal{R})=\sqrt{(I_{{\mathrm{y}}^{\prime }{\mathrm{y}}^{\prime }}+I_{{\mathrm{z}}^{\prime }{\mathrm{z}}^{\prime }}-I_{{\mathrm{x}}^{\prime }{\mathrm{x}}^{\prime }})/2},$$

$$b^{\prime }(\mathcal{R})=\sqrt{(I_{{\mathrm{x}}^{\prime }{\mathrm{x}}^{\prime }}+I_{{\mathrm{z}}^{\prime }{\mathrm{z}}^{\prime }}-I_{{\mathrm{y}}^{\prime }{\mathrm{y}}^{\prime }})/2},$$

$$c^{\prime }(\mathcal{R})=\sqrt{(I_{{\mathrm{x}}^{\prime }{\mathrm{x}}^{\prime }}+I_{{\mathrm{y}}^{\prime }{\mathrm{y}}^{\prime }}-I_{{\mathrm{z}}^{\prime }{\mathrm{z}}^{\prime }})/2}.$$

$$\Sigma /{\Sigma }_0(NGC2419)(\mathrm{R})=\frac{{\mathrm{I}}_{\mathrm{V}}(\mathrm{R})\cdot \mathrm{{\rm Y}}}{{\mathrm{I}}_{\mathrm{V},{\mathrm{R}}_{\mathrm{h}}=19pc}\cdot \mathrm{{\rm Y}}},$$

$${\sigma }_{\mathrm{LoS}}^2(R)=\frac{1}{N_{\mathrm{bound}}}\sum _{i=1}^{N_{\mathrm{bound}}}(v_{\mathrm{x},\mathrm{i}}-{\bar{v}}_{\mathrm{x}}{)}^2,$$

$$R_j=R_{\mathrm{h},\mathrm{f}}\tan [(j+0.5)0.5\pi /(n_R+1)],$$

$$\beta (r)\equiv 1-({\sigma }_{\theta }^2+{\sigma }_{\varphi }^2)/2{\sigma }_r^2,$$

$$t_{\mathrm{cr}}=R_{\mathrm{h},\mathrm{f}}/{\sigma }_{0,\mathrm{f}}.$$

Full text

1 Introduction

Star clusters form in collapsing cloud clumps within giant molecular clouds (GMCs, Tilley & Pudritz, 2004), and most stars are born in embedded star clusters (Lada & Lada, 2003; Kroupa, 2005, 2008; Marks & Kroupa, 2011). The timescale of evolution from a pure gaseous protocluster through the embedded cluster phase and to an exposed very young cluster is approximately $1~{}\,{\rm Myr}$. The gas is expelled from the embedded cluster through multiple mechanisms including stellar winds, type II supernova explosions, ionizing radiation, radiation pressure and outflows (Hopkins et al., 2013; Dib et al., 2013; Krumholz & Matzner, 2009; Murray et al., 2010; Hansen et al., 2012). The residual gas leaves with an outflow velocity of approximately $12\,{\rm km\,s}^{-1}$ for the observed Treasure Cheast cluster (Smith et al., 2005) and $25-30\,{\rm km\,s}^{-1}$ for the star clusters observed in the Antennae galaxies (Whitmore et al., 1999; Zhang et al., 2001). An outflow velocity of $10\,{\rm km\,s}^{-1}$ is commonly used in simulations (Kroupa et al., 2001; Banerjee & Kroupa, 2012, 2015, 2017; Kroupa et al., 2018). The star formation stops after gas expulsion. The gravitational potential well of the exposed star cluster is much shallower compared to that of the embedded cluster (e.g., Goodwin, 1997; Bastian & Goodwin, 2006; Goodwin & Bastian, 2006; Boily & Kroupa, 2003; Geyer & Burkert, 2001; Baumgardt & Kroupa, 2007; Marks et al., 2012; Pfalzner & Kaczmarek, 2013; Brinkmann et al., 2017; Banerjee & Kroupa, 2017), and the binding energy for the stars becomes smaller. Hence, the stars in the star cluster move onto elongated orbits.

NGC 2419, with ${\rm L_{\rm V}}\approx 4.4\times 10^{5}L_{\odot}$ (Baumgardt et al., 2009), is one of the brightest and most distant globular clusters (GCs) in the Milky Way (Bellazzini, 2007; Harris, 1996). The Galactocentric distance of NGC 2419 is $89.9~{}\,{\rm kpc}$. NGC 2419 is a metal-poor GC with a metallicity of $[Fe/H]\approx-2.32\pm 0.11$ dex and an alpha element abundance $[\alpha/Fe]=+0.2$ dex (Shetrone et al., 2001). The size of NGC 2419, with a projected half-mass radius of $R_{\rm h}\approx 19\,{\rm pc}$, is about five times larger than that of a normal GC at the same luminosity (Mackey & van den Bergh, 2005; Ibata et al., 2011b). Hence, both the background and the internal gravities of NGC 2419 are in the weak field regime, where the gravitational acceleration is much smaller than Milgrom’s gravitational constant, $a_{0}\approx 3.7{\rm pc/Myr^{2}}$. In the weak field regime, the gravitation is in the deep Milgromian dynamics (hereafter referred to as MOND, Milgrom, 1983) limit, providing a useful platform to test gravitational dynamics.

The velocity dispersion of stars in NGC 2419 has been measured by Baumgardt et al. (2009) and Ibata et al. (2011b). Baumgardt et al. (2009) found that a model with an isotropic velocity dispersion profile does not agree with the observed velocity dispersion profile in both Newtonian dynamics and Milgromian dynamics. Later, Sollima & Nipoti (2010) found that a radially anisotropic model is required to reproduce the velocity dispersion profile of NGC 2419 in Newtonian dynamics, while the radially anisotropic MOND models do not agree with the observed velocity dispersion profile of Baumgardt et al. (2009). Further simulations by Ibata et al. (2011b, a) have corroborated that the velocity dispersion of NGC 2419 is extremely radially anisotropic in both Newtonian and Milgromian dynamics. NGC 2419 has also been further extensively studied in both dynamics by other groups (Conroy et al., 2011; Brüns & Kroupa, 2011; Sanders, 2012b, a; Derakhshani & Haghi, 2014). In general, all the dynamical models are required to be strongly radially anisotropic in the two frameworks. Although an intermediate-mass black hole (IMBH) might also account for the very steep velocity dispersion profile in the central region of a star cluster, the surface density will have a shallow cusp in the center (Baumgardt et al., 2005), which is not the case in NGC 2419. The surface density profile of NGC 2419 is flat in the center according to several authors’ observations (Trager et al., 1995; Bellazzini, 2007; Ibata et al., 2011b).

To explain the velocity dispersion profile of NGC 2419, it is important to explore the possible mechanism for creating the required extremely radial anisotropy in both dynamics in the existing dynamical studies. A possible mechanism to generate such an anisotropy is partial relaxation (Lynden-Bell, 1967; Bertin & Trenti, 2003). A family of models undergoing partial relaxation has been investigated for a sample of GCs including NGC 2419 in Newtonian dynamics (Zocchi et al., 2012). A system undergoing partial relaxation is strongly radially anisotropic in its outer regions and thus may explain the kinematics of NGC 2419.

In this work, we shall propose another physically plausible origin of such a strong radial anisotropy. In the process of rapid gas expulsion, due to the weakening of the gravitational potential, the size of a star cluster expands violently and a large number of stars consequently move onto radial orbits (Brinkmann et al., 2017). The radial orbits generated by gas expulsion might lead to the radial anisotropy of the velocity dispersion. We consider the gas expulsion process in both Newtonian dynamics and MOND, since NGC 2419 is frequently used to test gravity in the two frameworks, as aforementioned.

In Milgromian dynamics, there is a phase transition from (quasi-)Newtonian dynamics to MOND in the process of gas expulsion, and the amount of phantom dark matter predicted by MOND increases with the weakening of gravitation.111Phantom dark matter is the effective but unreal matter generated from the MONDian gravitational potential by solving the standard Newtonian Poission’s equation (Milgrom, 1986; Wu et al., 2008). For a given star formation efficiency (SFE, hereafter denoted by $\epsilon$), after the gas expulsion, the gravitational potential well of a star cluster is deeper in MOND than that in Newtonian dynamics, i.e., the variation of the gravitational potential is smaller in MOND. It has been found that the size expansion of a star cluster undergoing gas expulsion is smaller in MOND than that in Newtonian dynamics with the same $\epsilon$ (Wu & Kroupa, 2018). Accordingly, the stars lose less potential energy in a MONDian system, and the orbits of these stars are less elongated. As a result, for a given $\epsilon$, the radial anisotropy of the star cluster is smaller in MOND. Moreover, a MONDian embedded cluster could leave a bound remnant with a much lower $\epsilon_{\rm crit}$ than a Newtonian system (Wu & Kroupa, 2018). An interesting question then arises as to whether or not the gas expulsion can result in an extremely radial anisotropy for the MOND models of NGC 2419.

In this work, the anisotropy profiles of star clusters undergoing gas expulsion is studied in both Newtonian dynamics and MOND. The paper is organised as follows: the $N$-body initial conditions (ICs) for the embedded clusters are constructed in the two dynamics, and a sudden gas expulsion process is applied to the embedded clusters by instantaneously reducing the masses of the $N$-body particles in Sec. 2. The models are evolved for $2\,{\rm Gyr}$ and the anisotropy profiles are compared in the two dynamics in Secs. 2.2 and 2.3. Since MOND allows a much lower $\epsilon_{\rm crit}$ to leave a bound object, we shall study the anisotropy profiles of the final bound objects with ultra-low $\epsilon$ (Sec. 3), which are impossible to leave bound cluster remnants in Newtonian dynamics (Brinkmann et al., 2017). We will summarise the results in Sec. 4.

2 Embedded clusters and gas expulsion

2.1 Models in Newtonian dynamics and in MOND

MOND has been tested for decades on the scale of galaxies, and the flat rotation curves for galaxies at large radii can be reproduced and in fact predicted from the observed density profile of luminous matter (e.g., Sanders & McGaugh, 2002; Famaey & McGaugh, 2012; Famaey & Binney, 2005; Famaey et al., 2007) by the empirical proposal of Milgrom (1983) and Bekenstein & Milgrom (1984) without dark matter particles. Moreover, on the scale of star clusters, a number of recent observations and simulations in the Milky Way suggest deviations from Newtonian dynamics. The line-of-sight (LoS) velocity dispersion profiles, $\sigma_{\rm LoS}(R)$, of these GCs are flat at large radii (Scarpa et al., 2007, 2011; Scarpa & Falomo, 2010; Lane et al., 2009; Lane et al., 2010; Hernandez & Jiménez, 2012; Hernandez et al., 2013; Durazo et al., 2017). In the case of MOND, the overall LoS velocity dispersion for an isolated and isotropic virialised system, $\sigma_{\rm LoS}^{2}=\frac{1}{3}\sigma^{2}=\frac{2}{9}(GM_{\rm b}a_{0})^{1/2}$, where the baryonic mass is $M_{\rm b}$ (Milgrom, 1994). The $\sigma_{\rm LoS}(R)$ profile for a virialised and isotropic MOND system within an external field has been fitted by an analytic function in Haghi et al. (2009); Kroupa et al. (2018) and Haghi et al. (2019, submitted), where $R$ is the 2-dimensional radius on the projected plane.

On large scales, MOND is confronted with several problems: the possible requirement of additional dark matter on the scale of clusters of galaxies (e.g., the lensing data from clusters, Natarajan & Zhao 2008; Angus et al. 2007, the Bullet Cluster 1E0657-56, Clowe et al. 2006, and the dark matter ring surrounding the rich galaxy cluster Cl0024+17, Jee et al. 2012) and the incorrect amplitude for the mass function of galaxy clusters (Angus & Diaferio, 2011; Angus et al., 2013, 2014). In addition, There is an inconsistency between MOND predictions and the observations of the third peak of the acoustic power spectrum of the CMB (Skordis et al., 2006; Zhao, 2008; Angus, 2009). MOND is a classical theory of gravitation and cannot readily be applied to cosmological problems without additional assumptions.

A remarkable success of MOND is that it precisely predicts the observed correlation between the baryonic matter and apparent dark matter in a galaxy. It has been pointed out that for example, the empirical mass-discrepancy-to-acceleration (MDA) correlation, cannot be exactly reproduced by the cold dark matter (CDM) hypothesis (Sanders, 1990; McGaugh, 2004; Famaey & McGaugh, 2012; Wu & Kroupa, 2015). Although there are several recent simulational results showing that a correlation between the baryonic and dark matter distribution in a galaxy can be found in the $\Lambda$CDM framework (Di Cintio & Lelli, 2016; Read et al., 2016; Keller & Wadsley, 2017; Navarro et al., 2017; Ludlow et al., 2017), the scatter of the relation from the $\Lambda$CDM simulations is generally much higher than that from the observations (see also Wu & Kroupa 2015). Therefore, it is important to continue studying MOND as a competitive alternative to the significantly challenged (Kroupa et al., 2010; Kroupa, 2012, 2015; Famaey & McGaugh, 2012) hypothesis of CDM.

2.1.1 $N$-body initial conditions and numerical tools

The initial parameters for the models of the embedded star clusters are the same as those in Wu & Kroupa (2018). The models are spherically symmetric and isotropic, with a Plummer density distribution (Plummer, 1911):

[TABLE]

where $M_{\rm b}$ and $r_{\rm P}$ are the total mass and the scale length of the Plummer model. Besides the same MONDian models presented in Wu & Kroupa (2018, models M1-M4 in Table 1), an equilibrium model (model N1) in Newtonian dynamics is constructed with the same density distribution as that of M1. We only construct one Newtonian model to compare with the MOND models, since the Poisson equation is linear and gravity is scale-free in Newtonian dynamics. The $N$-body initial conditions (ICs) are constructed using Lucy’s method (Lucy, 1974). The parameters of the Newtonian and MONDian models are summarised in Table 1. There are $10^{5}$ equal-mass particles in each model. For the MOND models we used the same interpolating function, $\mu(X)=\frac{X}{1+X}$, as that in Wu & Kroupa (2018), where $X=\frac{|{\bf g}|}{a_{0}}$ and ${\bf g}$ is the gravitational acceleration in MOND.

Since the number of particles is large, two-body relaxation can be ignored in the simulations. The evolution of models after the gas expulsion is calculated by using the collision-less particle-mesh $N$-body code NMODY (Londrillo & Nipoti, 2009). NMODY calculates gravitational accelerations and potentials by solving Poisson’s equation in both standard Newtonian and Milgromian dynamics (Nipoti et al., 2007; Nipoti et al., 2008; Nipoti et al., 2011; Wu & Kroupa, 2013; Wu et al., 2017). A grid resolution of $n_{r}\times n_{\theta}\times n_{\phi}=256\times 32\times 64$ is chosen for the simulations, where $n_{r},~{}n_{\theta}$ and $n_{\phi}$, respectively, are the number of cells in radial, polar and azimuthal dimensions. Here $n_{r}$ cells are used within a radial range of $0.15\,{\rm pc}$-$65.32\,{\rm kpc}$. The segmentation of grids is the same as that in Sec. 2.2 in Wu & Kroupa (2018). The global time steps are $dt=\frac{0.1}{\sqrt{\rm max|{\bf\nabla\cdot}{\bf g}|}}$, which ensures that there are 10 time steps for a circular orbit in the densest region of a system. The models are evolved for about $2\,{\rm Gyr}$ after gas expulsion.

We emphasise that a direct $N$-body algorithm does not exist for MOND, which is a mean-field description, because the generalised Poisson equation in MOND is non-linear. Thus only two-body relaxation free, i.e., collision-less, systems can at present be simulated in MOND. Applying star clusters to test MOND (Baumgardt et al., 2005) may thus not be straightforward. The linear standard Poisson equation of Newtonian gravitation, on the other hand, has an excellent direct $N$-body algorithm (Aarseth, 2003).

2.2 Gas expulsion in Newtonian dynamics and in MOND

Extensive numerical simulations have been performed showing that there is a critical value of the SFE, $\epsilon_{\rm crit}\approx 0.3$, for a single and equilibrium embedded cluster to leave a bound remnant after gas expulsion in Newtonian dynamics (Lada et al., 1984; Goodwin, 1997; Kroupa et al., 2001; Geyer & Burkert, 2001; Boily & Kroupa, 2003; Fellhauer & Kroupa, 2005; Baumgardt & Kroupa, 2007). In MOND, $\epsilon_{\rm crit}$ is much lower, and can be down to $2.5\%$ for an initially diffuse embedded cluster (Wu & Kroupa, 2018). In order to compare the kinematics of the final star cluster, i.e., the self-bound virialised object (see Sec. 2.3), in Newtonian dynamics and in MOND, we use $\epsilon=0.5$ and $0.4$ for the simulations in both dynamics. The above values of $\epsilon$ for an embedded cluster can leave a bound object after gas expulsion in both dynamics. Instantaneous gas expulsion is introduced by multiplying the mass of each particle by $\epsilon$ at the beginning of the simulations.

Since NGC 2419 is $89~{}\,{\rm kpc}$ away from the Galactic centre, the current tidal force from the Milky Way disc can be ignored in Newtonian dynamics. In MOND, the strong equivalence principle is violated due to the existence of the external field from the Galaxy. It is therefore necessary to consider both the internal and external gravitational accelerations when studying the dynamics of a system. At a Galactocentric distance of $89~{}\,{\rm kpc}$, the strength of the external field is $\approx 0.1a_{0}$ (Wu & Kroupa, 2015). Since the external field is much weaker than the internal field of the majority of NGC 2419222The internal gravitational acceleration of a system with a mass of $10^{6}M_{\odot}$ has been presented in Fig. 1 of Wu & Kroupa (2013). The internal acceleration approaches $0.1a_{0}$ at distances from the cluster centre larger than about $300~{}\,{\rm pc}$, which is more than $10R_{\rm h}$ for NGC 2419., the MOND models are simulated in isolation. The external field effect can be included by a truncation radius, i.e., the virial radius $r_{\rm vir}$,

[TABLE]

for the phantom dark matter halo (Wu & Kroupa, 2015, 2018). Here $g_{\rm ext}$ is the strength of the external gravitational acceleration. At the radius of $r_{\rm vir}$, the strength of the internal gravitational acceleration equals that of the external gravitational acceleration. Due to the external gravitational field, stars with sufficiently high energy are able to escape. The phantom dark matter halo is truncated at $r_{\rm vir}$, beyond which the external gravitational field from the Milky Way Galaxy dominates. $r_{\rm vir}$ is much smaller than the tidal radius in a MOND system (Wu & Kroupa, 2018). This indicates that the tidal field can be ignored in MOND. By the definition of $r_{\rm vir}$, the bound fraction of stars in a MOND final star cluster can be calculated from our simulations. The bound fraction of stars is the mass fraction of stars within $r_{\rm vir}$ in MOND, and within the tidal radius of the final star clusters in Newtonian dynamics.

2.2.1 Radial-orbit instability and shape of the final GCs

Since radial orbits are generated due to the shallower gravitational potential caused by gas expulsion, the star clusters might suffer from radial-orbit instability. For an initially spherical system, the radial-orbit instability alters the shape of the system into a triaxial morphology. Moreover, the velocity dispersion profile becomes more isotropic in the central region and more radially anisotropic in the outer region of the system (Barnes et al., 2005; Bellovary et al., 2008). The LoS velocity dispersion could differ along different projected angles of a star cluster, if the shape deviates significantly from a sphere. It is important to examine the radial-orbit instability of the star clusters undergoing gas expulsion.

The radial-orbit instability can be characterised by an anisotropy parameter, $\zeta\equiv\frac{2K_{\rm r}}{K_{\rm\theta}+K_{\rm\phi}}$, where $K_{\rm r},~{}K_{\rm\theta}$ and $K_{\rm\phi}$ are the kinematic energy on the radial, arzimuthal and polar directions in a spherical coordinate system (Polyachenko & Shukhman, 1981; Saha, 1992; Bertin et al., 1994; Trenti & Bertin, 2006). The radial-orbit instability appears in a spherical system when the value of $\zeta$ is larger than a critical value, $\zeta_{\rm crit}$. The values of $\zeta_{\rm crit}$ are in a range of around [1.6, 2.3] for spherical models with different distribution functions (e.g. May & Binney, 1986; Saha, 1991, 1992; Bertin & Stiavelli, 1989; Bertin et al., 1994). A spherical system is generically unstable when $\zeta\gtrsim 2.3$ in Newtonian dynamics. The radial-orbit instability of spherical Osipkov-Merritt MOND and Newtonian systems have been studied by Nipoti et al. (2011). The values of $\zeta_{\rm crit}\in[2.3,~{}2.6]$ in MOND systems and $\zeta_{\rm crit}\in[1.6,~{}1.8]$ in pure baryonic systems in Newtonian dynamics ($\xi_{s}$ in their table 1 is our $\zeta_{\rm crit}$). Nipoti et al. (2011) demonstrate that a pressure-supported MOND model is more radially unstable than an equivalent Newtonian system where a dark matter halo reproduces the potential which is equivalent to a MOND potential. But a MOND model is more stable than a pure Newtonian model without any dark matter. This is evident from the values of the critical instability parameter. We have examined the $\zeta$ values of all the final star clusters, and found that $\zeta\in[1.10,~{}1.72]$ for all the simulations in Table 1 (initial models $N1,~{}M1-M4$). This is consistent with the results of Nipoti et al. (2011), because none of the final star clusters have $\zeta$ exceeding the above range of values of $\zeta_{\rm crit}$. However, it is possible that the radial-orbit instability has been at work during the evolution of some of the simulated clusters. We address this question by calculating the shapes of all the final star clusters.

The shape of a final star cluster is described by the ratios of the principal axes, i.e., the intermediate-long axial ratio, $p$, and the short-long axial ratio, $q$. The $p$ and $q$ at a radius of $\mathscr{R}=\sqrt{x^{2}+(y/p)^{2}+(z/q)^{2}}$ can be calculated from the diagonalised normalised moments of the inertia tensor of the particle system, with an initial guess of $p=q=1.0$. The procedure follows (also see Wu et al., 2017):

• •

1. The ($I_{\rm xx}(\mathscr{R}),~{}I_{\rm yy}(\mathscr{R}),~{}I_{\rm zz}(\mathscr{R}),~{}I_{\rm xy}(\mathscr{R}),~{}I_{\rm xz}(\mathscr{R}),~{}I_{\rm yz}(\mathscr{R})$) of the normalised moments of inertia tensor, $I(\mathscr{R})$, are computed for a particle system within $\mathscr{R}$. Here

[TABLE]

$I_{\rm yy}(\mathscr{R})$ and $I_{\rm zz}(\mathscr{R})$ are defined in the same way. The component

[TABLE]

$I_{\rm xz}(\mathscr{R})$ and $I_{\rm yz}(\mathscr{R})$ are defined in the same way. $(x_{i},~{}y_{i},~{}z_{i})$ are the Cartesian coordinate of the $i_{th}$ particle within $\mathscr{R}$. The inertia tensor of each particle is normalised by its radius $\mathscr{R}_{i}$.

• •

2. The eigenvalues and eigenvectors of $I(\mathscr{R})$ are calculated by diagonalising the matrix of normalised inertia tensor. The eigenvalues are the new normalised $I_{\rm x^{\prime}x^{\prime}}(\mathscr{R}),~{}I_{\rm y^{\prime}y^{\prime}}(\mathscr{R})$ and $I_{\rm z^{\prime}z^{\prime}}(\mathscr{R})$, where $(x^{\prime},~{}y^{\prime},~{}z^{\prime})$ is the new reference frame which aligns with the principal axes of the particle system within $\mathscr{R}$. The normalised scale lengths of the principal axes of the particle system, $a^{\prime},~{}b^{\prime},~{}c^{\prime}$, can be obtained by

[TABLE]

The axial ratios of the particle system within $\mathscr{R}$ are $p=b^{\prime}/a^{\prime}$ and $q=c^{\prime}/a^{\prime}$.

• •

3. The whole final star cluster is rotated to align with the principal axes using the rotation matrix which diagonalises the normalised moments of inertia. In the new reference frame, the radius $\mathscr{R}$ is updated by using the new values of $p$ and $q$.

• •

4. We repeat the procedure from step 1 to step 3 until the axial ratios converge in the $j_{th}$ loop, in which $\Delta p_{j}=(p_{j}-p_{j-1})/p_{j-1}<10^{-4}$ and $\Delta q_{j}=(q_{j}-q_{j-1})/q_{j-1}<10^{-4}$.

The axial ratios of the final star clusters as functions of rescaled radial distance, $\mathscr{R}/R_{\rm h,f}$, are shown in Fig. 1. In Newtonian dynamics (left panel), for the final star cluster with $\epsilon=0.5$, $p$ is very close to 1.0 with $2\%$ deviation, while $q$ is about 0.9 in the innermost region ($\mathscr{R}<0.2R_{\rm h,f}$) and it grows quickly with the increase of radius. At a radius of $\approx 0.4R_{\rm h,f}$, $q\in[0.98,1.0]$. Such a final star cluster is thus almost perfectly spherical. For the final star cluster with $\epsilon=0.4$, the axial ratios are very similar to those of the star cluster with $\epsilon=0.5$. Therefore, the Newtonian models are almost spherical at $r>0.4R_{\rm h,f}$.

In MOND, the axial ratios of the final star clusters of models $M2-M4$ with $\epsilon$ of 0.4 and 0.5 behave very similarly to those in Newtonian dynamics. So do the axial ratios of the final star cluster of $M1$ with $\epsilon=0.5$. These MOND star clusters are almost spherical, especially in the outer regions with $\mathscr{R}>R_{\rm h,f}$. The shape of the final star cluster of $M1$ with $\epsilon=0.4$ behaves very differently to the Newtonian models and the other MOND models. In this star cluster, $p$ drops from near 1.0 to 0.79 within $\mathscr{R}\in[0.7,10.0]R_{\rm h,f}$, and $q$ declines from about 0.94 to 0.75 within $\mathscr{R}\in[0.5,10.0]R_{\rm h,f}$. This is a triaxial star cluster, which implies that the radial-orbit instability occurs in the $M1$ model with $\epsilon=0.4$ after gas expulsion. $M1$ is the most massive model, and the dynamics is expected to be the most similar to that of the Newtonian model. But this is not the case for the process of gas expulsion. Initially, models $N1$ and $M1$ are very similar, since the MOND effect is weakest in $M1$. When the gas is expelled from the star clusters with $\epsilon=0.4$, a large number of stars are moving along extremely radial orbits in both star clusters. These stars escape from the star cluster in Newtonian dynamics, but most of them remain bound in MOND within the truncation radius, $r_{\rm vir}$, where the strength of the internal and external field are comparable (see the bound fraction in Table 1). These stars account for the radial-orbit instability of the final star cluster of $M1$ in MOND. For the star cluster evolved from $M1$ with $\epsilon=0.5$, the potential does not change as much as that with $\epsilon=0.4$. For the $M2-M4$ models, since they are less massive than $M1$, the change of the MOND potential is smaller than that of Newtonian potential after gas expulsion. Fewer extremely radial orbits are generated in these MOND models. As a result, these final star clusters are nearly spherical. The radial-orbit instability in these models agree with the previous analysis of the instability parameter, $\zeta$.

We shall further study the surface density and kinematics within the final clusters in the following subsections. All the final star clusters are rotated to align with their principal axes.

2.2.2 Mass distribution after gas expulsion

At first, the major axis, i.e., $x$-axis, is chosen as the direction of the line-of-sight. The surface density distributions, $\Sigma(R)$, are presented in the upper panels of Fig. 2 for the Newtonian (upper left panel) models and the MONDian (upper right panel) models, where $R$ is the circular radial distance on the projected plane and the binning details are presented in Sec. 2.3. The $\Sigma(R)$ profiles are rescaled by the central surface densities, i.e., the value of $\Sigma$ at $R=0$. The values of the central surface densities are shown in Table 1, where $\Sigma_{0}$ is for the ICs and $\Sigma_{\rm 0,f}$ is for the final star clusters. After gas expulsion, the rescaled surface density profiles are very similar to the original rescaled Plummer profile in both Newtonian gravity and in MOND, which agrees well with the results in Boily & Kroupa (2003).

The V-band surface brightness profile of NGC 2419 measured by Ibata et al. (2011b, their Table 2) with errors is converted to $I_{\rm V}$ in units of $L_{\odot}\,{\rm pc}^{-2}$. A constant mass-to-light ratio, $\Upsilon$, is assumed here, as has been used for the modeling of NGC 2419 in Ibata et al. (2011b); McGaugh (2012), and the normalised surface density profile,

[TABLE]

is shown in the upper panels of Fig. 2 (red symbols). The observed surface density profile of NGC 2419 can be well reproduced by the rescaled Plummer’s profile. In general, the surface density profiles in both Newtonian dynamics and in MOND appear to be similar to that of the observed surface brightness profiles. The rescaling at the central regions of the star cluster makes it possible to directly compare the particle models and the observed surface brightness profile. One does not need to assume a mass-to-light ratio, which might introduce new uncertainties, for NGC 2419 when comparing the observations with models.

Since the final star clusters are not spherical in the central regions in both Newtonian and Milgromian dynamics, we also present the surface density profiles by choosing $y$-axis and $z$-axis as the LoS directions. The normalised surface density profiles on the $x-z$ plane and $x-y$ plane are shown in the upper panels of Figs. 3-4. We find that the normalised surface density profiles are very similar to the profiles shown in the upper panels of Fig. 2.

2.3 Kinematics

The middle panels of Fig. 2 display the LoS ($x$-axis direction) velocity dispersion, i.e.,

[TABLE]

where $N_{\rm bound}$ is the number of bound particles in the star cluster after gas expulsion and ${\bar{v}_{\rm x}}$ is the mean of the velocity component $v_{\rm x}$. Technically, different binning methods of the radius lead to slightly different profiles of $\sigma_{\rm LoS}(R)$. The radius $R$ is binned following the equation:

[TABLE]

where $n_{R}=20$ and $j=1,~{}2,~{}3,...,n_{R}$. The value of $\sigma_{\rm LoS}(R_{j+1})$ is the LoS velocity dispersion of particles within $R_{j}$ and $R_{j+1}$, and the related projected radius is the mean projected radius of particles within $R_{j}$ and $R_{j+1}$.

The $\sigma_{\rm LoS}(R)$ profiles are normalised by the effective velocity dispersion, i.e., velocity dispersion interpolated at the projected half-mass radius (denoted as $R_{\rm h}$ for the ICs and $R_{\rm h,f}$ for the 2 Gyr old star clusters after undergoing sudden gas expulsion initially at time T=0), $\sigma_{\rm 0}\equiv\sigma_{\rm LoS}|_{\rm r=R_{\rm h}}$ for the initial models and $\sigma_{\rm 0,f}\equiv\sigma_{\rm LoS}|_{\rm r=R_{\rm h,f}}$ for the final star clusters, respectively. The observed LoS velocity dispersion profile for NGC 2419 (Ibata et al., 2011b) is listed in Table 2 and plotted in the middle panels for comparison, where the normalisation parameter is $\sigma_{\rm Los}$ at the projected half-mass radius of NGC 2419 (i.e., 19pc). Moreover, the anisotropy profiles of the velocity dispersion,

[TABLE]

truncated at the virial radius, $r_{\rm vir}$ (Wu & Kroupa, 2018), are shown in the lower panels of Fig. 2, where $\sigma_{r}$, $\sigma_{\theta}$, and $\sigma_{\phi}$ are the components of the velocity dispersions in spherical coordinates $(r,~{}\theta,~{}\phi)$. The $\beta(r)$ profiles are isotropic for all the ICs.

In addition, the crossing time of the final star clusters is defined as

[TABLE]

The simulation time of $2~{}\,{\rm Gyr}$ is always longer than $300~{}t_{\rm cr}$ of the final star clusters in both dynamics listed in Table 1. The final star cluster models are fully re-virialised with a SFE between 0.4 and 0.5.

2.3.1 GCs in Newtonian dynamics

The Newtonian models predict larger values of $\sigma_{\rm LoS}/\sigma_{\rm 0,f}$ in the central regions for the final star clusters. The normalised $\sigma_{\rm LoS}(R)$ is larger for the N1 model with a smaller value of $\epsilon$. However, the normalised $\sigma_{\rm LoS}(R)$ profiles cannot be distinguished for N1 models with different $\epsilon$ in the outer region where $R>R_{\rm h,f}$. The $\sigma_{\rm LoS}(R)$ profiles for the Newtonian models agree well with those of the two observational samples (samples $A$ and $A+B$) of stars for NGC 2419 in Ibata et al. (2011b). The data are binned with an equal number of stars as the radius increases, and the radius of a bin is taken as the average radius of the observed stars in the binning range. Since the normalised $\sigma_{\rm LoS}(R)$ profiles for the final model star clusters are almost indistinguishable from the radius of the innermost bin of the observed data, any star cluster with $\epsilon\leq 0.5$ undergoing sudden gas expulsion should leave the same normalised $\sigma_{\rm LoS}(R)$ profile in Newtonian dynamics. For the case of a star cluster with $\epsilon<0.4$, a more rapidly declining $\sigma_{\rm LoS}(R)$ profile with increasing $R$ for radii smaller than $R_{\rm h,f}$ will be left, due to the fact that more radial orbits are generated through a larger fraction of mass loss. Observations show that the SFE is less than $0.3$ for the dense clumps in giant molecular clouds (Lada & Lada, 2003; Higuchi et al., 2009; Kainulainen et al., 2014; Megeath et al., 2016). Consequently, an embedded cluster model with $\epsilon>0.5$ disagrees with the observational constraints and we do not simulate models with $\epsilon>0.5$.

The $\beta(r)$ profiles (lower left panel) show that the final star clusters are radially anisotropic, especially in the outer regions of the star clusters where $R>R_{\rm h,f}$. The radial anisotropy increases when $\epsilon$ is smaller. The velocity dispersion for the final star cluster with $\epsilon=0.4$ is extremely radially anisotropic. For a model with $\epsilon=0.5$, $\beta(R_{\rm h,f})\approx 0.4$, while for $\epsilon=0.4$, $\beta(R_{\rm h,f})\approx 0.6$. These $\beta(r)$ profiles result in almost indistinguishable normalised $\sigma_{\rm LoS}(R)$ profiles, which agree well with the observed data of NGC 2419, especially for the star sample A+B in Ibata et al. (2011b). The process of gas expulsion predicts that the $\beta(r)$ profiles for sufficiently isolated GCs are strongly radially anisotropic. Note that the initial embedded star clusters are isotropic. The $\beta(r)$ profiles of the final star clusters might be less radially anisotropic if the embedded star clusters are strongly tangentially anisotropic and then undergo gas expulsion.

In addition, the normalised LoS velocity dispersion along the $y$-axis and $z$-axis for the same final star clusters are shown in the lower left panels of Figs. 3-4. The LoS velocity dispersion profiles are almost the same along the three principal axes’ directions. This further confirms again that gas expulsion can naturally explain the kinematics of NGC 2419, without introducing any assumptions on special directions as the line-of-sight.

Other physical processes, such as stellar evolutionary mass loss and two-body relaxation in a star cluster, will further drive the GCs to become more radially anisotropic. In the former case, the mass of a GC is further reduced and the gravitational potential becomes shallower. More stars move on elongated orbits in the evolving potential. In the latter case, the GC expands in the outer region and collapses in the core region. More radial orbits are expected in this process. However, a systematic study of the impact of these physical processes on the kinematics of GCs is still lacking. In addition, the tidal field from the Galaxy might allow the stars with high energy, i.e., stars on radial orbits, to escape. As a result, the velocity dispersion becomes less radially anisotropic. However, our models represent star clusters in the outer halo of a galaxy where the tidal field can be ignored. Therefore, a rapidly declining $\sigma_{\rm LoS}(R)$ profile is expected for GCs, especially for the GCs in the outer halo of a galaxy. Given that gas expulsion is physically realistic for all the GCs (see Marks et al., 2012; Zonoozi et al., 2011, 2014, 2017; Haghi et al., 2015), it is necessary to consider this process when one uses tidal heating models to explain the existence of observed flat $\sigma_{\rm LoS}(R)$ profiles at large radii for the Galactic GCs (Scarpa et al., 2007, 2011; Scarpa & Falomo, 2010; Lane et al., 2009; Lane et al., 2010; Hernandez & Jiménez, 2012; Hernandez et al., 2013; Durazo et al., 2017). However, the process of gas expulsion has been ignored in the above existing studies, therefore an initially strongly radially anisotropic velocity dispersion has not yet been taken into account. The details on the $\sigma_{\rm LoS}(R)$ and $\beta(r)$ profiles changed by stellar evolutionary mass loss, two-body relaxation and the tidal field in the inner region of the Galaxy under the consideration of gas expulsion in Newtonian dynamics, will be systematically studied in future projects by applying direct $N$-body simulations.

2.3.2 Newtonian models for NGC 2419

The surface density and velocity profiles shown in Fig. 2 are normalised by constants. Therefore, it is possible to reproduce the surface density profile and the LoS velocity profile of NGC 2419 in physical units by choosing specific normalisation parameters for the initial models.

The distance modulus of NGC 2419, $(m-M)_{\rm 0}=19.71\pm 0.08$ (Di Criscienzo et al., 2011). The extinction is $E(B-V)=0.08\pm 0.01{\rm mag}$ (Di Criscienzo et al., 2011). The Solar absolute V-band magnitude is $M_{\rm V,\odot}=4.83$. The observationally determined V-band mass-to-light ratio of NGC 2419, $\Upsilon=M/L_{\rm V}$, lies in a wide range of $[1.2,~{}2.55]$ (Bellazzini et al., 2012; Ibata et al., 2011b; Baumgardt et al., 2009). The $\Upsilon$ obtained from dynamical models are $1.9\pm 0.16$ (Ibata et al., 2011b) and $2.05\pm 0.5$ (Baumgardt et al., 2009), respectively. On the other hand, the $\Upsilon$ derived from stellar evolution models is $1.2-1.7$ and the best fitting value is $1.5\pm 0.1$ with $1\sigma$ error (Bellazzini et al., 2012).

To reproduce the mass, size and stellar velocities of NGC 2419, the initial mass and Plummer radius of the embedded cluster models need to be rescaled (the method is described by Madejsky & Bien, 1993). We present three sets of embedded cluster models which are possible to reproduce the mass distribution and kinematics of NGC 2419 in Table 3 and in Fig. 5. We find that the surface density and LoS velocity dispersion profiles of the final star clusters evolving from models $N3$ and $N4$ agree the best with those of NGC 2419. Model $N2$ predicts final star clusters with a lower central velocity dispersion profile, which does not agree with the observed $\sigma_{\rm LoS}(R)$ data within the $1\sigma$ error range. Moreover, the overall masses of the three models are within $7.48\times 10^{5}$-$1.08\times 10^{6}M_{\odot}$, and they agree with the dynamical mass suggested by Ibata et al. (2011b) and Baumgardt et al. (2009). The mass-to-light ratios, $\Upsilon$, of the models are within $1.70$-$2.46$, which are consistent with the observational V-band $M/L_{\rm V}$ (Bellazzini et al., 2012; Ibata et al., 2011b; Baumgardt et al., 2009). Consequently, it is possible to reproduce the mass distribution and kinematics of NGC 2419 by the mechanism of gas expulsion in Newtonian dynamics.

2.3.3 GCs in MOND

In the case of MOND, the normalised $\sigma_{\rm LoS}(R)$ profiles along the $x$-axis (middle right panel of Fig. 2) for the final star clusters decrease more slowly compared to the Newtonian models. The normalised $\sigma_{\rm LoS}(R)$ profile of the most compact model, i.e., the final star cluster evolved from the embedded cluster model M1, declines most steeply with radius. For this embedded cluster, when $\epsilon$ is larger (i.e., $\epsilon=0.5$), $\sigma_{\rm LoS}/\sigma_{\rm 0,f}$ is smaller in the centre and is slightly larger at large radii compared to that in Newtonian dynamics. The $\beta(r)$ profiles for the final star cluster of $M1$ with $\epsilon=0.4$ and $0.5$ are the most radially anisotropic. We find that the $\beta(r)$ profile for the final star cluster of model M1 is more radial than that of model $N1$ for the same SFE. Is this in contrast to the fact that the $\sigma_{\rm LoS}(R)/\sigma_{\rm 0,f}$ profile of the final star cluster of M1 does not decline with radius as fast as that of model $N1$? The $N1$ and $M1$ models are the most massive models in Newtonian dynamics and in MOND. The $\sigma_{\rm LoS}(R)/\sigma_{\rm 0,f}$ and $\beta(r)$ profile of the final M1 star cluster is therefore expected to be the most similar to that of the N1 model. As aforementioned, the more radially anisotropic velocity dispersion in MOND for the final M1 star cluster comes from a larger fraction of bound stars (see Table 1), including some extremely radial orbits which are escaping stars in N1 but are bound in M1. The very radial orbits only exist in the very compact and massive MOND models such as in M1. In this model, the potential is quasi-Newtonian at the beginning and evolves to an intermediate-MOND state after gas expulsion. More extremely radial orbits can be generated and also bound to the system. However, since the gravitational potential is deeper for the final star cluster of M1, all of the components of velocity dispersion ($\sigma_{\rm r},~{}\sigma_{\rm\theta}$ and $\sigma_{\rm\phi}$) are larger than those in Newtonian dynamics. This explains why the stronger radial anisotropy for the final star cluster of $M1$ cannot reproduce a rapidly declining $\sigma_{\rm LoS}(R)/\sigma_{\rm 0,f}$ profile like that of model $N1$.

When the embedded clusters become less compact, the slopes of the normalised $\sigma_{\rm LoS}(R)$ profiles decline with $R$ more gently. The values of $\sigma_{\rm LoS}(R)/\sigma_{\rm 0,f}$ for all the final star clusters are smaller in the centres and larger at large radii, compared to NGC 2419. The innermost bin of the observed normalised $\sigma_{\rm LoS}(R)$ is above the profile predicted by models and the last two bins of the observed data are below all the predicted $\sigma_{\rm LoS}(R)/\sigma_{\rm 0,f}$ curves in MOND. If MOND as a mean-field theory is applicable in star clusters, this indicates that gas expulsion in MOND is not the sole mechanism for the observed kinematics of NGC 2419. That is why the $\sigma_{\rm LoS}(R)/\sigma_{\rm 0,f}$ profiles of the models undergoing gas expulsion appear very different from those in Ibata et al. (2011b), where the models are constructed initially in equilibrium in the framework of MOND.

The LoS velocity dispersion along other principal axes, i.e., $y$- and $z$-axis, are illustrated in the lower right panels of Fig. 3 and Fig. 4, respectively. For the final star clusters evolved from $M1$, the central $\sigma_{\rm LoS}/\sigma_{\rm 0}$ profiles along the $y$- and $z$-axis are smaller than those along the $x$-axis for both values of $\epsilon$, especially for the $\epsilon=0.4$ model along the $z$-axis. Furthermore, in the outer region where $R=10~{}R_{\rm h,f}$, the values of $\sigma_{\rm LoS}/\sigma_{\rm 0}$ are larger along the $y$- and $z$-axis than those along the $x$-axis. In short, the normalised LoS velocity dispersion profiles are flatter along the intermediate and the short axes. This is related to the triaxial shape of the final star cluster. For other models, the $\sigma_{\rm LoS}/\sigma_{\rm 0}$ profiles along the $y$- and $z$-axis are very similar to those along the $x$-axis.

Owing to the deeper gravitational potential in MOND, for a given $\epsilon$, a final star cluster does not expand as much as in Newtonian dynamics. The stellar orbits in the final star clusters are less elongated. Therefore, the anisotropy profiles are less radial in MOND compared to Newtonian dynamics. Moreover, thanks to the fact that the effective mass of phantom dark matter is larger, relative to the baryonic mass, for the more diffuse systems, the size expansion is smaller for these models (see Wu & Kroupa 2015 and Wu & Kroupa 2018). Fewer radial orbits are generated in the more diffuse systems. This results in a more slowly declining $\sigma_{\rm LoS}(R)$ profile and a less radially anisotropic $\beta(r)$ profile for the more diffuse systems. The $\beta(r)$ profiles in the lower right panels of Fig. 2 confirm this.

On the other hand, the proper motion for the Virgo stellar stream (VSS) indicates that NGC 2419 may be associated with the orbit of the VSS, which is a highly eccentric Galactic orbit with a pericentre of $\approx 11~{}\,{\rm kpc}$ (Casetti-Dinescu et al., 2009). The external field effect from the Milky Way might account for the formation of an extreme radial anisotropy for NGC 2419 in MOND. That is, when a star cluster is moving on an eccentric Galactic orbit inwards and outwards, the strength of the external field from the Galaxy increases and decreases, respectively. As a result, the star cluster expands and contracts on such an orbit. Radial orbits in the star cluster are generated in this process, which might give rise to a radial anisotropy. The external field for outgoing orbits has been studied in Wu & Kroupa (2013). More sophisticated models including the effects of space-varying external field and gas expulsion will be studied for NGC 2419 in the future.

Moreover, since the $\sigma_{\rm LoS}(R)$ profiles for the final star clusters are not unique in MOND, it is possible to reproduce the observed flat $\sigma_{\rm LoS}(R)$ profiles at large radii for distant GCs by using initially diffuse embedded cluster models and small values of $\epsilon$ undergoing gas expulsion. Such diffuse initial conditions would result from star formation throughout a giant molecular cloud, perhaps comparable to the Cyg OB2 association (Knödlseder, 2000). Distant GCs with moderate or mild radial anisotropy are thus allowed to exist in MOND. In contrast, all the GCs are expected to be strongly radially anisotropic after gas expulsion in Newtonian dynamics, unless another physical mechanism is taken into account.

To summarise, gas expulsion alone does not account for the observed velocity dispersion profile for NGC 2419 in MOND. Other physical mechanisms, such as, the time-varying external field effect, reproduce the $\sigma_{\rm LoS}(R)$ profile for NGC 2419 (Wu & Kroupa, 2013). However, star clusters surviving gas expulsion in MOND can have a slowly declining $\sigma_{\rm LoS}(R)$ profile with a mild or moderate radial anisotropy. These models can explain the observed flat $\sigma_{\rm LoS}(R)$ profiles at large radii for the distant Galactic GCs, while flat $\sigma_{\rm LoS}(R)$ profiles cannot be reproduced by gas expulsion-only models in Newtonian dynamics. One has to invoke other physical mechanisms to account for these profiles within the framework of Newtonian gravity.

3 MOND models with ultra-low SFE

In Newtonian dynamics, the anisotropy profile of a final star cluster is more radially anisotropic when the SFE is lower for the same embedded cluster model. The question whether a smaller $\epsilon$ would make the $\sigma_{\rm LoS}(R)$ profile of a MOND model more compatible with the data from NGC 2419 then arises. Since models with ultra-low SFE can leave bound remnants after sudden gas expulsion in MOND (Wu & Kroupa, 2018), we shall further study the $\sigma_{\rm LoS}(R)$ and $\beta(r)$ profiles for such systems. Such low SFEs are typical for whole molecular clouds which spawn an OB association sufficiently massive to remain bound in MOND. The values of $\epsilon$ are chosen to be $30\%,~{}20\%$ and $10\%$ in the MOND models. In Newtonian dynamics, an embedded cluster with such values of SFE cannot survive sudden gas expulsion. For comparison, the embedded cluster models with $\epsilon$ of $0.5$ and $0.4$ are plotted in Fig. 6 as well. We now make a pure model analyses by tuning $\epsilon$.

The normalised $\sigma_{\rm LoS}(R)$ profiles are shown in Fig. 6. In model M$1$, $\sigma_{\rm LoS}(R)$ of the star clusters surviving with $\epsilon$ of $50\%$ and $40\%$ have larger normalised values in the centres, i.e., $\sigma_{\rm LoS}/\sigma_{\rm 0,f}\approx 1.4$ at $0.1~{}R_{\rm h,f}$, compared to the initial models, i.e., $\sigma_{\rm LoS}/\sigma_{0}\approx 1.2$. The velocity dispersion profiles of model M$1$ with these SFEs decline slowly within $0.5~{}R_{\rm h,f}$ and then drop rapidly with radius within $2~{}R_{\rm h,f}$. The normalised $\sigma_{\rm LoS}(R)$ profiles for these models decrease more slowly than those of the ICs again at radii $R>2R_{\rm h,f}$. Consequently, a compact embedded cluster with a large value of $\epsilon$ may result in a bound star cluster with a rapidly declining $\sigma_{\rm LoS}(R)$ profile. The velocity dispersion profiles are almost flat at all radii for $\epsilon$ of $30\%$ and $20\%$, and there is an increasing normalised $\sigma_{\rm LoS}(R)$ profile with increasing $R$ for model M$1$ when $\epsilon$ decreases down to $10\%$. The reason is that the bound mass for this surviving star cluster is small ($\approx 6\times 10^{4}M_{\odot}$, see Wu & Kroupa 2018), while the size of the star cluster is as large as $R_{\rm h,f}\approx 200\,{\rm pc}$. The final “star cluster”, formed from a giant molecular cloud complex with an overall SFE typical of molecular clouds ($\epsilon\in[0.2\%,~{}20\%]$ and with a luminous-weighted average of $8\%$ for a large sample of GMCs in the Galaxy, see Murray 2011), is an ultra-diffuse system, and accordingly, the relaxation time scale is longer than $2~{}\,{\rm Gyr}$. The growing $\sigma_{\rm LoS}(R)$ profile with $R$ indicates a not yet fully virialised surviving star cluster. This agrees with the evolution of Lagrangian radii for this star cluster (Wu & Kroupa, 2018).

The values of the central normalised velocity dispersion, $\sigma_{\rm LoS}(R\approx 0.1R_{\rm h,f})/\sigma_{\rm 0,f}$, are smaller when the models are less massive or more diffuse. Noteworthily, the normalised $\sigma_{\rm LoS}(R)$ profiles of models M$2$-M$4$ almost overlap those of the ICs. This indicates that although the radii of the clusters and the values of $\sigma_{\rm 0,f}$ evolve after gas expulsion, the radial distribution of velocity dispersion of these models does not change. The shapes of $\sigma_{\rm LoS}(R)$ of models M$2$-M$4$ with $\epsilon=0.1$ are flat, which is similar to model M$1$. In addition, when the model is more compact, the overall $\sigma_{\rm LoS}(R)$ profile is flat for a larger value of $\epsilon$. For instance, model M1 with $\epsilon\leq 0.3$ leaves a flat $\sigma_{\rm LoS}(R)$ profile, while in models M$3$ and M$4$, the $\sigma_{\rm LoS}(R)$ profiles are flat only when $\epsilon$ is reduced down to $0.1$. This is also because, i) the size expansion is the largest for model M$1$ and the smallest for models M$3$-M$4$ with the same value of $\epsilon$ and ii) the bound fraction is the smallest for model M$1$ and the largest for models M$3$-M$4$ for the same SFE (Wu & Kroupa, 2018). The virialisation timescales for star cluster models with an ultra-low SFE are longer than the simulation time of 2 Gyr. Hence the flat $\sigma_{\rm LoS}$ profiles might be related to the not yet virialised status of the final star cluster. Accordingly, the smaller SFE does not form a star cluster with a more rapidly declining $\sigma_{\rm LoS}(R)$ profile for the same embedded cluster model in MOND. In contrast, when the value of SFE decreases, the $\sigma_{\rm LoS}(R)$ profile declines more slowly with radius or even becomes flat after being evolved for a few Gyr. The observed $\sigma_{\rm LoS}(R)$ profile of NGC 2419 cannot be reproduced by reducing the value of SFE in MOND. Although gas expulsion is not the sole mechanism for the observed kinematics in NGC 2419 in MOND, the slower decrease and the flat $\sigma_{\rm LoS}(R)$ profiles predicted by MOND are consistent with the existence of the observed flat $\sigma_{\rm LoS}(R)$ profiles for other GCs in the Galaxy (Scarpa et al., 2007, 2011; Scarpa & Falomo, 2010; Lane et al., 2009; Lane et al., 2010; Hernandez & Jiménez, 2012; Hernandez et al., 2013; Durazo et al., 2017). Although these GCs are not outer halo GCs and the external field from the Galaxy at their locations is about $1a_{0}$ (Lane et al., 2010), the slowly declining and the flat $\sigma_{\rm LoS}$ profiles after gas expulsion for low SFEs are expected for them. Due to the fact that the external field enables the high energy stars to escape from the GCs, the strongly radially anisotropic orbits are lost. We will systematically study the tidal field and external field effects for GCs after undergoing gas expulsion in a future project.

In general, the $\beta(r)$ profiles for all the surviving star clusters become radially anisotropic at the end of the simulations (Fig. 7). The surviving star clusters are the most radially anisotropic for the massive and compact model, namely, model M$1$, with all values of SFE, while they are the least radially anisotropic for the least massive and diffuse model, model M$4$.

For model M$1$, the $\beta(r)$ profiles are very similar for all the surviving star clusters with SFE between 0.1 and 0.5. The $\beta(r)$ profiles, starting from near [math] in the central region, increase steeply to around $0.5$-$0.7$ at a radius of $2R_{\rm h,f}$, and then the profiles rise slowly out to large radii. The values of $\beta$ at $R=R_{\rm h,f}$ for model M2 are between $0.4$ and $0.6$, for model M$3$ they are in the range of $[0.2,~{}0.5]$ and decrease to $0.1-0.4$ for model M$4$. Given that model M$1$ is dominated by quasi-Newtonian gravity, while model M$4$ is dominated by deep MOND, the radial anisotropy is weaker for models dominated by MOND gravity. In a deep MOND system such as model M$4$, the internal gravity is much stronger than a pure Newtonian system with exactly the same density distribution. Therefore, the particles are more tightly bound to a deep MOND system. Since a MOND system does not expand as much as a Newtonian system when the two systems lose the same fraction of mass in gas, less radial orbits are generated in MOND after the gas expulsion.

4 Conclusion and discussions

In this work, we presented the first analysis of the line-of-sight velocity dispersion profiles and the anisotropy profiles of star clusters surviving gas expulsion in Newtonian and in Milgromian dynamics. A series of simulations with the SFE ranging from $10\%$-$50\%$ are performed for embedded clusters with initial masses in the range of $[10^{5},10^{7}]M_{\odot}$. We summarise and discuss our main results here.

The LoS velocity dispersion of the surviving star clusters are studied in both Newtonian dynamics and MOND, with $\epsilon=0.5$ and $0.4$. In Newtonian dynamics, the $\sigma_{\rm LoS}(R)/\sigma_{\rm 0,f}$ profiles for the final star clusters are almost indistinguishable in the region where $R>R_{\rm h,f}$. In the inner region where $R<R_{\rm h,f}$, the models with smaller values of $\epsilon$ have larger values of $\sigma_{\rm LoS}/\sigma_{\rm 0,f}$. On the other hand, the value of the observed $\sigma_{\rm LoS}(R)$ profile for NGC2419 reaches $\approx 7\,{\rm km\,s}^{-1}$ in the centre and drops to $1-2\,{\rm km\,s}^{-1}$ at large radii where $R>3\,{\rm R_{e}}$ ($\,{\rm R_{e}}$ is the effective radius of the star cluster, Baumgardt et al., 2009; Ibata et al., 2011b, a). The existing dynamical models usually require a highly radially anisotropic profile for this star cluster (Ibata et al., 2011b, a). Here we have found that gas expulsion naturally provides an explanation for the observed $\sigma_{\rm LoS}(R)$ and the model-required $\beta(r)$ profiles of NGC 2419 in Newtonian dynamics. Furthermore, all the isolated Newtonian GCs should have a strongly radially anisotropic distribution of velocity dispersion at least in the outer regions of the star clusters ($R>R_{\rm h}$), if gas expulsion with a low SFE is physically realistic for these GCs.

In contrast, GCs are allowed to have mildly or moderately radially anisotropic profiles in MOND undergoing gas expulsion with the same values of the SFE. In general, the central $\sigma_{\rm LoS}(R)/\sigma_{\rm 0,f}$ values are lower and the outer $\sigma_{\rm LoS}(R)/\sigma_{\rm 0,f}$ values are higher for the same model in MOND as in Newtonian dynamics. Radial-orbit instability occurs in the most massive MOND model with a SFE of 0.4, while the corresponding Newtonian model with the same model parameters is stable. The radial-orbit instability in the MOND post-gas-expulsion model comes from the bound high energy stars in the Milgromian gravitational potential. Such high energy stars escape from the final star cluster in Newtonian gravitation. Moreover, the $\sigma_{\rm LoS}(R)/\sigma_{\rm 0,f}$ profiles decline more slowly for the MOND models. The central $\sigma_{\rm LoS}/\sigma_{\rm 0,f}$ is smaller and the values of $\sigma_{\rm LoS}/\sigma_{\rm 0,f}$ in the outer regions are larger for the more diffuse models and also for models with smaller $\epsilon$. This stems from the gravitational potential well being deeper for the star clusters in MOND with the same density distribution as the Newtonian models, and the more diffuse models and the models with lower $\epsilon$ are in a deeper MOND regime.

In addition, star clusters can survive much lower $\epsilon_{\rm crit}$ in MOND (Wu & Kroupa, 2018). The $\sigma_{\rm LoS}(R)$ and $\beta(r)$ profiles are studied for the star clusters undergoing gas expulsion with $\epsilon=30\%$, $20\%$ and $10\%$. For the more compact embedded cluster, i.e., model M1, the $\sigma_{\rm LoS}(R)$ profile is almost flat for such a low SFE. When $\epsilon$ is down to $10\%$, the $\sigma_{\rm LoS}(R)$ profile is growing with radius. The possible reason is that the surviving star cluster is not fully virialised within the simulation time of $2~{}\,{\rm Gyr}$. When the model for the embedded cluster becomes more diffuse, the critical value of $\epsilon$ to generate a rapidly falling $\sigma_{\rm LoS}(R)$ profile decreases. For instance, the $\sigma_{\rm LoS}(R)$ profiles are falling fast with incrasing radius for $\epsilon\approx 0.2$ in models M3 and M4. In these models, only when $\epsilon$ is 0.1, the $\sigma_{\rm LoS}(R)$ profiles are flat. To summarise, the observed $\sigma_{\rm LoS}(R)$ profile of NGC 2419 cannot be reproduced by reducing the SFE undergoing gas expulsion in MOND.

The anisotropy profiles, $\beta(r)$, for these MOND models show that the more compact the embedded cluster is, the more radially anisotropic the final star cluster is. For the same embedded cluster with a compact density distribution (i.e., model M1), the $\beta(r)$ profiles for the final star clusters for different $\epsilon$ are very similar to each other. When the embedded cluster model becomes more diffuse, the final star clusters are more radially anisotropic for smaller values of $\epsilon$. Mildly ($\beta(r)\approx 0.2$ at large radii) and moderately ($\beta(r)<0.4$ at large radii) radial anisotropy profiles can be generated through gas expulsion in MOND, if the embedded cluster is a diffuse system and the gravity is dominated by deep MOND. This is because the particles in a MOND system are more tightly bound, and the gravitational potential does not change as much as that in quasi-Newtonian gravity after gas expulsion. There are fewer radial orbits generated by the sudden removal of gas. Therefore, MOND predicts a larger variety of velocity dispersion profiles and anisotropy profiles through the process of gas expulsion than Newtonian gravitation, which agrees with the variety of LoS velocity profiles for GCs observed in the Galactic halo.

We note in particular that while gas-expulsion-only models fail to account for the observed $\sigma_{\rm LoS}(R)$ profile of NGC 2419 (Fig. 2) in MOND, if this cluster is on a very eccentric orbit then the time-varying external field effect will lead to a reproduction of the observed profile (Wu & Kroupa, 2013). At the same time, the gas expulsion process in MOND models can also account for flat outer $\sigma_{\rm LoS}(R)$ profiles, while the sole gas expulsion Newtonian models cannot. It has been suggested that GCs may reside in their own dark matter halos which would possibly account for the observed flat $\sigma_{\rm LoS}(R)$ profiles. Moreover, the possibility of existence of a dark matter halo around NGC 2419 is discussed in Ibata et al. (2013). The dynamical models which reproduce the surface mass density and kinematics of NGC 2419 cannot fully rule out a dark matter halo around this star cluster, but the dark matter particle hypothesis is challenged by other data (Kroupa, 2012, 2015; Famaey & McGaugh, 2012). It remains unknown whether this possibility is consistent with models of star cluster formation (Moore, 1996). Another physical mechanism is required in Newtonian dynamics to account for the observed flat $\sigma_{\rm LoS}(R)$ profiles.

5 Acknowledgments

The authors thank the anonymous referee for the very valuable and helpful comments. The authors thank Luca Ciotti, Pasquale Londrillo and Carlo Nipoti for sharing their NMODY code. The authors acknowledge Ortwin Gerhard for sharing the code for generating the self-consistent N-body ICs using Lucy’s method. XW thanks for support through grants from the Natural Science Foundation of China (Numbers NSFC-11503025 and NSFC-11421303), a Anhui Natural Science Foundation grant (Number 1708085MA20), and “the Fundamental Research Funds for the Central Universities”. XW thanks for support from the “Hundred Talents Project of Anhui Province”.