Open Cluster Dynamics under the Influence of Outflow-Ambient Interactions

The interaction between gas and stars can generate either friction or anti-friction, depending on the strength of the stellar outflows. For stars that have no outflows, the analytic approximations for DF are adopted from Ostriker (1999):

where $v_{*}$ is the stellar velocity, $\rho$ is the ambient gas density, $\mathcal{M}\equiv v_{*}/c_{s}$ is the Mach number in the gas, and $\Lambda\equiv b_{\mathrm{max}}/b_{\mathrm{min}}$ is the Coulomb factor. In this work, we assume that the $a_{\rm DF}$ direction is always aligned with the relative motion direction, because the timescales of hydrodynamic phenomena $\tau_{\rm hyd}$ are typically much shorter than the stellar orbital dynamic timescales $\tau_{\rm dyn}$,

in which $l_{\rm hyd}$ and $l_{\rm dyn}$ are the typical spatial scale of the hydrodynamic process and the stellar dynamics, respectively. The overdense structures causing DF is typically a few times the stellar size in absence of outflows (see e.g., Ostriker 1999), and the typical orbital scales of stars in OCs are no less than about $10^{-1}~{}\mathrm{pc}$ (see e.g., Nilakshi et al. 2002) (unless in close encounter events, which are typically rapid and dwarf DF effects by the ambient gas). Therefore, it is almost always eligible to assume that the star-gas interaction structures settle at the steady or quasi-steady states, and the direction of DF stays along the stellar motion relative to the gas. We use the Bondi radius of each star, $b_{\mathrm{min}}\equiv 2GM/c_{s}^{2}$, for the minimum impact parameter, and we approximate $b_{\mathrm{max}}$ with the size upper limit of a giant molecular cloud ($100~{}\mathrm{pc}$ in this work; e.g., Solomon et al. 1987; Miville-Deschênes et al. 2017; Sun et al. 2018).

Whenever a star has outflows, the gravitational feedback coming from the interactions between such outflows and the ambient gas will likely accelerate the star. We adopt the analytic approximation described in Gruzinov et al. (2020) and Li et al. (2020),

where $R_{0}\equiv[\dot{m}_{\rm w}v_{\rm w}/(4\pi\rho v_{*}^{2})]^{1/2}$ is the standoff distance (where the total pressure of the incoming medium equals that of the outflow), $\dot{m}_{\rm w}$ is the outflow wind mass-loss rate, $v_{\rm w}$ is the wind radial velocity and $u\equiv v_{*}/v_{\rm w}$. $R_{\mathrm{s}}=R_{0}[3(1-\theta\cot\theta)/\sin^{2}\theta]^{1/2}$ is the contact discontinuity location at the polar angle $\theta$ (where the gas over-densities are located) in spherical coordinates for which the star is at the origin, and the axis $\theta=0$ is in the opposite direction of $v_{*}$. In general, for typical OCs, $v_{*}$ spans $0.5-5~{}\mathrm{km}~{}\mathrm{s}^{-1}$ (e.g., Tarricq et al., 2021; Soubiran et al., 2018; Mathieu, 1984), which is significantly lower than the $v_{\rm w}$ of almost all stellar outflows. Even asymptotic giant branch (AGB) stars, which are known for their slow outflows, still have $v_{\rm w}\gtrapprox 30~{}\mathrm{km}~{}\mathrm{s}^{-1}$ (Habing & Olofsson, 2004; Ramstedt et al., 2008). Therefore, one can take the $u\rightarrow 0$ limit of eq. (3) and adopt a simple analytic approximation, $a_{\rm NDF}\approx 8.18G\rho R_{0}$ (see also Li et al., 2020).

Similar to the DF case, the NDF direction stays along the direction of star to gas relative motion. While the stellar dynamic timescales are estimated similar to eq. (2), the estimation of hydrodynamic timescales approximately reads,

where the spatial scale $l_{\rm hyd}$ should stand for the distance from the star to the bow shock, and $v_{\rm w}$ is the stellar wind velocity. As one can easily observe, the inequality $\tau_{\rm hyd}\ll\tau_{\rm dyn}$ still holds in the NDF cases.

There are multiple types of stars in an OC, and their interactions with the ambient gas are different. Unless otherwise noted, we use the anti-friction recipes for main sequence stars, red giant branch (RGB) stars, and neutron stars (NSs), while we assume that white dwarfs (WDs) and BHs have no outflows and thus obey $a_{\rm DF}$ in eq. (1). Admittedly, WDs and BHs in binaries can accrete and have consequent activities (e.g., disk winds, decretions, and jets), yet detailed discussions on those phenomena only add to the complications and obscure the effects concerned, and are beyond the scope of the current paper. For main sequence stars and giants, the wind properties are calibrated for stars with solar metallicity. The mass-loss rates are approximated by Reimers formula (e.g., Kudritzki & Reimers 1978; see also Choi et al. 2016),

