{'fraction_non_alphanumeric': 0.03969478737997256, 'fraction_numerical': 0.03873028120713306, 'mean_word_length': 4.259497238191861, 'pattern_counts': {'":': 0, '<': 35, '<?xml version=': 0, '>': 11, 'https://': 2, 'lorem ipsum': 0, 'www.': 2, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 2, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
The internal proper motion kinematics of NGC346: past formation and future evolution

We investigate the internal kinematics of the young star-forming region NGC 346 in the Small Magellanic Cloud. We used two epochs of deep F555W and F814W Hubble Space Telescope ACS observations with an 11-year baseline to determine proper motions, and study the kinematics of different populations, as identified by their color-magnitude diagram and spatial distribution characteristics. The proper motion field of the young stars shows a complex structure with spatially coherent patterns. NGC 346 upper-main sequence and pre-main sequence stars follow very similar motion patterns, with the outer parts of the cluster being characterized both by outflows and inflows. The proper motion field in the inner ~10 pc shows a combination of rotation and inflow, indicative of inspiraling motion. The rotation velocity in this regions peaks at ~3 km/s, whereas the inflow velocity peaks at ~1 km/s. Sub-clusters and massive young stellar objects in NGC 346 are found at the interface of significant changes in the coherence of the proper motion field. This suggests that turbulence is the main star formation driver in this region. The similar kinematics observed in the metal-poor NGC 346 and the Milky Way star-forming regions suggest that the differences in the cooling conditions due to the different amounts of metallicity and dust density between the SMC and our Galaxy are too small to alter significantly the process of star clusters assembly and growth. The main characteristics of our findings are consistent with various proposed star cluster formation models.


INTRODUCTION
Regions of massive star formation (SF) are important to study for the insight they provide into the origin of the stellar mass distribution, the formation of massive stars (M > 8 M ), and the formation and evolution of associations and star clusters. Nevertheless, the process whereby stars form in giant molecular clouds is still a poorly constrained problem.
The first few million years, during which stars are still forming and the region still contains a significant amount of gas, set the conditions for the subsequent evolution of the stellar system. This initial stage is regulated by gas and stellar dynamics, stellar evolution, and radiative transfer. However, the complex interplay between these quantities is poorly understood by theory and lacks sufficient observational constraints (i.e. Elmegreen 2007;Price & Bate 2009;Portegies Zwart et al. 2010;Krumholz 2014;Kruijssen et al. 2019).
This shortfall of information leaves many essential and critical properties uncertain, such as the duration and the efficiency of the SF process, and hence the probability that a cluster will survive the rapid ejection of gas caused by its first supernova explosions (Carpenter et al. 2000;Lada & Lada 2003;Chevance et al. 2020). The internal kinematics of young star clusters and associations carries the signature of the process that led to the systems' formation. Radial velocity and proper motion (PM) studies can gather information on the dynamical state of gas and stars, and test different models of SF (Elmegreen 2002;Parker et al. 2014). Currently, theories describing the onset and development of SF can be divided into two competing scenarios: a rapid process, which proceeds on a dynamical timescale, comparable to the freefall collapse of the molecular cloud (e.g., Elmegreen 2000;Dobbs et al. 2011;Hartmann et al. 2012;Grudić et al. 2018;Jeffreson & Kruijssen 2018); and a slow one that persists for multiple freefall timescales (e.g. McKee 1989; Tan et al. 2006;Krumholz & McKee 2020) Astrometric measurements by ESA's Gaia spacecraft (Gaia Collaboration et al. 2016) have already provided spatial and kinematic clustering properties for several star-forming regions in the Milky Way, including Orion (Großschedl et al. 2018;Kounkel et al. 2018;Getman et al. 2019;Kuhn et al. 2019), Taurus (Luhman 2018;Galli et al. 2019), ρ Oph (Cánovas et al. 2019), Serpens (Herczeg et al. 2019), NGC 6530 , and IC 5070 (Kuhn et al. 2020).
In this work we extend this type of investigation for the first time to a low metallicity environment (Z = 0.2Z , Russell & Dopita 1992;Rolleston et al. 2003;Hunter et al. 2007), such as the Small Magellanic Cloud (SMC), to explore whether different global conditions (e.g., as large-scale dynamics and metallicity) can affect the duration of SF.
As the brightest and largest SMC star-forming region, NGC 346, located in the northern part of the galaxy bar, is one of the best-studied extra-galactic young clusters. In an attempt to better understand the formation mechanism and early evolution of NGC 346, here we take advantage of the Hubble Space Telescope's (HST) exquisite astrometric capabilities and longevity to measure the PM displacements of high, intermediate, and low-mass (∼ 8 to ∼ 1 M ) stars and infer the system's internal kinematics.
The left panel of Fig. 1 shows the footprint of the HST mosaic superimposed on the Digitized Sky Survey (DSS) image, while the right panel shows the HST color-composite image. Here NGC 346 appears as a vast agglomeration of bright white stars, still embedded in diffuse ionized gas (Hα and Nii, shown in blue), extending over ∼ 50 pc from the northwest to the southeast. The star-forming region contains a "plethora" of massive stars (Niemela et al. 1986;Massey et al. 1989;Evans et al. 2006;Dufton et al. 2019) that are ionizing the relatively isolated, large Hii region N66. O and B stars are embedded in an extended halo of pre-main sequence (PMS) stars (Nota et al. 2006), organized in several sub-clusters, clumps, and asterisms (Sabbi et al. 2007;Gouliermis et al. 2014).
In Cignoni et al. (2011) we studied the SF history of the region and found that SF started, with remarkable synchronization, about ∼ 6 Myr ago. Then, progressing inward, it peaked ∼ 3 Myr ago, and is now continuing at a lower rate. Similar conclusions were recently reached by Dufton et al. (2019).
Extended CO clouds are still associated with NGC 346 (Rubio et al. 2000;Müller et al. 2015;Neelamkodan et al. 2021), and Spitzer observations revealed a multitude of young stellar objects (YSOs, Simon et al. 2007;Sewi lo et al. 2013) associated with the sub-clusters and asterisms identified by Sabbi et al. (2007). Rubio et al. (2018) analyzed HK spectra of three YSOs and concluded that these are likely ∼ 25 − 26 M Class I type, suggesting  Table 1. List of the observations analyzed in this paper. The first and sixth columns report the proposal IDs, the observing dates are listed in columns two and seven, the used filters are in columns three and eight, the covered regions in columns four and nine. Exposure times are in columns five and ten.
that, at least in the central cluster, star formation is still ongoing.
The bright clump of red stars to the north of NGC 346 is the core of the star cluster BS90. NGC 346 and BS90 are considered non-interacting, and their proximity is likely a simple visual alignment (Bica & Schmitt 1995). Using the synthetic color-magnitude diagram (CMD) method initially introduced by Tosi et al. (1991), Sabbi et al. (2007) concluded that the cluster is 4.5 ± 0.1 Gyr old, with metallicity Z = 0.02, E(B − V ) = 0.08, distance modulus (m − M ) 0 = 18.9, and a total mass ∼ 8 × 10 4 M . The halo of BS90 partially overlaps in projection with NGC346, and the two clusters are embedded in the SMC field, making it impossible to separate the different stellar populations using only their coordinates, or the characteristics of their CMDs.
The paper is organized as follows: section 2 contains a description of the observations and the photometric analysis. We discuss the characteristics of the stellar populations found in the region in Section 3. We present PM measurements in Section 4, and in Section 5 we analyze the kinematics of the various populations. In Section 6 we focus our attention on NGC 346 internal kinematics, and We present our conclusions in Section 7.

OBSERVATIONS AND DATA REDUCTION
(HST)'s Wide-Field Channel (WFC) of the Advanced Camera for Surveys (ACS; pixel scale ∼ 50 mas pixel −1 ) initially observed NGC 346 through the filters F555W and F814W in July 2004 (GO-10248, PI A. Nota). To measure the star-forming region's internal kinematics, we repeated these observations in July 2015, using the same orientation and filter set (GO-13680, PI. E. Sabbi). In both epochs, images were taken at three different pointing positions, each covering an area of 200 × 200 . In the first epoch, for each filter, a single long exposure was taken at the nominal center of the cluster. The other two pointings were spaced by ∼ 1 off the cluster's center towards the northeast and the southwest, respectively. A four-step dither pattern was applied in each filter for both the northern and the southern pointings to the right. In both filters and epochs, we marked all the sources with positional errors < 1 mas and photometric errors < 0.1 in orange. Sources with photometric rms < 0.2 in both filters and with positional rms in F814W < 3 mas in both epochs are in blue.
to remove hot pixels, better sample the point spread function (PSF), and fill in the detector gap. In the second epoch, the central pointing consisted of threedithered exposures in F555W and two-dithered exposures in F814W. As in the first epoch, we acquired four exposures per filter for both the northern and southern pointing. However, no dither stepping was applied to F555W filter observations covering the southern portion of the mosaic. Table 1 summarizes the list of the observations. The observations analyzed in this paper can be downloaded from the Mikulski Archive for Space Telescopes (MAST) via 10.17909/y1cf-m826 We carried out the analysis of the data directly on the bias-subtracted, flat-field and charge-transfer-efficiency (CTE) corrected flc exposures produced by the stan- Figure 3. The mF 555W − mF 814W , mF 814W CMD of the stars detected in the NGC 346 region, with photometric rms < 0.1 in both F814W and F555W filters. Sources with positional rms < 0.1 mas/yr in both X and Y directions and in both epochs have been highlighted as follow: RGB stars in the magnitude range 17.9 < mF 814W < 20.4 and color mF 555W − mF 814W > 0.5 are shown in red, and UMS stars in the magnitude ranges 17.9 < mF 814W < 20.9 and color mF 555W − mF 814W < 0.1 are shown in blue. LMS stars fainter than mF 814W < 21.9 and to the left of the line connecting mF 555W − mF 814W = 0.7, mF 814W < 21.9 to mF 555W − mF 814W = 1.3, mF 814W < 24.0 are marked in green. PMS stars fainter than mF 814W > 21.2, to the right of the line connecting mF 555W − mF 814W = 0.75, mF 814W < 21.2 to mF 555W − mF 814W = 1.5, mF 814W < 24.0, with photometric rms < 0.2 in both F814W and F555W, and positional rms < 0.3 mas/yr in both X and Y are shown in purple. dard calibration pipeline CALACS v10.2.4. Compared to drizzled ( drc) images, flc data have the advantage of not being re-sampled, thus providing a more direct representation of the astronomical scene. However, flc exposures are still affected by geometric distortion. To take into account this effect, we created a distortion-free reference frame using the geometric distortion correction for the ACS/WFC detectors described in Anderson & King (2006) 1 , and we related the photometry and astrometry of each exposure to that frame.
We analyzed all images using the Fortran routine hst1pass 2 (Anderson et al. 2022, in prep.). The program performs a single pass of finding and measures each star in each exposure by fitting a library of spatially variable empirical PSFs 3 , ignoring any contribution from neighbors.
Seasonal and orbital thermal changes can cause variations in the optical path length of HST up to a few microns within the timescale of an orbit (Bély et al. 1993;Lallo et al. 2006). These changes affect the telescope's focus and translate to small but measurable differences in the PSF from one exposure to another, possibly affecting the precision with which source positions are determined. To take into account these effects, in each image, we measured the average residuals after the subtraction of the empirical PSF-library from the bright (S/N > 150), isolated (out to a 10-pixel radius), and non-saturated stars. We then applied, for each image, the correction to the PSF that minimizes its residuals and ran the photometry on all the sources with S/N > 3 and without companions within 3 pixels. This approach allowed us to create for each exposure a catalog of X and Y coordinates in pixels, magnitudes, quality of the fit (qfit, as defined in Anderson et al. 2008) and χ 2 of the fit.
The three pointing positions allowed us to observe the same stars both close and far from the readout amplifiers, thus verifying that in both epochs the CTE residuals were negligible both in stellar magnitudes and positions. We derived the photometric zero points by matching our observations to the NGC 346 photometric catalog published by Sabbi et al. (2007), and available for download from VizieR 4 We used a six-parameter linear transformation based on the cross-identified, well-measured, and nonsaturated stars to match all the exposures in a common photometric and astrometric reference frame for each epoch+filter combination. We obtained the initial reference frame positions by averaging the single-image positions. For both filters, the median 1-D rms in the position of the bright stars (18 < m < 21) is less than 0.5 mas (=0.01 pxl) in the first epoch and 0.7 mas (0.015 pxl) in the second epoch. Fig. 2 shows the position rms as a function of magnitude.

STELLAR CONTENT
The inspection of the CMD in Fig. 3 reveals a complex history of SF, characterized by episodes of different intensity. A detailed analysis of the region SF history is presented in Cignoni et al. (2011). Here, we will briefly discuss the properties of the CMD that will be useful for the subsequent analysis.
Above m F 814W 21, the CMD shows two welldefined sequences: the blue (−0.3 ≤ m F 555W − m F 814W ≤ 0.4) upper main sequence (UMS), and the red giant branch (RGB), with (m F 555W − m F 814W ≥ 0.8). The UMS includes intermediate and high-mass stars (M 3 M ) that formed between 2 and ∼ 600 Myr ago. The youngest component belongs to NGC 346, while the remaining stars are associated with the SMC field.
The RGB consists of evolved low-mass stars ( 2 M ) characteristic of a stellar population older than ∼ 1 Gyr. The compact group of stars around m F 814W 18.5 and m F 555W − m F 814W 1.0 is the red clump (RC) and corresponds to the core Helium burning phase for stars in the mass range between ∼ 1 and ∼ 2.5 M , with a main-sequence lifetimes ranging between ∼ 1 and 10 Gyr. The older stellar population includes both SMC and BS90 stars. A few bright (m F 814W < 18) evolved red super-giants connect the UMS to the top of the RC.
The majority of the stars below m F 814W ∼ 22 are low-mass main-sequence (LMS) stars. This sequence includes stars that formed between a few tens of Myr and several Gyr ago and belong either to the SMC field or BS90. NGC 346 stars below ∼ 3 M are still in the pre-main sequence (PMS) phase, and populate the faint (m F 814W > 20) cloud of red (m F 555W − m F 814W > 1.0) sources to the right of the LMS.

PROPER MOTIONS
The first step in measuring PM displacements is to define the reference system. We considered basing the absolute astrometric reference frame on the Gaia Early Data Release 3 (EDR3; Gaia Collaboration et al. 2021) catalog. However, below m F814W = 18 (which corresponds to our saturation threshold), there are fewer than 600 stars in EDR3 with position errors < 1 mas and PM errors < 1 mas/yr in common with our catalog. Furthermore, this list includes both the young and old stellar populations, that could have different relative motions. Thus, to create reliable astrometric reference frames, we would have had to project the stars to the positions they had in 2004. Given that the average PM uncertainty for these stars is 0.9 mas/yr, the reference frame of the first epoch would not meet the required astrometric accuracy. We, therefore, decided to follow the same approach used by Bellini et al. (2014) and to create a network of reference sources directly from our catalog. The PM displacement of each source is then measured relative to this reference network. We measured the PM of NGC 346 stars with respect to the RGB stellar population since the latter can be easily identified in the CMD (Fig. 3) using color and magnitude selections. We used the F814W exposures acquired in the first epoch to determine the average position of the RGB stars and create the reference network. In total, we identified more than 1400 RGB stars with photometric error < 0.1 magnitudes in both filters and epochs and with position errors < 1.0 mas in both epochs. We then applied the six-parameter global transformations needed to translate stellar positions in each exposure onto the reference frame system. To limit the number of mismatches, we considered only stars whose coordinates in the master-frame match to within 2.5 pixels (0.125 ).
To mitigate the effect of small, uncorrected systematics in, e.g., the HST geometric-distortion solutions and PSF models, we applied local corrections derived from the average residuals between the transformed positions of the closest N=30 RGB reference stars. The median distance of the furthest reference star is 600 pixels (∼ 30 ).
For each epoch, we computed the final position of each star as the sigma-clipped mean of the transformed positions of all images in the F814W filter at that epoch. We repeated the same analysis for the data taken in the F555W filter. In each filter, we defined the PM of each star as the difference between its position in the second and the first epoch, divided by the 11 years temporal baseline. We estimated the PM errors by adding the positional errors in each epoch in quadrature and dividing by the temporal baseline.
For the analysis of the PMs, we considered only sources with photometric errors < 0.1 magnitude in both filters and epochs and with X and Y positional errors < 1.0 mas in both epochs. As shown in Fig. 2, only sources in the magnitude range 18 < m < 23 met both our astrometric and photometric requirements. Fig. 4 shows the PM 1-D rms of the PMs as a function of magnitude. PM measurements in the F555W filter are affected by larger errors because of the smaller number of dithers in the second epoch's southern pointing, and the higher background level, caused by the ionized gas. Therefore, we decided to measure PMs primarily using the F814W filter, although we verified that the F555W observations provide consistent results. After the selections in photometry and position errors, the average PM rms in the F814W filter is 0.06 ± 0.02 mas/yr and 0.07 ± 0.02 mas/yr in the F555W. For consistency with previous studies of NGC 346 (e.g. Sabbi et al. 2008;Cignoni et al. 2011;Gouliermis et al. 2014), in this work we adopted 60.4 kpc as the nominal distance of the system. Thus, the average PM rms 5 in the F814W filter is 17.2 ± 5.7 km/s and 20.1 ± 5.7 km/s in the F555W Figure 7. Distribution of low-mass stars, as selected from the RGB and LMS VPDs. Upper-Left Panel : mF 555W − mF 814W , mF 814W CMD centered on the RGB evolutionary phase. Stars that in the RGB VPD are within 0.03 mas/yr from the over-density peak at µx = −0.04 mas/yr and µy = −0.02 mas/yr are shown in purple, while stars more than 0.08 mas/yr away from the peak are in pink. Upper-Middle Panel : Spatial distribution of the sources found within 0.03 mas/yr from the over-density peak in the RGB VPD superimposed on the F555W HST image. Upper-Right Panel : Spatial distribution of the RGB stars found more than 0.08 mas/yr away from the over-density peak in the RGB VPD. Lower-Left Panel : Portion of the mF 555W − mF 814W , mF 814W CMD centered on the LMS evolutionary phase. Stars that in the LMS VPD are within 0.03 mas/yr from µx = −0.04 mas/yr and µy = −0.02 mas/yr are shown in dark green, while stars more than 0.08 mas/yr away are in light green. The spatial distribution of the two groups of stars is shown in the Lower-Middle Panel and Lower-Right Panel respectively.
filter. While the PM uncertainties for most individual stars are too large to resolve their internal motion within the cluster, the averaged PMs for a population of stars or over a spatial subregion of the cluster are accurate enough to study the internal cluster kinematics, since averaging decreases the uncertainties by N −1/2 .
The vector-point diagram (VDP, Fig. 5) of the PMs in the NGC 346 region shows a complex structure with indications of multiple peaks in the higher density part of the plots and an elongated structure toward the lower right corner. It is clear from the inspection of Fig. 5, that this plot alone is not sufficient to separate the different populations.

KINEMATICS OF THE STELLAR POPULATIONS
The internal dynamics of the SMC is strongly affected by its recent history of interaction with the Large Magellanic Cloud, with the SMC being in the process of tidal disruption (Zivick et al. 2018). The numerous attempts to derive the internal motions of the various SMC components using either line-of-sight (Stanimirović et al. 2004;Harris & Zaritsky 2006;Evans & Howarth 2008;Dobbie et al. 2014) or PM (Kallivayalil et al. 2006(Kallivayalil et al. , 2013van der Marel & Sahlmann 2016;Zivick et al. 2018) measurements have thus far provided contradicting results, highlighting how, especially on small scales, younger and older stars in the field of the SMC may have different kinematics.
To better understand the behavior of the various stellar populations in the NGC 346 region, in the following analysis we will take advantage of the fact that for each star we have three primary bits of information: • The position projected on the sky (Fig. 1); • The position in the CMD (Fig. 3); • The position in the VPD (Fig. 5).
As a starting point, we used the CMD to separate RGB, LMS, and UMS stars, and then analyzed the motion of the three stellar components separately. The VPDs for the three populations are shown in Fig. 6.

Low-Mass Stars in the SMC Field
In this section we will examine the motion of the stars with mass M 2 M in the SMC field. As discussed in section 3, this group of stars can be divided in RGB and LMS stars. The left panel of Fig. 6 shows the VDP of the RGB stars. Since the reference frame and the transformations are based on the RGB stars, the weighted average of this population in the VDP is by construction centered on zero (µ x = −0.0061 ± 0.0003 mas/yr, µ y = −0.0028 ± 0.0002 mas/yr).
The properties of the RGB stars selected using their positions in the VPD of the PMs are shown in the three upper plots of Fig 7. In the CMD (upper-left panel) sources within 0.03 mas/yr (corresponding to a velocity dispersion of ∼ 8.5 km/s) from the overdensity peak found at at µ x = −0.04 mas/yr and mu y = −0.02 mas/yr are characterized by a small color dispersion. The RC of this populations is compact and well defined, suggesting that this group of stars has similar ages, chemical composition, distance, and it is affected by comparable dust extinction. Spatially, these sources are all concentrated around the nominal position of the 4.5 Gyr old BS90 cluster, as shown by the density contours in the upper-middle panel.
In contrast the stars that are found more than 0.08 mas/yr away from the density peak are likely part of the old field of the SMC. In this case the RC shows a broader dispersion in magnitude, suggesting that the population could be distributed along a broader line of site, and include stars that formed over a larger period of time. This is further supported by the fact that the brighter portion of the RC terminates in the V-like shape formed by the red and blue edges of the helium-burning phase typically observed in the CMDs of stellar populations younger than 1 Gyr, suggesting that this component formed stars for several billion years. These sources cover the entire studied region (as shown in the upper right panel). The impact of BS90 on the PM measurements is negligible: in fact after removing the RGB stars found within 75 from the center of BS90, the weighted average of the RGB in the VPD remains µ x = 0.0020 ± 0.0001 mas/yr, µ y = 0.0002 ± 0.0002 mas/yr.
The middle panel of Fig. 6 shows the PM VPD of the LMS stars. In addition to stars from the old field and BS90, this plot also includes the SMC's younger (∼ 50 − 100 Myr to 1 Gyr) component. The VPD is centered at µ x , µ y = 0, 0 mas/yr, but the presence of BS90 can still be recognized in the second higher density contour. The properties of the LMS stars are highlighted in the three lower panels of Fig 7. As for the RGB stars, the sources found within 0.03 mas/yr from µ x = −0.04, µ y = −0.02 are concentrated around the center of BS90 and in the CMD are distributed along a tight and well-defined sequence, suggesting that they all belong to a coeval stellar population, with negligible distance variation. On the contrary, the remaining stars show a uniform spatial distribution and a much broader LMS, as expected for field stars, affected by different reddening, distance, and possibly metallicity. These results imply that BS90 has a velocity of v x = −11.5 ± 0.5, v y = −5.7 ± 0.7 km/s (corresponding to v RA = 5.34 ± 0.70, v Dec = 11.67 ± 0.50 km/s) relative to the more broadly distributed field RGB population.

Intermediate and High-Mass Stars
The distribution of the UMS stars in the VPD (Fig. 6, right panel) is quite different from those observed for the RGB and LMS stars. This is the only diagram that shows an elongated structure in the upper left part (µ x −0.025 mas/yr, µ y 0.03 mas/yr), associated with the NGC 346-N66 system.
To better separate NGC 346 stellar content from the young SMC field, we used the isophotal contour level 5.48×10 −17 erg/cm 2 /Å/s (Fig. 8 -left panel). The central panel shows the VPD of the stars outside the isophotal level. The majority of these stars most likely belong to the SMC field. The right panel shows the VPD of the stars inside the isophotal contour. The plot here shows two separate peaks, embedded in an extended halo. The majority of these sources likely belong to NGC 346.
The VPD of the SMC young field is clearly shifted to the right of NGC 346 and peaks at µ x = 0.043 mas/yr, µ y = −0.040 mas/yr, with an extended and broad tail towards the lower-right part of the diagram. Contrary to what is found by Zivick et al. (2018), both the Kolmogorov-Smirnov test and the Cramér-von Mises criterion reject the hypothesis that young (UMS) and old (RGB) field stars have similar kinematics. The discrepancy between Zivic's and our results is likely due to the recent star formation history of the NGC 346 region. Cignoni et al. (2011) reported an elevated star formation rate over the past 100 Myr compared to the rest of the SMC. We suggest that the elevated and localized activity occurred in the past ∼ 100 Myr likely dominates the kinematics of the young field stars, and is the culprit for the observed difference.  The distribution of the PMs appears double-peaked and surrounded by an extended halo. To confirm that the two peaks are not an artifact of a too aggressive binning and take into proper account measurement errors, we ran an extensive Markov chain Monte Carlo (MCMC) Bayesian inference comparing a single Gaussian with the sum of two and three functions. Both the Akaike (AIC; Akaike 1974), and the Bayesian information criterion (BIC; Schwarz 1978) favor three Gaussians (two narrow peaks + a broader component) in the Y-direction. In the X-direction, the two narrow peaks are almost aligned; thus, we tested only the one and two peak hypotheses. The favored results are shown in the two histograms in Fig. 9.
Although the test run above is not a formal fitting of the distribution of motions, the parameters of that define the Gaussian functions favored by the Bayesian inference can be used to select the stars that likely belong to the different moving groups. In particular, Fig. 10 shows the spatial distribution of the stars found within 0.5 σ from the Fig. 9 Gaussian curve peaks. The com- Figure 10. Spatial distribution of the stars found within 0.5 σ from the three Gaussian peaks found in the side histogram of Fig. 9. The left panel shows the distribution of the stars associated with the broader green Gaussian, the central panel the distribution of the stars associated with the narrow magenta Gaussian, and the right panel the distribution of the stars found under the blue Gaussian. The arrows indicate the average motion of the stars found in the higher density regions, after the mean motion of NGC 346 has been subtracted. ponent characterized by the broader velocity dispersion (σ x2 = 18.61 ± 2.71, σ y2 = 12.27 ± 2.29 km/s) is also the most extended and covers the entire NGC 346 complex (left panel). With the respect to the NGC 346 systemic velocity, these stars are moving from northwest to southeast with a mean velocity v x2 = 7.80 ± 1.87, v y2 = −12.47 ± 2.16 km/s. The stars associated with the magenta (central panel) and the blue (right panel) Gaussian curves are clustered within the inner 8.8 pc from the photometric center of NGC 346. With respect to the NGC 346 median motion, most of the stars in the magenta component seem to move along an arch from South to North. After the subtraction of the NGC 346 systemic velocity, the mean velocity in the Y-direction is v y3 = −1.37 ± 0.23 km/s, Figure 12. PM VPD of the PMS stars, measured in the F814W filter. The median error in the µx and µy directions is shown in the upper right corner and the velocity dispersion σ y3 = 3.09 ± 0.30 km/s. The blue component moves in the opposite direction to the green one, from the center to the northwest. The mean velocity in the Y-direction is redshifted compared to NGC 346 median velocity by v y1 = 10.16±0.38 km/s, with velocity dispersion σ y1 = 7.30 ± 0.33 km/s. Combined, the two components along the X-direction are on average blueshifted compared to NGC 346 systemic velocity by v x1 = −2.19 ± 0.13 km/s and σ x1 = 6.48 ± 0.20 km/s. Fig. 10 reveals a complicated internal kinematics for NGC 346. The system appears dynamically hot, with large velocity dispersion and sudden changes in directions.
The comparison between the VPDs of the SMC young field and NGC 346 UMS stars (Fig. 8) shows a certain level of overlap between the motions of the two populations. To address the impact of field contamination, we measured the stellar density of UMS stars between 17.9 < m F 814W < 20.4, that meet our photometry and astrometry requirements, at a distance greater than R > 150 (corresponding to R 44 pc) from the center of NGC 346. This yield a stellar density of < 0.01 star per squared arcsecond. We then divided the region inside the isophotal contour shown in Fig. 8 in a 200 × 200 pixel grid, and we measured the density of the UMS stars in each cell. Assuming a uniform distribution and 100% completeness for the UMS stars of the SMC young field, we estimate that the field contamination in the various cells ranges from 5 to 30%, with an average value of 12%. We noticed that this is likely an upper limit, since the NGC 346 region is likely more incom-plete, due to the higher crowding and background level, than the outer parts of the image.
We statistically removed the contamination of the SMC young field from NGC 346 systemic velocity with the respect of the RGB, which became v sys,x = −9.87 ± 0.23 km/s and v sys,y = 11.82 ± 0.37 km/s, and in equatorial coordinates it corresponds to v sys,RA = −11.64 ± 0.37 km/s and v sys,y = 10.07 ± 0.23 km/s. As shown in Fig. 11 (left panel), after statistically removing the field contamination, the center, and possibly the entire upper part, of the NGC 346 complex, seems to rotate clockwise, once the systemic velocity of NGC 346 is subtracted. The lower half of the system appears to move away from the complex as if NGC 346 is breaking into two parts.

Pre Main Sequence Stars
Because of NGC 346 young age, stars below ∼ 3 M are still in the PMS evolutionary phase. With the exception of a few small clumps, these stars are mainly concentrated within the isophotal contour shown in Fig. 8.
Most of the PMS stars do not meet our PM selection criteria, likely because of the elevated level of stellar crowding and high background. However, in Fig. 3 CMD, PMS stars fainter than m F 814W > 21.2 represent the cleanest sample of NGC 346 stars. We therefore decided to relax the selection criteria (blue stars in Figs. 2 and 4) to perform a second, independent, analysis, although affected by larger uncertainties, of NGC 346 internal kinematics. Fig. 12 shows the VPD for the PMS stars with photometric rms < 0.2 in both the F555W and F814W filters, and with positional rms in F814W < 3 mas in both epochs (purple sources in the CMD shown in Fig. 3). As in the case of NGC 346 UMS stars, the plot shows two clear peaks, surrounded by an extended halo, in the upper left part of the diagram. A third, weaker peak is visible towards the lower right part of the plot.
As we did for NGC 346 UMS stars, we divided the PMS stars in a 200 × 200 pixel grid and for each cell we measured the average motion of the PMS stars, corrected for the cluster's systemic velocity. Since this group of stars was selected using the CMD, where the field contamination is negligible, no further correction was needed. Although the plot appears noisier than in the case of the UMS stars, the two distributions are qualitatively very similar. Both diagrams suggest that the upper part of the NGC 346 is rotating. The southern end of the system, on the other end, displays an entirely different motion, resembling an outflow. We derived for each UMS star in NGC 346 the radial (v ρ ) and tangential (v θ ) velocity components to study the profile of the cluster mean velocity and velocity dispersion. To start, we consider as the rotation center NGC 346's photometric center (X = 3501.43, Y = 3568.69, corresponding to R.A. = 14.7706984, Dec = −72.1775412, ICRS), whose position in Fig. 11 is marked by the green cross symbol. We subtracted the system's mean PM from each star, and then we derived the v ρ and v θ components for each source. In this reference system, within the inner R ∼ 11 − 13 pc, both v ρ and v θ ∝ r, with fitted slope ∼ 0.3 and ∼ 0.2 km/s/pc, respectively, suggesting that the star-forming regions is rotating as a rigid body, and expanding.

Velocity and Velocity Dispersion Profiles
Given the evidence for rotation, we then proceeded to identify the center of rotation, which need not coincide with the photometric center. We looped through a 20 × 20 pixel grid, and for each grid element, after having subtracted the mean PM of the stars found within 10 pc, we derived the sigma-clipped mean v θ value of the stars found within the same radius. We considered as center of rotation the grid point that gave the largest |v θ |, found at RA = 14.7575892, Dec = −72.1779791, ICRS (X = 3502.35, Y = 3845.09), ∼ 4.5 pc west of the photometric center. In Fig. 11 the position of the rotation center is marked with a red X symbol. Fig. 13 shows v θ and v ρ as a function of distance from NGC 346's rotation center. Despite the large dispersion, v θ (upper-right panel) appears to increase between ∼ 0.5 and 3.5 pc (where it reaches the maximum value v θmax = −3.2 km/s) at 0.2 km/s/pc, and then progressively decrease to almost zero at 10 pc (v θ10 pc = −0.5 km/s).
Within the inner 4 pc v ρ is negative and almost constant, indicative of an inflow in the inner regions. Instead, between ∼ 3.7 and ∼ 8 pc it increases with the distance with fitted slope ∼ 0.2 km/s/pc, indicative of expansion of the outer regions (Fig. 13, lower-right panel). Similar expansions have been observed also in some young Milky Way star clusters  and predicted by numerical simulations of the stellar feedback impact on young star cluster evolution (Grudić et al. 2022).
The combination of rotation and inflow in the central regions is suggestive of inspiraling motion. Fig. 14 shows the spatial distribution and direction of motion of NGC 346's individual UMS stars in the new reference system. In this plot the majority of the stars above Y = −500 pixels is spiraling from the north east towards the rotation center. A logarithmic spiral is overdrawn to guide the eye. The distribution of YSOs and sub-clusters aligns well with this motion.
We speculate that in proximity of the center, the stream of spiraling stars bend behind NGC 346 and re-emerge below Y = −500 as one of the two tails of stars that are moving toward the southeast, away from NGC 346. The dashed red line in Fig. 14 suggests the possible path followed by these sources.
A narrow stream of stars (indicated by the dashed blue line) and YSOs seems to connect the bright small cluster to the east of NGC 346 (Sc-13 in the classification by Sabbi et al. 2007) to the outer tail. These streams, inflows, and outflows of stars are reminiscent of the filament-like structures that are feeding the growth of star-forming regions in the hierarchical collapse models.
Hydrodynamic simulations of turbulent molecular clouds predict that rotation is a common characteristic of embedded massive (> 1000 M , Lee & Hennebelle 2016;Mapelli 2017;Ballone et al. 2020) star clusters, with the rotation curves increasing from the center, peaking approximately at 1-2 half mass radii, and then decreasing in the outer regions, where they can become slightly retrograde (e.g. Tiongco et al. 2017). Over time, the rotation of star cluster is expected to decrease because of stellar mass loss and two-body relaxation (Einsel & Spurzem 1999;Ernst et al. 2007;Kim et al. 2008;Hurley & Shara 2012;Tiongco et al. 2017).
Evidence for internal rotation in NGC 346 has been also found over the central ∼ 1 arcmin 2 by Zeidler et al.
(2022 -submitted) using line-of-sight velocities obtained with the ESO Multi-Object Spectroscopic Explorer (MUSE, Bacon et al. 2010). Similarly to what we observed in the PM distribution, the MUSE dataset shows two statistically significant velocity groups, and a possible third one. Like the magenta and blue groups shown in Fig. 10, one velocity group is more centrally concentrated, while the other one appears elongated along the northwest-southeast direction. As in the case of the PMs, the more centrally concentrated group is rotating with an angular velocity of Ω = −0.41 ± 0.07 Myr −1 , which translates to v 5 pc rot,RV = −1.98 ± 0.34 km/s and v 10 pc rot,RV = −3.95 ± 0.67 km/s, at a radial distance from the center of 5 and 10 pc, respectively, comparable to the results of this work. Unfortunately the MUSE field of view is too small to provide any information about the motion of the NGC 346 outskirts.

DISCUSSION & CONCLUSIONS
The star-forming complex NGC 346 presents several challenges to cluster formation models. The sound crossing time is t S = 50 pc/c S ∼ 2 × 10 8 Myr; thus, the two edges of the NGC 346 parental cloud were not in thermal contact with each other on a freefall collapse timescale and should therefore present considerable age differences. However, all NGC 346 stars fall in the age range 2-6 Myr (Evans et al. 2006;Sabbi et al. 2007;Cignoni et al. 2013;Dufton et al. 2019), raising the question of how star formation developed in such a synchronized manner over such a wide scale. The Upper Scorpius OB association (Preibisch & Zinnecker 1999), the Serpens molecular clouds (Herczeg et al. 2019), and the NGC 7000 and IC 5070 Nebulae (Kuhn et al. 2020) pose similar challenges to all the models of cloud collapse that do not include an accelerating star formation rate.
The spatial and kinematic clustering in young stellar systems can provide valuable constraints on the conditions that led to the onset of star formation (e.g., Elmegreen 2002;Parker et al. 2014;Kuhn et al. 2019). Therefore, we used two epochs of deep HST/ACS observations with a temporal baseline of 11 years ( 10.17909/y1cf-m826) to measure the PM displacements of the stars in the star-forming region NGC 346, and reconstruct the cluster kinematics.
• Our analysis of NGC 346 stellar PMs reveals a complex pattern of inflows and outflows. The starforming region's inner ∼ 10 pc is rapidly rotating, following a spiraling movement towards the center from the outer north-east end. Increasing rotation with the distance from the center is commonly found in the most up-to-date hydrodynamic simulations of turbulent molecular clouds (Lee & Hennebelle 2016;Mapelli 2017;Ballone et al. 2020).
• NGC 346 is characterized by a clumpy structure, with at least 15 different sub-clusters and asterisms (Sabbi et al. 2007;Hennekemper et al. 2008), that often host recently formed YSOs (Simon et al. 2007;Sewi lo et al. 2013). When we compared the position of sub-clusters and massive YSOs with NGC 346 motion patter, we found that overdensities tend to concentrate at the interface of significant changes in the coherence of the motion field, where we expect to see substantial gas friction and compression. This behaviour appears to be in good agreement with the predictions of rapid star formation in turbulent molecular cloud scenarios (e.g., Klessen et al. 2000;Bonnell & Bate 2002;Bonnell et al. 2003).
• NGC 346 shows a large velocity dispersion, and it is increasingly expanding with the distance from the center. Similar behavior is an expected outcome, for example, of global hierarchical collapse models (Vázquez-Semadeni et al. 2017, where the collapse of the parental cloud starts as small-scale events within a larger structure and culminates after a few Myr as a series of filamentary flows that accrete onto massive central clumps. These "river-like" structures drive the system accretion motions "down the gravitational potential" and are responsible for the late appearance of massive stars and their segregation in the central clumps. In agreement with this picture, the PMs studied in this paper revealed longitudinal inflows and outflows that extend over tens of parsecs. These motions could be at the origin of NGC 346's elevated degree of mass segregation reported by Sabbi et al. (2008).
• The inspection of the motion field of both UMS and PMS stars indicates that the upper part of NGC 346 is rotating along a coherent spiral-like pattern (Fig. 14). Both Cignoni et al. (2011) and Dufton et al. (2019) noted that the massive stars in the center of the system likely formed a couple of million years later than the outer part of the cluster. We speculate that the spiraling motion has been feeding this more recent episode of SF.
• On a broader scale, our PMs analysis also showed a difference in the relative motions of young (< 600 Myr) UMS and old (> 1 Gyr) RGB field stars. From the analysis of line-of-sight velocities of stars and gas Stanimirović et al. (2004); Harris & Zaritsky (2006); Evans & Howarth (2008); Dobbie et al. (2014) proposed for the SMC possible age-dependent internal kinematics. On the other hand, Zivick et al. (2018) did not find significant differences between the PMs of UMS and RGB stars in 30 HST fields. We believe that the elevated SF rate in the NGC 346 region over the past ∼ 100 Myr is the main reason behind the difference between our and Zivick's findings. In particular, our results show how on a small scale, the specific recent SF history of a region likely dominates the local kinematics. This test highlights the need for a uniform-coverage high-precision astrometric study of the entire SMC to reconstruct the galaxy's internal dynamic and evolution through its history of interactions with the LMC and the MW.
Our new measurements for the NGC 346 region represent the first attempt to infer the conditions that induced star formation in a metal-poor giant molecular cloud. We demonstrated that the velocity field of the star-forming region is complex and well matches the conditions predicted by recent hydrodynamic simulations and hierarchical collapse models. The similarities be- Figure 14. Spatial distribution of NGC 346 UMS stars. The arrows indicates the direction of the motion with the respect of the rotation center, while their length is proportional to their velocity. As in Fig. 10 the yellow stars mark the position of the massive YSOs. The red logarithmic spiral, described by the equation x = cos(θ)e 0.21θ y = sin(θ)exp 0.21θ , is not a fit, and has been superimposed only to guide the eye. The dashed line connect the stars in the center to the tail of stars that appear to escape NGC 346 from the south west. The dashed blue line highlight the stream of stars that connect the small cluster to the west to the southern tail.
tween our results and those found for star-forming regions in the Milky Way suggest that the differences in the cooling conditions due to the different amounts of metallicity and dust density between the SMC and our Galaxy are too small to alter significantly the process of star clusters assembly and growth. By extending this type of study to a broader range of star formation rates and stellar densities and by comparing the results to theoretical predictions, we can gain new insights into star formation and open a new pathway toward validating or repudiating specific models.
We thank an anonymous referee for carefully reading this manuscript and helping us improve the quality of the paper. We are grateful to M. Gieles, L. Strolger, and R. Klessen for the useful suggestions and discussions. The HST observations used in this paper are associated with program Nos. 10248 and 13680. The specific observations analyzed in this paper can be accessed from MAST at the Space Telescope Science Institute via 10.17909/y1cf-m826. Support for the program 13680 was provided through a grant from the Space Telescope Science Institute. This work is based on observations obtained with the NASA/ESA Hubble Space Telescope, at the Space Telescope Science Institute, which is operated by the Association of Universities for Research in Astronomy, Inc., under NASA contract NAS 5-26555. This work has made use of data from the European Space Agency (ESA) mission Gaia (https://www. cosmos.esa.int/gaia), processed by the Gaia Data Processing and Analysis Consortium (DPAC, https://www. cosmos.esa.int/web/gaia/dpac/consortium). Funding for the DPAC has been provided by national institutions, in particular the institutions participating in the Gaia Multilateral Agreement.