where we use the approximate scaling laws $L\propto M^{3.5}$ (Kuiper, 1938) and $R\propto M^{0.8}$ (Demircan & Kahraman, 1991), with $L$ as the stellar luminosity, $R$ as the stellar radius, and $M$ as the stellar mass. Such mass-loss rates are on the high end of typical values, yet stars in an OC are generally young and tend to have stronger stellar winds. The stellar wind velocity has significant variations and uncertainties from star to star. We adopt a simplified power-law fitting for main sequence stars,

which yields observed values of $v_{\rm w}=400~{}\mathrm{km}~{}\mathrm{s}^{-1}$ for solar-mass stars (e.g., Brooks et al., 2015), and $v_{\rm w}\approx 2000~{}\mathrm{km}~{}\mathrm{s}^{-1}$ for type O8 stars (e.g., Bernabeu et al., 1989). The uncertainties of wind properties are even greater for RGB stars; for simplicity, we assume all of them have $v_{\rm w,RGB}\approx 30~{}\mathrm{km}~{}\mathrm{s}^{-1}$ and $\dot{m}_{\rm w,RGB}\approx 3\times 10^{-9}~{}M_{\odot}~{}\mathrm{yr}^{-1}$ (Habing & Olofsson, 2004; Ramstedt et al., 2008). NSs can have significant outflows due to complex acceleration mechanisms, such as winds originating from young and intensely hot neutron stars are propelled by photons and neutrons near the stellar surface (Salpeter & Shapiro, 1981; Duncan et al., 1986). We assume $v_{\rm w}\approx 3\times 10^{4}~{}\mathrm{km}~{}\mathrm{s}^{-1}=0.1c$ ($c$ for the speed of light), and $\dot{m}_{\rm w}\approx 3\times 10^{-9}~{}M_{\odot}~{}\mathrm{yr}^{-1}$ (e.g., Kato, 1983; Paczynski, 1990; Mor et al., 2023).

Stars in an OC evolve for hundreds of millions of years before the cluster evaporates (Portegies Zwart et al., 2010; Krumholz et al., 2019; Krause et al., 2020). Given that the main sequence lifetime of an 8 $M_{\odot}$ star with solar metallicity is about 35 Myr (see stellar model PARSEC v1.2s ¹¹1https://people.sissa.it/$\sim$sbressan/parsec.html; Bressan et al., 2012; Tang et al., 2014; Chen et al., 2015), before the cluster evaporates all massive stars will have evolved off the main sequence. Very massive stars ($M\gtrapprox 25~{}M_{\odot}$) are not considered in this work. We assume that each star has finished its main sequence phase and entered its RGB stage at the time we consider significant gas-star interactions to occur. The evolution time of the RGB phase is approximated by a piecewise function in line with the solar metallicity and solar composition stellar models (Z=0.017, Y=0.279) of PARSEC v1.2s:

After the RGB stage, we assume that stars of $M<8\;M_{\odot}$ will undergo their AGB phase and end up as WDs. Because the AGB stage experiences massive outflows that lose about $50\%$ of the stellar mass within a relatively short period ($\Delta t\lessapprox 10^{5}~{}\mathrm{yr}$), the impact of AGB outflows due to anti-friction can be accounted as a pulse momentum injection, by integrating $\mathrm{d}v_{*}/\mathrm{d}t=a_{\rm NDF}$ assuming $u\ll 1$ and $\dot{m}_{\rm w}\approx\Delta M/\Delta t$. The increments regarding the velocity is then better estimated by the change in $v_{*}^{2}$,

where $m_{p}$ is the proton mass. $\Delta M$ and $\Delta v^{2}$ are the total mass-loss and the increment of $v_{*}^{2}$ during the period $\Delta t$, respectively. We assume that the consequent WD has $M_{\rm WD}=\mathrm{min}\{0.5M,M_{\odot}\}$ (here $M$ is for the progenitor RGB star mass), and $\Delta M$ is given by the decrease of mass before becoming a WD. For stars with $M>8\;M_{\odot}$, we estimate their $\Delta v^{2}$ similarly, assuming that the expansion period of the supernova ejecta is $\Delta t=10^{4}~{}\mathrm{yr}$, $v_{\rm w}=0.1c$, and $\Delta M$ is also deduced from $M_{\rm NS}=\mathrm{min}\{0.1M,1.5M_{\odot}\}$ and the progenitor mass $M$. Note that we consider this pulse momentum injection by AGB outflow in all types of gas-star interactions and intentionally ignore the “kicks” in NS (WD) momenta induced by the asymmetries of supernovae (mass-loss in WD progenitor), so as not to obscure the NDF by such kicks (while NDF only adds to the effects of the kicks). Other epochs of stellar evolution are ignored as they have negligible impact on the stellar cluster dynamics.

The cluster samples used in the studies on open clusters are usually generated in initial virial equilibrium that can be achieved through various numerical methods (Aarseth et al., 1974), such as N-body integration, the Monte Carlo method, and the fluid dynamical approach (e.g., Aarseth, 1973; Spitzer & Hart, 1971a, b; Larson, 1970a, b). In this work, the initial sample is generated by McLuster, which is an open-source tool that can set initial conditions for N-body simulations or generate artificial star clusters for direct study (Küpper et al., 2011). The virial equilibrium cluster samples generated using McLuster are detailed in Table 1. All generated samples use the broken power-law initial mass function over the stellar mass range $0.2<(M/M_{\odot})<25$ (Kroupa et al., 1993; Kroupa, 2002; Reid et al., 2002).

with the $W_{0}$ parameter in the King (1966) model, specifying model concentration, set to 5.0, and the half-mass radius along the line of sight set to $2.0~{}\mathrm{pc}$.

With the help of McLuster, different initial conditions can be set when generating cluster samples by modifying the default parameter values. For example, “-Q” sets the virial ratio to ensure the cluster is in virial equilibrium, “-S” adjusts the degree of mass segregation, and “-b” sets the binary fraction. The method used for generating samples with specified degrees of mass segregation and binary fraction is based on the approaches of Baumgardt et al. (2008); Kroupa (1995a, 2008). The fiducial sample is a 1500-star ensemble S1500 without segregation and binaries (see Figures 1–2), which will be used for studying different physical scenarios (§3 and 4). S500 and S3000 are two samples with different numbers of stars, used for comparison with S1500 in the fiducial model study. With the same stellar number as S1500, S1500-seg is a fully mass-segregated cluster, where the most massive star occupies the lowest energy orbit, and the least massive star occupies the highest energy orbit. Note that the previously mentioned samples do not include initial binaries, as introducing binaries would complicate internal interactions within the cluster. For example, collisions could occur between the outflows of stars in binary systems (e.g., Wang & Li, 2022). Such complicated interactions would require further study using hydrodynamic simulations, which is beyond the scope of this current paper and reserved for future works.

Following Kroupa et al. (2001), this work generates a sample with a $100\%$ binary fraction through random pairing (S1500b) to study simplified physical scenarios (considering only gas-star interactions). The binary fraction is defined as

where $N_{\rm bin}$ is the number of binary and $N_{\rm single}$ is the number of single stars. The method of identifying binaries follows Parker et al. (2009), which is calculating pairs with negative relative energy.

The fiducial model studies an OC immersed in the cold neutral medium (CNM), one of the phases of the interstellar medium (ISM) in the Galaxy, which has $\rho\approx 30~{}m_{p}~{}\mathrm{cm}^{-3}$ and $T\approx 100~{}\mathrm{K}$ (Draine, 2011). The OC then evolved for $\Delta t_{\rm evo}=200~{}{\rm Myr}$ within the bounding box. An important assumption is that star formation in the cluster is a single-age event, and once the cluster is formed, all the gas is expelled. Different models exploring physical parameters (especially those about the ambient gas properties) use the same initial fiducial sample S1500 (see Table 1-2). Simulations for each model, including the fiducial model specific to the CNM, are conducted with three types of gas-star interactions, indicated as “NDF”, “DF” and “None” as the suffixes of the models, to emphasize the modes of the gas-star interactions. Note that the stars that do not launch winds (e.g., WDs) always experience the DF effect even in the model marked as “NDF”.

Statistics about cluster sizes require a clear definition of the cluster center. This work uses the center of mass of the bounded stars, calculated through these two steps:

1. 1.

   Select all stars within the $L_{\rm box}=500~{}\mathrm{pc}$ simulation bounding box and calculate their total energy (the summation of kinetic and gravitational energy) with respect to the center of mass. Stars with negative total energy are considered bounded.

2. 2.

   Calculate the location and velocity of the center of mass for the bounded stars, which are used as the origin points in the configuration and momentum spaces in subsequent analyses.

Similar to the center of mass, the “center of light” is also desired for the analyses of stellar luminosity distribution, which is defined as the luminosity-weighted average coordinates for the bounded stars identified in the Step 1. Considering the observational conditions, unless specifically emphasized, all half-mass and half-light radii in this paper represent 2D values projected along the line of sight, calculated by averaging over three lines of sight as samples,

where $\{r_{1/2}^{xy},r_{1/2}^{yz},r_{1/2}^{zx}\}$ are the projected half-mass or half-light radii calculated on the $x-y,\,y-z,\,z-x$ planes, respectively.

For better presentations of fluctuations in the cluster characterizations, the uncertainties about statistical results in this work are estimated by the non-parametric bootstrap resampling method, which is a widely accepted and commonly used statistical method in astrophysics (see also Efron, 1979; Feigelson & Babu, 2012; Malofeeva et al., 2023). In detail, the bootstrap method generates datasets of the same length as the original data through resampling with replacement, repeated 1000 times. The standard deviations in each bin were then computed based on the results from these 1000 bootstrapped datasets.

Statistics directly related to observables are presented in the left column of Figure 4. All three models exhibit the same trend of dynamical evaporation, indicated by their increasing half-mass and half-light radii (in each time interval without star deaths), and decreasing velocity dispersions. The half-mass radii for the CNMF-NDF model are larger than the gas-free CNMF-None model by about $15\%$, while those for the CNMF-DF model are smaller by approximately the same percentage. Acceleration forces (NDF) tend to increase the total energy of stars, relocating them to shallower regions in the potential well and reducing the magnitude of the velocity dispersion. The deceleration forces (DF), in contrast, tend to result in larger velocity dispersions and smaller cluster sizes.

One can observe “cliffs” at $t_{\rm evo}\approx 39~{}{\rm Myr}$ in both the half-mass and half-light radii for the unsegregated NDF model in Figure 4. These cliffs are the consequences of both the expulsion of massive stars ($M\gtrapprox 7~{}M_{\odot}$) due to NDF acceleration , and the termination of massive stars on the main sequence as well. They are more pronounced in the half-light radius, as the luminosity is more sensitive to stellar mass ($L\propto M^{3.5}$). In contrast, the DF models slow down the massive stars more efficiently, delaying their evaporation and resulting in a slower increase in cluster size compared to the case without gas-star interaction.

To illustrate the influence of the star numbers, the results for samples S500 and S3000 with different numbers of stars are presented in the right column of Figure 4. The increasing cluster size and decreasing velocity dispersion over time are qualitatively similar to S1500, while the cluster size grows much faster with NDF in S500 and slower in S3000. With decreasing star numbers, the evaporation timescale comes closer to the crossing timescale ($t_{\rm evap}\propto N_{\rm cl}t_{\rm cross}/\ln N_{\rm cl}$; see also Binney & Tremaine 2008), making the cluster more susceptible to the acceleration caused by NDF, and vice versa.

The formation and evolution of binaries within OCs has long been a topic of concern (e.g., Kroupa, 1995b, a, c; Li & Mao, 2018). Figure 5 illustrates the evolution of the binary fraction. For both S1500 and S1500b, the three runs for the CNMF model do not exhibit significant differences in $F_{\rm b}$. For S1500b, binary systems are rapidly disrupted in the first 1 Myr, during which the soft binaries are destroyed. Subsequently, the clusters reach an equilibrium where the binary fraction becomes roughly stable. This behavior is consistent with Kroupa (1995b, a); Kroupa et al. (1999); Parker et al. (2009).

Actual stellar clusters may have different extents of mass segregation during their formation and early evolution, that more massive stars tend to be distributed closer to the cluster center (e.g. Chandrasekhar, 1943; Chandrasekhar & von Neumann, 1942). We study a completely segregated cluster sample S1500-seg generated by McLuster (see also Table 1 and Figure 6), and present the evolution of properties in Figure 7.

The differences in the half-light radii and half-mass radii between NDF and the other two models are more significant than for the fiducial model S1500. After $200~{}{\rm Myr}$ evolution, the half-light radii for the NDF model are about twice as large as in the DF and gas-free models. The NDF model “overturns” the segregation by raising more massive stars to larger distances due to stronger acceleration, while the DF and gas-free models only add to the segregation during the evolution. With more numerical experiments not presented here, we confirm the general trend that the more mass-segregated the clusters are initially, the more important the NDF becomes along the dynamical evolution of OCs.

The gas-star interactions–both NDF and DF–are the results of stellar outflows and motion. Although the CNMF model runs assume no bulk motion of the whole cluster, the real situation is usually more complicated. The centers of mass of an OC and a gas cloud can move on their own tracks when evolving in galaxies, and the relative bulk motion between the two objects should be expected. We hence introduce models CNM3 and CNM5 to study such bulk motions at center of mass speeds $v_{\rm c}=3~{}\mathrm{km}~{}\mathrm{s}^{-1}$ and $5~{}\mathrm{km}~{}\mathrm{s}^{-1}$, respectively (see also Table 2).

Due to the morphologies of the stellar wind-ambient interaction regions, the magnitudes of NDF acceleration decreases with faster stellar velocities (eq. 3). This relation qualitatively yields the phenomena in Figure 8, where the increase in total energy caused by NDF is slightly suppressed with higher $v_{\rm c}$. Regarding the increase in half-mass and half-light radii, a similar suppression is observed in Figure 9. Concerning velocity dispersion, the NDF suppression results in stars settling in deeper regions of the potential well, leading to larger velocity dispersion than in the CNMF model. While these patterns can be qualitatively anticipated for higher $v_{\rm c}$ velocities, detailed analyses should be based on cluster-specific calculations rather than extrapolations.

Outflowing stars are expected to be pushed further away from the cluster center once NDF is in action. For typical main sequence stars immersed in the CNM, however, such effects are generally obscured by the overall trend of OC dynamical evaporation. When an outflowing star travels at a low initial speed in a uniform medium, the asymptotic relation near $t\rightarrow\infty$ between the travel distance $l$ and the time $t$ roughly reads, using the approximate expression for eq. (3) and integrating $\mathrm{d}^{2}l/\mathrm{d}t^{2}\approx a_{\rm NDF}$,

In the meantime, the velocity scales as $v\propto t^{1/2}$ asymptotically. These estimations also give the lower limit of time required to drive a star out of the cluster even without the gravitational potential well. When the potential well is present, direct exclusion of a solar-mass main sequence star only via NDF is not possible.

For objects with intensive outflows, however, the accelerations are much more significant dynamically. Although $a_{\rm NDF}\propto v_{\rm w}^{1/2}$ increases sub-linearly with $v_{\rm w}$, the NS with $v_{\rm w}=0.1c$ can still escape from the stellar ensemble well before the evaporation of OCs. Figure 10 illustrates the efficient expulsion of NS (upper panel) by NDF (solid lines), compared to the DF (dotted lines) and gas-free (dashed lines) cases where the ordinary “evaporation” process of an OC and the increment of a star’s velocity by AGB-outflow-induced pulse are present. NS are driven out to $\gtrapprox 10^{2}~{}\mathrm{pc}$ from the center of the cluster within only about $50~{}{\rm Myr}$, shorter than the lifetime of a typical OC. After escaping from the OC potential, the $l\propto t^{3/2}$ scaling and eq. (14) holds semi-quantitatively for each of the NS at large distances ($\gtrapprox 10~{}\mathrm{pc}$). Assuming that the gaseous disk half-thickness is about $0.5~{}{\rm kpc}$ in the Galaxy, such expulsion process will eventually take about $100~{}{\rm Myr}$ to remove an NS from the disk. Such quick expulsions will likely lead to the scarcity of NS in OCs, and probably in the whole gaseous galactic disk. Figure 11 illustrates a reduced expulsion of NSs attributed to the attenuation of the NDF effect by the bulk motion. The NS expulsion is ubiquitous for gas-immersed OCs unless the cluster bulk motion is too fast relative to the gas. Previous works often attribute the NS exclusion from clusters to the “kicks” during asymmetric supernovae explosions (e.g., Fragione & Banerjee, 2020). With the NDF effects, however, NS exclusion can still take place without these “kicks”. Note that even in the case of intensive outflows from stars, the winds of the stars cannot combine to drive a wind from the OC. This is because the standoff distance of the contact discontinuity,

is much smaller than the distance between member stars (e.g., the minimum distance in the initial OC of the fiducial case is about $0.01~{}\mathrm{pc}$).

In contrast to NS, WDs do not launch outflows in our models. This causes the absence of NDF during their WD epochs, and the pulsive acceleration in their AGB stages is clearly insufficient to exclude them from the OCs. However, the deficiency of WDs in young OCs has been proved by observations (e.g. Kalirai et al., 2001). Once the WD outflows due to binary accretion processes are included (which requires detailed modeling in feedback and the non-isotropic outflow patterns, e.g. Li et al. 2020, and is also postponed to future works), WDs are also expected to experience expulsion in their dynamical evolution. To illustrate another possible exclusion mechanism of WDs, we set up a model that enables pulses in velocity for newly formed WDs with a magnitude of $2~{}\mathrm{km}~{}\mathrm{s}^{-1}$ and random directions, following Fellhauer et al. (2003). As one can see from Figure 12, such pulses can indeed cause the exclusion of WDs regardless to the condition of gas-star interactions.

As a star departs from its cluster, its path is computed considering the Galaxy’s external potential. Incorporating the Galactic potential in models of OC offers a more realistic scenario. Following Allen & Santillan (1991), the Galactic potential $\phi_{\rm Gal}$ is composed of three parts: the central part $\phi_{1}$, the disc $\phi_{2}$, and the halo $\phi_{3}$, i.e., $\phi_{\rm Gal}=\phi_{1}+\phi_{2}+\phi_{3}$. These three parts of the potential can be expressed in cylindrical coordinates ($R,z$) as follows:

where

The parameters for the above equations are $M_{1}=1.41\times 10^{10}~{}M_{\odot}$, $b_{1}=387~{}\mathrm{pc}$, $M_{2}=8.56\times 10^{10}~{}M_{\odot}$, $a_{2}=532~{}\mathrm{pc}$, $b_{2}=250~{}\mathrm{pc}$, $M_{3}=10.7\times 10^{10}~{}M_{\odot}$, $a_{3}=1.2\times 10^{4}~{}\mathrm{pc}$ and $r_{\mathrm{d}}=10^{5}~{}\mathrm{pc}$. The OC is positioned on a circular orbits at a speed of $220~{}\mathrm{km}~{}\mathrm{s}^{-1}$, specifically $(v_{\rm R},v_{\phi},v_{\rm z})=(0,220,0)~{}\mathrm{km}~{}\mathrm{s}^{-1}$, at a distance of $8.5~{}{\rm kpc}$ from the Galactic center.

The impacts of the gas-star interactions in this galactic scenario are illustrated by Figure 13. While the overall trend is consistent with Figure 4 that both the half-mass and half-light radii increase over time, the tidal fields dominate the dynamical evolution and generally dwarf the differences in the half-mass radii and velocity dispersions for the three gas-star interaction models. The half-light radii, which are more sensitive to massive stars, are greater for NDF and DF models compared to the gas-free case. Dinnbier & Kroupa (2020) studied the dynamical evolution of the tidal tails of OCs in detail and provided a semi-analytical model, concluding that such oscillations in the velocity dispersions solely depend on the local galactic dynamical properties, especially the orbital frequency and the epicycle frequency, rather than any cluster properties (e.g., the initial velocity dispersion or mass).

It is expected that the tidal field could be one of the significant effects of the evolution of internal dynamical features within OC. However, this work focuses on the effect of gas-star interactions. While incorporating the tidal field would render the scenario more realistic, it also complicates the physical mechanisms, making it challenging to analyze the effects of different mechanisms clearly. Therefore, the tidal field will not be considered in subsequent scenarios, meaning the external Galaxy’s potential will be excluded. In future work, the influence of both tidal fields and gas interactions on the dynamical characteristics of stars or clusters should be included.

In a virialized self-gravitating system, stars staying at larger distances from the cluster center generally move at slower velocities. When the NDF and DF accelerations are sufficiently low and gradual, their impact on the cluster configurations is mostly adiabatic, and the system remains virialized during the evolution. When the ambient gas density is raised to $3\times 10^{8}~{}m_{p}~{}\mathrm{cm}^{-3}$ in the model AGND, however, the NDF acceleration for a solar-mass star moving at $1~{}\mathrm{km}~{}\mathrm{s}^{-1}$ becomes $a_{\rm NDF}\approx 10^{-9}~{}\mathrm{cm}~{}\mathrm{s}^{-2}$, which is about $10\times$ the gravitational acceleration by a $10^{3}~{}M_{\odot}$ cluster at $10~{}\mathrm{pc}$. A cluster evolving under such strong NDF or DF accelerations no longer remains adiabatic; the dispersal takes place directly within the first $3~{}{\rm Myr}$ resulting in a quickly expanding spatial size and higher velocity dispersions (Figure 14). This phenomenon also appears in the models M7 and M5 that stand for typical molecular clouds, while the models for more diffuse gas (WNM0, WNM1) are qualitatively similar to the CNMF case.

Another effect that emerges in various ambient models is the stratification of stars. The acceleration by NDF on stars tend to increase their energy, which is similar to some “buoyancy” in self-gravitating OCs. As more massive stars launch more powerful stellar winds (§2.2), they tend to move quicker to the “surfaces” of OCs under NDF, and vice versa. We notice that such “buoyancy” will not make the young OCs “inversely segregated” initially, as the timescales required to fully “inversely segregate” the cluster (about $10^{7}-10^{8}~{}\mathrm{yr}$) (see e.g. Figure 15) is longer than the time period of star formation ($\lessapprox 10^{7}~{}\mathrm{yr}$) (Palla & Stahler, 1999; Palla et al., 2005).

Since the stellar luminosity scales as $L\propto M^{3.5}$, the half-light radii are considerably greater than half-mass radii in NDF-affected clusters. For those models that have dense ambients (e.g. M4, M5), the half-light radii are $\gtrapprox 3\times$ greater than the half-mass radii after only about $10~{}{\rm Myr}$, indicating an obvious radial stratification in terms of stellar types. When attempting to infer the age of an observed OC by the observations in astrometry and dynamics, one should take care of the possible encounters with dense clouds, whose effect in puffing up the OCs is similar to the intrinsic evaporation of clusters. Statistics on the radial distributions of mass and light in such clusters, therefore, may be helpful in reducing this type of parameter degeneracies.

According to eq. (3), the magnitude of NDF acceleration depends much more sensitively on $\rho$ than on $T$. Figure 16 illustrates the general trend that the NDF effects increase with denser ambient gas, raising the total energy $E$ and easing the evaporation of cluster stars. The scaling relation in eq. (3) is sub-linear with respect to $\rho$ ($\propto\rho^{1/2}$). Assuming a relatively invariant velocity distribution over the same period of time, one can infer the scaling relation $\mathrm{d}E/\mathrm{d}t\propto\rho^{1/2}$, and subsequently $\Delta|1/r|\propto\rho^{1/2}$. This pattern is qualitatively seen in Figure 14, yet we note that such simple scaling no longer applies quantitatively when there are significant increases in velocity dispersions.

Models M5 and M7 show that OCs in relatively dense ambients will quickly puff up or even disperse within $20~{}{\rm Myr}$, a timescale that is comparable or shorter than the crossing time of an OC over a dense molecular cloud. If one adopts the typical sizes (about $20~{}\mathrm{pc}$) of dense molecular clouds ($\rho\gtrapprox 10^{4}~{}m_{p}~{}\mathrm{cm}^{-3}$) (Cernicharo, 1991; Bergin & Tafalla, 2007) with an OC crossing time $\approx 40~{}{\rm Myr}$ (assuming bulk velocity $\approx 0.5~{}\mathrm{km}~{}\mathrm{s}^{-1}$), Model M5 indicates that a cluster should puff up to $r_{1/2}\gtrapprox 10~{}\mathrm{pc}$ after this crossing. Similar results also occur in Model M7, even if one limits the interaction time to the typical lifetime of dense molecular clouds (about a few Myr). The encounter with molecular clouds can lead to significant changes in its speed due to NDF or DF, scaled with respect to the data of solar-like stars,

For AGB stars with massive outflows (see also eq.(8) in §2.1), the increment in stellar velocity caused by NDF could raise further to $\Delta v_{*}^{\rm NDF}\approx 0.2~{}\mathrm{km}~{}\mathrm{s}^{-1}$ after $10^{5}~{}\mathrm{yr}$.

AGN disks are also considered as places where stars can form, whose lifetimes are at the order of $10^{0}-10^{3}~{}{\rm Myr}$ (Marconi et al., 2004; Martini & Weinberg, 2001). If a $10^{3}~{}M_{\odot}$ OC forms and stays inside the AGN disk for $\gtrapprox 5~{}{\rm Myr}$, it could evaporate almost completely before it is no longer considered as a cluster anymore, even if the destruction processes by orbital motion shears and tidal forces are taken into account. In other words, the likelihood that one can find OCs formed in AGN disks should not be significant.

As elaborated in §3.3 and 4.1, higher stellar velocities and lower ambient gas densities suppress the NDF acceleration. In various models, the expulsion of NS via NDF behaves accordingly. Figure 17 illustrates that, even for Model WNM1 whose ambient density is merely $3~{}m_{p}~{}\mathrm{cm}^{-3}$, such expulsion mechanism is still working, although the time required is considerably longer. This effect is, nevertheless, still more perceptible than the CNM5 case (see Figure 11). For Model WNM0, the expulsion mechanism weakened compared to WNM1, but all NSs still tend to move away from the OC.

The binary fraction corresponding to various models does not show any considerable changes compared to the fiducial CNMF model, except for the AGND model with DF only (Figure 18). The binary fraction increased drastically in the first $4~{}{\rm Myr}$ of evolution. Strong deceleration caused by DF shrinks the separations between stars and reduces stellar velocities, facilitating the formation of binaries. Several studies have explored this phenomenon. Rafikov (2013, 2016) proposed that gas-assisted inspirals efficiently merge supermassive BH binaries within a Hubble time. The NDF would nonetheless overturn this scenario, widening the binary separations. Wang & Li (2022) demonstrated in simulations of AGB star-outflowing pulsars that a dense and slow outflow exerts a positive torque on the binary, causing $\gtrapprox 10\%$ orbit expansion. Extrapolating these mechanisms to the OCs immersed in dense gas can cause a severe inhibition of binary occurrence rates.

While gas DF has garnered most of the attention, NDF has been relatively overlooked. Opposite to DF, NDF accelerates stars along their motion direction. Such difference in the effects between NDF versus DF can “flip” the signs of multiple physical mechanisms. For instance, one may expect that the OC evaporation is suppressed or even inhibited by DF when encountered with gas, while the actual scenario could be right the opposite due to NDF. In the densest gases including dark molecular clouds, AGN disks, and AGN tori (whose average column density is as high as about $10^{24}~{}{\rm cm}^{-2}$; see, e.g., Zhao et al. 2021), dispersal of OCs can be almost instant compared to their dynamical lifetimes as clusters. Interpretations about observations on cluster kinetics should take the possibility of NDF into account. Pflamm-Altenburg & Kroupa (2009) argued that massive stellar clusters over a certain critical mass can act as cloud condensation nuclei and accrete gas from the ambient medium. Such critical mass is required for the cluster to accrete gas and induce density instability in the ISM roughly satisfies $M_{\rm cl}/T_{\rm g}b_{\rm cl}>3.33\times 10^{3}~{}M_{\odot}~{}{\rm K}^{-1}~{}\mathrm{pc}^{-1}$ (where $M_{\rm cl}$ is the cluster mass, $T_{\rm g}$ is the gas temperature, and $b_{\rm cl}$ is the Plummer parameter; see also Plummer 1911). Admittedly, massive clusters can suppress the NDF effect by accreting from the ambient gas, which is also beyond the scope of the current study and will be addressed in future works. In the scenarios considered in this work, nevertheless, the cluster mass is insufficient to initiate this type of accretion.

Observation missions like Gaia (Gaia Collaboration et al., 2016) provide abundant astrometric data for stars in our galaxy, facilitating the study of OCs’ dynamical evolution (e.g., Maurya et al., 2023). Because of the radial stratification in stellar masses (§4.1), one can prospectively tell the effects of gas-star interactions by measuring the half-mass and half-light radii, as well as the velocity dispersions. The combination of high precision photometry with astrometry can further help the researchers to find out the radial distribution functions of different stellar types. These measurements can be a direct characterization of gas-star interaction history that are potentially important for identifying an OC’s “invasion” into dense gases during its evolution history.

Due to substantial outflows, NSs gradually drift away from their associated OC under the influence of NDF in all simulations. Conversely, WD stars influenced by NDF experience a pulse momentum injection from the main sequence to the AGB stage and subsequently to the WD stage. However, upon transitioning into the WD phase, WDs decelerate due to DF as they lack outflows. This mechanism ensures that part of WDs remain confined within the OC and experience conventional evaporation. If the kick received by a WD is strong enough, such as through an asymmetric mass loss process (Fragione & Banerjee, 2020), it is more likely to be expelled from the OC (see in Figure 12).

It is commonly assumed that kicks by asymmetric supernovae drive NS to high velocities relative to the cluster centers, and eventually expel them from OCs. Lai (2004) reviewed various physical mechanisms leading to kicks, and Bortolas et al. (2017) suggested that NSs are generally scattered away from SgrA* due to such kicks. Contenta et al. (2015) indicated that the natal kick of NSs can alter the cluster’s lifetime by almost a factor of 4. Admittedly, the kicks serve as an effective expulsion mechanism of NS. However, kicks only entail a one-time momentum injection, whose intensity and direction of momentum injection is stochastic and can well result in the NS’s deceleration. The NDF effects are, in contrast, a much more steady and robust mechanism that almost ensures the expulsion without gravitational scatterings. Such effects can further lead to the reduction of NS binaries and mergers within OCs, which is potentially relevant to the spatial distribution of NS-related mergers and gravitational wave events.

Due to the limitations in physical modeling and computation, this work does have some caveats and issues that should be addressed in future works. For example, simplistic isotropic stellar wind models are adopted for each type of star, which is suitable only for a fraction of NS. Many NS–especially those with accretion disks–may launch directed or bipolar outflows, which could affect the intensity and direction of the NDF forces (Li et al., 2020). The WD model in this study lacks wind and is susceptible to DF only. While stand-alone WD stars may not have outflows, accretion from their companions may cause disk winds and jets. This study ignores the time lag in forming stars within the OC and does not consider the interaction between star and star-forming gas. However, the expulsion of residual star-forming gas does not occur instantaneously after the completion of the star formation process in realistic scenario. Observations reveal that star formation efficiencies within OCs vary widely, ranging from several percent to 30 percent for dense clumps within molecular clouds (Lada & Lada, 2003; Higuchi et al., 2009; Megeath et al., 2016), and from 0.1 percent to a few percent for their associated giant molecular clouds (Evans et al., 2009; Murray, 2011). This implies that newly formed or forming stars (such as protostars and pre-MS stars, which have outflows) would interact with star-forming gas, resulting in the influence of NDF on the initial phase space distribution of stars within young OC. Moreover, the stellar winds from these stars may blow unprocessed gas out of the cluster, inhibiting the ongoing star formation process. Refined models in the future should take the complexities of stellar outflows and the star formation process into account. Besides, this work neglects the interactions between the outflows of stars in binary systems and between binary systems and the outflows of single stars. These effects will be investigated in future work using more refined hydrodynamic simulations.

Although our study considers the highly dense scenario of an AGN disk (AGND scenario), we have omitted multiple physical conditions relevant to AGN in order to isolate and emphasize the impact by NDF. Incorporating the shear effect is crucial for understanding stellar formation and evolution within AGN disks. However, this topic is specific to disks and beyond the scope of this paper, which should be addressed in a separate study. Future studies on AGN disks will require additional considerations, such as orbital motion and vertical stratification of gas. The limited thickness of AGN disks differs substantially from the assumed uniform gas ambient in this paper, especially considering the vertical density gradients. Furthermore, isolated stars would experience the DF effects in the highly dense environment of a galaxy’s circumnuclear regions. Future works shall also explore these possibilities of NDF-affected stars.

In addition, this work focuses on relatively small and sparse OCs with about $10^{3}$ stars. In high-density environments, particularly in AGN disks, super star clusters are often present instead of OCs, whose star formation efficiencies may be typically high. Investigating larger ensembles of stars such as globular clusters can impose additional challenges. Globular clusters contain about $10^{6}$ or even more stars, while each of the stars might still be dynamically important. Accurate treatments of frequent close encounters are beyond the capacity of ordinary orbital integrators with tree-based or particle-mesh-based gravity solvers, and special algorithms are required for sufficient accuracy. During close encounters, shock structures caused by gas-star interactions may collide, resulting in highly complex gas morphologies and interaction patterns. Understanding how such complex physical scenarios affect the dynamical evolution of stars within globular clusters could be explored once proper algorithms are prepared for complicated computations. Even without proper treatments of frequent close encounters in more massive clusters, the dynamical evolution of compact objects (NS, WD, and BH) are still qualitatively feasible within the current framework for globular clusters and super star clusters. Detailed discussions are nontheless beyond the scope of the current paper focusing on OCs, and a subsequent paper is being composed specifically addressing this issue.