The Unusual Initial Mass Function of the Arches Cluster

Astrometry and photometry of the Arches cluster were obtained from observations with the infrared channel of the WFC3 on the HST for four epochs between 2010 and 2016. The 2010 epoch contains images in the F127M, F139M, and F153M filters (GO-11671, PI: A. M. Ghez), while the 2011, 2012, and 2016 epochs have images only in the F153M filter (GO-12318 and GO-12667, PI: A. M. Ghez; GO-14613, PI: J. R. Lu). A detailed description of the 2010–2012 observations is provided in Hosek et al. (2015, hereafter H15). The 2016 observations were designed to mimic the earlier F153M epochs in order to maximize the astrometric precision between the data sets. These observations have a field of view (FOV) of 120'' × 120'', providing coverage of at least 30% of the cluster area within successive circular annuli of width 0.25 pc out to 3 pc (Figure 1).

Zoom In Zoom Out Reset image size

We extract high-precision astrometry and photometry using the FORTRAN codes img2xym_wfc3ir, a version of the img2xym_WFC package for WFC3-IR (Anderson & King 2006), and KS2, a generalization of the software developed for the Globular Cluster Treasury Program (Anderson et al. 2008; see also Bellini et al. 2018). A detailed description of this procedure and the analysis of the subsequent astrometric and photometric errors is provided in Appendix A of H15. In short, point-spread function (PSF) fit astrometry and photometry are extracted using a grid of spatially varying PSF models across the field. No significant differences in measurement precision were found for the 2016 epoch compared to the previous epochs, with average astrometric and photometric errors of 0.16 mas (1.3 × 10⁻³ pixel) and 0.008 mag, respectively, for the brightest nonsaturated stars (error on the mean calculated from 21 individual measurements). The photometry is calibrated to the Vega magnitude system using the improved KS2 zero points derived in Hosek et al. (2018, hereafter H18), which uses significantly more stars than the original zero-point derivation in H15.

The stellar positions in each epoch are transformed into a master astrometric reference frame using a second-order polynomial transformation in both X and Y (12 free parameters). The master frame is constructed such that there is no net motion of the cluster, as only high-probability cluster members (≥0.7) in the H15 catalog are used as reference stars. An iterative process is used to match stars, calculate initial proper motions, and rematch stars using those proper motions to identify stars across the epochs. The star matching is done by position, using a search radius of 0.5 pixel (006). Proper motions are calculated for stars detected in at least three F153M epochs using a linear fit to the X and Y positions as a function of time, weighted by their astrometric errors. The final star catalog contains ∼45,000 stars with proper-motion errors 3 times smaller than those of H15 on account of the increased time baseline, reaching a precision of ∼0.03 mas yr⁻¹ at the bright end (Figure 2).

Zoom In Zoom Out Reset image size

The K-band spectroscopy of a sample of Arches cluster members was obtained using OSIRIS (Larkin et al. 2006) with laser guide star adaptive optics (Wizinowich et al. 2006) on the Keck I telescope on 2014 May 16. The Kbb filter was used with the 010 pixel scale, which provides spectral coverage of 1.965–2.381 μm at R ∼ 3800 over a 16 × 64 FOV. A single field was observed near the core of the cluster (J2000: α = 17:45:50.7, δ = −28:49:23.4; Figure 1) at a position angle of 28° using 10 dithered exposures of 900 s for a total integration time of 9000 s. This field was chosen to maximize the number of non-Wolf-Rayet (WR) stars (F153M ≥ 14.5 mag; see Section 3.4) while avoiding the densest inner region of the cluster. The spectroscopic sample contains five stars, as described in Table 1.

Table 1.  OSIRIS Spectroscopic Sample

Namea	R.A.b	Decl.b	Spectral Typec	F127M	F153M	AKsd	Teff	log g
	(J2000)	(J2000)	(literature)	(mag)	(mag)	(mag)	(K)	(cgs)
47	17:45:50.68	−28:49:24.39	O4–5 Ia	17.01 ± 0.01	14.85 ± 0.01	2.43	${34750}_{-1500}^{+3000}$	${3.50}_{-0.15}^{+0.30}$
44	17:45:50.62	−28:49:24.77	⋯	17.18 ± 0.01	14.91 ± 0.01	2.43	${34500}_{-1500}^{+3000}$	${3.75}_{-0.25}^{+0.15}$
53	17:45:50.64	−28:49:24.14	O4–5 Ia	17.10 ± 0.01	14.91 ± 0.01	2.43	${37000}_{-2000}^{+2000}$	${3.50}_{-0.10}^{+0.30}$
55	17:45:50.73	−28:49:24.54	O5.5–6 I–III	17.09 ± 0.01	14.97 ± 0.01	2.40	${34500}_{-1500}^{+3000}$	${3.85}_{-0.15}^{+0.25}$
60	17:45:50.74	−28:49:21.08	O4–5 Ia	17.20 ± 0.01	15.03 ± 0.01	2.39	${36000}_{-1500}^{+2500}$	${3.60}_{-0.15}^{+0.20}$

Notes.

^aAs defined in the catalog from Figer et al. (2002). ^bMeasured in 2010 F153M epoch. ^cFrom Clark et al. (2018). ^dDerived using the extinction map in Section 3.2.

Download table as:  ASCIITypeset image

The OSIRIS data cubes were reduced using version 4.1.0 of the OSIRIS data reduction pipeline⁷ (ODRP; Krabbe et al. 2004). The ODRP corrects for dark current, electronic biases and crosstalk, and cosmic rays and properly extracts the wavelength-calibrated spectrum at each spaxel (spatial pixel). The science data cubes were averaged together using the "Mosaic Frames" module to create the master science data cube. One-dimensional science spectra were extracted using a 3 × 3 aperture box centered on the spaxel with the highest integrated flux for the star. This aperture size was chosen to maximize the signal-to-noise while minimizing contamination from nearby stars.

After extraction, the raw science spectra need to be corrected for contamination from sky features such as continuum, OH emission lines, and telluric absorption lines. The standard set of calibration observations (sky frames and telluric standards) were obtained at the telescope, but we found that the sky features were better corrected using the Skycorr⁸ (Noll et al. 2014) and molecfit⁹ (Kausch et al. 2015; Smette et al. 2015) software packages. Skycorr removes sky emission lines by fitting physically related OH line groups in a reference sky spectrum and scaling them to match the science spectrum (e.g., Davies 2007). The sky continuum is measured by a linear interpolation of the wavelength channels without line emission and then combined with the OH line model to produce the final sky spectrum that is subtracted from the science spectrum. In this case, a reference sky spectrum for each star is extracted using a box annulus formed by a 5 × 5 and 7 × 7 spaxel box centered on the star itself and then rescaled to science spectrum aperture size. Once Skycorr has removed the sky emission and continuum, the telluric absorption lines are modeled using molecfit, which uses a radiative transfer code and an atmospheric profile based on the date and location of the observations to predict atmospheric lines caused by molecules such as H₂O, CO₂, and CH₄. The telluric model is then divided out of the science spectrum to produce the final reduced science spectrum.

However, as discussed by Lockhart et al. (2018), OSIRIS introduces a shape to the stellar continuum due to its varying sensitivity as a function of wavelength that cannot be modeled by molecfit. This requires the extra step of creating an OSIRIS "flat" free of sky, telluric, and stellar-flux contributions. We construct this flat using the observed telluric standards, empirically subtracting the sky, and using molecfit to remove the telluric lines. In the A0 V spectrum, the only remaining feature is the Br-γ line. To remove this line, we combine the A0 V and G2 V spectra using the technique described in Do et al. (2009), replacing the A0 V spectrum between 2.155 and 2.175 μm with the spectrum of the G2 V star after it has been divided by the solar spectrum. Finally, we smooth the resulting spectrum using a median filter (kernel size = 51 pixels) to create the OSIRIS flat. The science spectra are divided by this flat and normalized to produce the final science spectra (Figure 3).

Zoom In Zoom Out Reset image size

Cluster membership probabilities are calculated using the proper motions derived in Section 2.1 and the Gaussian mixture-model technique described in H15. This approach provides the flexibility needed to fit the complex kinematics of the cluster and field populations while taking the proper-motion errors into account. To reduce outliers, an error cut of 1.42 mas yr⁻¹ (one-third of the difference between the average cluster and field population proper motions in H15) is adopted, resulting in a membership catalog of 29,895 stars. This is significantly larger than the sample analyzed in H15 (∼6000 stars) because we adopt a proper-motion error cut that is 2.2× larger, do not impose a magnitude error cut, and generally have improved proper-motion errors due to the extra epoch of data. As a result, a five-Gaussian mixture model is required to fit the cluster and field populations (Figure 4), as opposed to the four-Gaussian model used in H15. This is confirmed by the Bayesian information criterion (BIC; see Equation (20) and Section 5.2 for description), which significantly favors the five-Gaussian model.

Zoom In Zoom Out Reset image size

Individual cluster membership probabilities are calculated as

$$\begin{eqnarray}&&{P}_{\mathrm{pm}}^{i}=\displaystyle \frac{{\pi }_{c}{P}_{c}^{i}}{{\pi }_{c}{P}_{c}^{i}+{\displaystyle \sum }_{k}^{K}{\pi }_{k}{P}_{k}^{i}},\end{eqnarray} \tag{ 2 }$$

where π_c and π_k are the fraction of total stars in the cluster and kth field Gaussian, respectively, and ${P}_{c}^{i}$ and ${P}_{k}^{i}$ are the probability of the ith star being part of the cluster and kth field Gaussian, respectively. A table describing the parameters of the Gaussian mixture-model fit is provided in Appendix A.

Red clump (RC) stars are used to correct for differential extinction across the field. The intrinsic magnitudes and colors of these stars do not vary significantly with age or metallicity, making them useful "standard crayons" with which to measure extinction (Girardi 2016). While not associated with the Arches cluster itself, RC stars are numerous in the Galactic bulge and have a density distribution that is sharply peaked at the GC (Wegg & Gerhard 2013). Thus, we assume that the extinction of the RC stars is similar to that of the cluster, so an extinction map derived using RC stars can be used for cluster stars. This approach was validated in H15, who showed that an RC extinction map significantly reduced the differential extinction in proper-motion-selected Arches members.

We improve the extinction map presented in H15 by using a refined sample of RC stars identified using an unsharp-masking technique (e.g., De Marchi et al. 2016) and adopting a revised version of the optical/near-infrared extinction law presented in H18 (Appendix B). The advantage of the unsharp-masking technique is that it increases the contrast of high-density features, such as the RC population, while reducing low-frequency noise. We select RC stars using the criteria described in H18: we calculate a best-fit line to the high-density RC feature in the color magnitude diagram (CMD) after unsharp masking and identify stars within ΔF153M = 0.3 mag of the best-fit line as the RC population (see Figure 7 from H18). This width is selected to encompass the RC feature and is likely caused by the distribution of stellar distances, metallicities, and ages within the population, all of which alter their location in the CMD. In addition, we consider only stars with P_clust ≤ 0.02 in order to eliminate cluster members from the sample (which is necessary since the populations overlap in CMD space) and require a photometric error better than 0.05 mag in both the F127M and F153M filters in order to remove field interlopers that scatter into the selection space. Ultimately, 875 RC stars are used in the final extinction map.

The Arches extinction map is created using a spatial interpolation of the RC star sample with a fifth-order bivariate spline¹⁰ (Figure 5). All pixels with r_cl < 0.25 pc are removed from the map, since high stellar crowding prevents an adequate number of RC stars from being detected at these radii. Ignoring the extreme values at the edge of the field where the interpolation becomes invalid, the extinction map values are in the range 1.9 mag < A_Ks < 2.65 mag, with a median extinction of A_Ks = 2.38 mag for stars with P_pm ≥ 0.5. We will adopt this as an initial estimate for the average extinction of the cluster and include a term in the IMF analysis to capture residual differential extinction in the cluster due to errors in the extinction map (Section 4).

Zoom In Zoom Out Reset image size

Observational completeness is determined using artificial-star planting and recovery tests. We plant a total of 675,000 artificial stars and run them through the same detection pipeline as the real stars. These stars are generated in three sets. The first set contains 400,000 artificial stars with magnitudes drawn from the observed CMD, perturbed by a random amount drawn from a Gaussian distribution with a width equal to the photometric uncertainty. These stars are planted uniformly across the field. The second set contains 175,000 artificial stars that are assigned to a grid of magnitudes and colors in order to cover sparsely populated regions of the CMD (e.g., the brightest and faintest observed magnitudes) in order to improve the confidence of the completeness corrections in these regions. These stars are also given a uniform spatial distribution. The final set of 100,000 artificial stars is generated based on the brighter stars in the observed CMD (F153M ≤ 18 mag) and planted according to the radial profile of the Arches cluster from H15. This increases the confidence of the completeness correction near the cluster center, where the effects of stellar crowding are strongest.

After the artificial stars are extracted by the detection pipeline, their photometric and astrometric errors are lower than the real data errors because they do not account for PSF uncertainty. Following H15, a magnitude-dependent error term is added in quadrature to the artificial-star errors so their distribution matches those of the real-star errors. Proper motions are then calculated and photometry differentially dereddened for the artificial stars in the same manner as the real stars. To be successfully recovered, an artificial star must be detected within 0.5 mag of its planted magnitude and 0.5 pixels of its planted position in at least three of the four F153M epochs and the F127M epoch and have a proper-motion error ≤1.42 mas yr⁻¹. The resulting F127M and F153M completeness curves as a function of differentially dereddened magnitudes in different cluster radius bins (0 pc ≤ R ≤ 3 pc in steps of 0.25 pc) are shown in Figure 6.

Zoom In Zoom Out Reset image size

For the IMF analysis, we calculate the completeness for each star based on its cluster radius and position in the CMD. Within a given radius bin, the CMD is binned in steps of 0.15 mag in F153M (range: 24.5–12.3 mag) and 0.2 mag in F127M–F153M (range: 0–5 mag). The completeness in each bin is assigned to the lowest value from the F127M and F153M completeness curves at the respective F153M and F127M magnitudes at the center of the bin. At the average color of the cluster, the F153M curve sets the completeness limit.

Starting with the cluster membership catalog described in Section 3.1 (29,895 stars), we apply a series of cuts in order to produce a high-quality sample for the IMF analysis. We require the following.

• 1.  
  P_pm ≥ 0.3, in an effort to reduce the number of field stars in our sample.
• 2.  
  A minimum of 30% completeness, as determined in Section 3.3. Due to the limited HST completeness at small cluster radii, we only consider stars with r_cl > 0.25 pc. We thus achieve a depth of F153M ≤ 21.18 mag, corresponding to M ≥ ∼1.8 M_⊙. We note that the analysis is not sensitive to this choice of the completeness limit; adopting a minimum completeness of 50% does not significantly impact the results, other than changing the lower mag limit to F153M ≤ 20 mag (M ≥ ∼2.5 M_⊙).
• 3.  
  A minimum of 30% area coverage within successive circular annuli of width 0.25 pc. As discussed in H15, this is achieved for r_cl ≤ 3.0 pc.
• 4.  
  All F153M measurements for a given star to agree with its median F153M magnitude within 0.5 mag. This was found to remove situations where a faint star is misidentified as a nearby bright star.
• 5.  
  WR stars removed from our sample, given the uncertainty in their stellar models and thus stellar masses. We use the population of spectroscopically identified WR stars in the Arches cluster (Figer et al. 2002; Martins et al. 2008; Clark et al. 2018) to determine their F153M magnitudes at the average cluster extinction of A_Ks = 2.38 mag. The faintest of these stars, star B1 in Clark et al. (2018), is found to have a differentially dereddened magnitude of F153M = 14.1 mag (observed F153M = 14.01 mag; the star is less extinguished than the cluster average), so we adopt a conservative magnitude cut of F153M ≥ 14.5 mag.

Finally, a photometric color cut is used to remove obvious field contaminants from the sample. High-probability cluster members (P_pm ≥ 0.6) are corrected for differential extinction as described in Section 3.2, and a 3σ clipping algorithm is used to calculate the average F127M–F153M color and standard deviation as a function of F153M magnitude. For the entire sample, stars with differentially dereddened colors larger than 2σ to the blue or 3σ to the red of the cluster sequence are automatically assigned P_pm = 0, while all others are unchanged. This color cut is more conservative to the red in order to account for the fact that some stars may have intrinsic reddening due to circumstellar disk material due to the cluster's young age (e.g., Stolte et al. 2015).

After these cuts, we are left with a sample of 981 stars with $\sum {P}_{\mathrm{pm}}=638.0$. The CMD of this sample before and after the differential extinction correction is shown in Figure 7, and a summary of the cuts and their impact on the sample size is given in Table 2.

Zoom In Zoom Out Reset image size

Table 2.  Sample Selection

Selection Description	Criterion	Nstars	$\sum {P}_{\mathrm{pm}}$
Original sample		29,895	1290.7
		Cut from Sample
Membership	Ppm ≥ 0.3	28,237
Completeness	≥0.3	539
F153M mag diff.	≤0.5 mag	45
WR stars	F153M ≥ 14.5 mag	16
Color cut	See Section 3.4	78
Final sample		980	636.7

Download table as:  ASCIITypeset image

Despite these efforts, some field contamination inevitably remains in our sample. This is due to stars with similar proper motions and colors as the cluster, so their membership probabilities are artificially inflated. In Section 5.1, we derive revised cluster membership probabilities after the IMF analysis using the best-fit cluster and field model in order to take full advantage of the photometric information. We find that the number of cluster stars based on P_pm is ∼6% larger than the number of cluster stars based on the revised membership probabilities, and we thus conclude that the sample contains approximately this amount of field contamination.

Effective temperatures and surface gravities are derived for the spectroscopic stars by comparing the spectra to non-LTE CMFGEN model atmospheres (Hillier & Miller 1998; Hillier & Lanz 2001). Non-LTE treatment is required due to the high temperatures of the stars and the presence of significant stellar winds, as evidenced by the Br-γ emission inferred from the weak Br-γ photospheric absorption line. Uncertainties in the stellar parameters are estimated by adjusting the models until they no longer provide good fits to the main diagnostic lines. Throughout the analysis, we assume a terminal velocity (V_inf) of 2000 km s⁻¹, since this cannot be constrained from the spectra.

The best-fit model spectra are shown in Figure 8, and the corresponding T_eff and log g values are reported in Table 1. Here T_eff is constrained to within ±3000 K or better and is determined primarily from the He ii/He i line ratios, as well as the absorption component of the He i 2.113 μm line. Stars 47, 55, and 60 were recently classified as O4–5 Ia stars and star 53 as an O5.5–6 I–III star by Clark et al. (2018). Our derived temperatures are consistent with the observed T_eff–versus–spectral type relation for galactic O-type stars within the uncertainties (Martins et al. 2005). The log g values are less well constrained since they rely on the weak Br-γ lines and thus are not used in the IMF analysis.

Zoom In Zoom Out Reset image size

For a cluster model Θ, we adopt a likelihood function with four components,

$$\begin{eqnarray}&&\begin{array}{l}{ \mathcal L }({{\boldsymbol{k}}}_{\mathrm{obs}},{N}_{\mathrm{cl}},{N}_{\mathrm{WR}},{T}_{\mathrm{eff}}| {\rm{\Theta }})\,=\\ \quad p({{\boldsymbol{k}}}_{\mathrm{obs}}| {\rm{\Theta }})\cdot p({N}_{\mathrm{cl}}| {\rm{\Theta }})\cdot p({N}_{W}| {\rm{\Theta }})\cdot p(\{{T}_{\mathrm{eff}},{m}_{\mathrm{obs}}\}| {\rm{\Theta }}),\end{array}\end{eqnarray} \tag{ 4 }$$

where $p({{\boldsymbol{k}}}_{\mathrm{obs}}| {\rm{\Theta }})$ is the probability of obtaining the observed distribution of stars in CMD space, with ${{\boldsymbol{k}}}_{\mathrm{obs}}$ representing the set of observed F153M magnitudes and F127M–F153M colors; $p({N}_{\mathrm{cl}}| {\rm{\Theta }})$ is the probability of detecting the number of observed cluster stars N_cl; $p({N}_{W}| {\rm{\Theta }})$ is the probability of detecting the observed number of WR stars; and $p(\{{T}_{\mathrm{eff}},{m}_{\mathrm{obs}}\}| {\rm{\Theta }})$ is the probability of measuring the observed T_eff values for the spectroscopic stars given their F153M magnitudes m_obs.

To calculate $p({{\boldsymbol{k}}}_{\mathrm{obs}}| {\rm{\Theta }})$, we must first calculate the CMD probability distribution for the cluster model and the field. The intrinsic CMD probability distribution for cluster stars generated by the model Θ, $p{({{\boldsymbol{k}}}_{\mathrm{int}}| {\rm{\Theta }})}_{\mathrm{cl}}$, is calculated according to the procedure described in Section 4. Here ${{\boldsymbol{k}}}_{\mathrm{int}}$ is the distribution of synthetic star magnitude and colors in the model cluster. To reduce the impact of stochastic effects in the synthetic CMD, the model cluster is generated with a total mass of 5 × 10⁶ M_⊙ (∼500 times more massive than the expected mass of the Arches), regardless of the M_cl designated by the model. To calculate the observed CMD probability distribution for the model cluster, we apply the observational completeness and make the same magnitude cuts as the observed sample (Section 3.4),

$$\begin{eqnarray}&&p{({{\boldsymbol{k}}}_{\mathrm{int}}| {\rm{\Theta }})}_{\mathrm{cl},\mathrm{obs}}=\displaystyle \frac{{\displaystyle \sum }_{r=0}^{{N}_{r}}p{\left({{\boldsymbol{k}}}_{\mathrm{int},r}| {\rm{\Theta }}\right)}_{\mathrm{cl}}* C(r)}{{\displaystyle \sum }_{k=0}^{{N}_{k}}{\displaystyle \sum }_{r=0}^{{N}_{r}}p{\left({{\boldsymbol{k}}}_{\mathrm{int},r}| {\rm{\Theta }}\right)}_{\mathrm{cl}}* C(r)},\end{eqnarray} \tag{ 5 }$$

where $p{({{\boldsymbol{k}}}_{\mathrm{int},r}| {\rm{\Theta }})}_{\mathrm{cl}}$ and C(r) are the intrinsic model cluster CMD probability distribution and observational completeness at a cluster radius r, N_r is the number of radius bins, and N_k is the total number of magnitude–color bins in the CMD itself.

In addition to the synthetic cluster, we construct a CMD probability distribution for the field stars. We select all stars with P_pm ≤ 0.03; apply the same differential extinction correction, magnitude, and color cuts as the IMF analysis sample; and normalize to calculate the field CMD probability distribution $p({{\boldsymbol{k}}}_{\mathrm{obs},f})$,

$$\begin{eqnarray}&&p({{\boldsymbol{k}}}_{\mathrm{obs},f})=\displaystyle \frac{{{\boldsymbol{k}}}_{\mathrm{obs},f}}{{\displaystyle \sum }_{k=0}^{{N}_{k}}{{\boldsymbol{k}}}_{\mathrm{obs},f}},\end{eqnarray} \tag{ 6 }$$

where ${{\boldsymbol{k}}}_{\mathrm{obs},f}$ is the observed field CMD. Note that we do not apply a completeness correction, since the CMD is already "observed" and thus already inherently included, and that $p({{\boldsymbol{k}}}_{\mathrm{obs},f})$ is not dependent on the cluster model.

With the cluster and field CMD probability distributions in place, we can calculate the probability of observing the ith star given its color and magnitude ($p({k}_{\mathrm{obs},i}| {\rm{\Theta }})$). We infer that the field membership probability for a given star is P_f = 1 − P_pm. To incorporate observational error, we assume that ${k}_{i}={k}_{i}^{{\prime} }+{\epsilon }_{i}$, where _i is drawn from a normal distribution centered at zero and with standard deviation drawn from the set of observational errors σ_k,i. Thus,

$$\begin{eqnarray}\begin{array}{rcl}p({k}_{\mathrm{obs},i}| {\rm{\Theta }}) & = & {\int }_{-\infty }^{\infty }({P}_{\mathrm{pm}}\ast p{\left({{\boldsymbol{k}}}_{\mathrm{int}}| {\rm{\Theta }}\right)}_{\mathrm{cl},\mathrm{obs}}+{P}_{f}\ast p({{\boldsymbol{k}}}_{\mathrm{obs},f}))\\ & & \ast \,\displaystyle \frac{1}{\sqrt{2\pi {\sigma }_{k,i}}}{e}^{\tfrac{-{\left({k}_{i}^{{\prime} }-{k}_{i}\right)}^{2}}{2{\sigma }_{k,i}^{2}}}{{dk}}_{i}^{{\prime} }.\end{array}\end{eqnarray} \tag{ 7 }$$

The final CMD likelihood is calculated by multiplying the individual likelihoods for the observed stars together,

$$\begin{eqnarray}&&p({{\boldsymbol{k}}}_{\mathrm{obs}}| {\rm{\Theta }})=\displaystyle \prod _{i=1}^{{N}_{\mathrm{obs}}}p({k}_{\mathrm{obs},i}| {\rm{\Theta }}),\end{eqnarray} \tag{ 8 }$$

where N_obs is the number of stars in the sample.

The second component of the likelihood, $p({N}_{\mathrm{cl}}| {\rm{\Theta }})$, is calculated from the number of cluster stars we would predict to observe given the cluster model. Returning to the intrinsic synthetic cluster CMD ${{\boldsymbol{k}}}_{\mathrm{int}}$, we perturb the photometry of each star by a random amount drawn from the photometric error of the observations at its magnitude and then apply the magnitude cuts and observational completeness. Following Lu et al. (2013), we linearly scale the number of stars in the simulated cluster after it is convolved with the observational completeness (N_sim) to the cluster model mass in order to obtain the expected number of observed stars N_e,

$$\begin{eqnarray}&&{N}_{e}={N}_{\mathrm{sim}}\ast \left(\displaystyle \frac{{M}_{\mathrm{cl}}}{5\ast {10}^{6}}\right),\end{eqnarray} \tag{ 9 }$$

where M_cl is the cluster model mass. The probability of obtaining the observed number of cluster stars ${N}_{\mathrm{cl}}=\sum {P}_{\mathrm{pm}}$ is calculated from a Poisson distribution:

$$\begin{eqnarray}&&p({N}_{\mathrm{cl}}| {\rm{\Theta }})=\displaystyle \frac{{N}_{e}^{{N}_{\mathrm{cl}}}{e}^{-{N}_{e}}}{{N}_{\mathrm{cl}}!}.\end{eqnarray} \tag{ 10 }$$

The purpose of applying the observational errors to ${{\boldsymbol{k}}}_{\mathrm{int}}$ for this calculation is to account for any potential Malmquist bias that is introduced by our magnitude cuts. Note that this is not done in Equation (5) for the CMD component of the likelihood, since the observational errors are already accounted for in Equation (7).

The third component of the likelihood is based on the predicted number of WR stars in the cluster model, which serves as a constraint on the cluster age (e.g., Lu et al. 2013). The brightest stars in the inner region of the cluster (r_cl <0.75 pc) were cataloged by Figer et al. (2002), and later spectroscopic studies identified 13 WR stars among this sample (Martins et al. 2008; Clark et al. 2018). In the cluster model, we calculate the number of predicted WR stars within this radius range and, similarly scaling that number to cluster model mass, calculate the probability of obtaining the observed number of WR stars,

$$\begin{eqnarray}&&p({N}_{W}| {\rm{\Theta }})=\displaystyle \frac{{N}_{{W}_{0}}^{{N}_{W}}{e}^{-{N}_{{W}_{0}}}}{{N}_{W}!},\end{eqnarray} \tag{ 11 }$$

where N_W = 13 and is the number of WR stars in the observations within r_cl < 0.75 pc, and ${N}_{{W}_{0}}$ is the number of WR stars predicted by the scaled cluster model in that same radius range.

The final component of the likelihood comes from from the T_eff measurements from the spectroscopic sample. For each star, we calculate $\overline{{T}_{{\mathrm{eff}}_{0}}}$ and ${\sigma }_{T{\mathrm{eff}}_{0}}$, which represent the median T_eff and its standard deviation for all stars in the cluster model with $({m}_{\mathrm{obs}}-{\sigma }_{{m}_{\mathrm{obs}}})\leqslant m\leqslant ({m}_{\mathrm{obs}}+{\sigma }_{{m}_{\mathrm{obs}}})$ and $({\mathrm{col}}_{\mathrm{obs}}-{\sigma }_{{\mathrm{col}}_{\mathrm{obs}}})\,\leqslant \mathrm{col}\leqslant ({\mathrm{col}}_{\mathrm{obs}}+{\sigma }_{{\mathrm{col}}_{\mathrm{obs}}})$, where m_obs, ${\sigma }_{{m}_{\mathrm{obs}}}$, col_obs, and ${\sigma }_{{\mathrm{col}}_{\mathrm{obs}}}$ are the F153M magnitude and F127M–F153M color of the observed star and its respective errors. The likelihood of measuring T_eff for the star is then

$$\begin{eqnarray}&&p({T}_{\mathrm{eff}},{m}_{\mathrm{obs}}| {\rm{\Theta }})=\displaystyle \frac{1}{{\sigma }_{\mathrm{tot}}\sqrt{2\pi }}\ast {e}^{-{\left({T}_{\mathrm{eff}}-{T}_{{\mathrm{eff}}_{0}}\right)}^{2}/(2{\sigma }_{\mathrm{tot}}^{2})},\end{eqnarray} \tag{ 12 }$$

where T_eff and ${\sigma }_{{T}_{\mathrm{eff}}}$ are the measured effective temperature and associated error of the star and ${\sigma }_{\mathrm{tot}}=\sqrt{{\sigma }_{{T}_{\mathrm{eff}}}^{2}+{\sigma }_{T{\mathrm{eff}}_{0}}^{2}}$. The likelihood of the spectroscopic sample is calculated by multiplying the individual likelihoods together,

$$\begin{eqnarray}&&p(\{{T}_{\mathrm{eff}},{m}_{\mathrm{obs}}\}| {\rm{\Theta }})=\displaystyle \prod _{i=1}^{{N}_{\mathrm{spec}}}p({T}_{{\mathrm{eff}}_{i}},{m}_{{\mathrm{obs}}_{i}}| {\rm{\Theta }}),\end{eqnarray} \tag{ 13 }$$

where N_spec is the number of stars in the spectroscopic sample.

We derive the best-fit cluster model using the Bayes theorem,

$$\begin{eqnarray}&&P({\rm{\Theta }}| {{\boldsymbol{k}}}_{\mathrm{obs}},{N}_{\mathrm{cl}},{N}_{\mathrm{WR}},{T}_{\mathrm{eff}})=\displaystyle \frac{{ \mathcal L }({{\boldsymbol{k}}}_{\mathrm{obs}},{N}_{\mathrm{cl}},{N}_{\mathrm{WR}},{T}_{\mathrm{eff}}| {\rm{\Theta }})P({\rm{\Theta }})}{P({{\boldsymbol{k}}}_{\mathrm{obs}},{N}_{\mathrm{cl}},{N}_{\mathrm{WR}},{T}_{\mathrm{eff}})},\end{eqnarray} \tag{ 14 }$$

where $P({\rm{\Theta }}| {{\boldsymbol{k}}}_{\mathrm{obs}},{N}_{\mathrm{cl}},{N}_{\mathrm{WR}},{T}_{\mathrm{eff}})$ is the posterior probability for the given model Θ, ${ \mathcal L }({{\boldsymbol{k}}}_{\mathrm{obs}},{N}_{\mathrm{cl}},{N}_{\mathrm{WR}},{T}_{\mathrm{eff}}| {\rm{\Theta }})$ is the likelihood equation, P(Θ) is the priors on the model free parameters, and $P({{\boldsymbol{k}}}_{\mathrm{obs}},{N}_{\mathrm{cl}},{N}_{\mathrm{WR}},{T}_{\mathrm{eff}})$ is the sample evidence. To sample the parameter space to find the best-fit model, we use Multinest, a publicly available multimodal sampling algorithm shown to be more efficient that Markov chain Monte Carlo algorithms when exploring complex parameter spaces (Feroz et al. 2009). We adopt an evidence tolerance of 0.5, a sampling efficiency of 0.8, and 1000 live points to run the analysis. The algorithm is run using the python wrapper module PyMultinest (Buchner et al. 2014).

We test the accuracy of this procedure by running the analysis on simulated clusters of known properties. A discussion of how the simulated clusters are created and the results of the tests are provided in Appendix D. We find that the analysis is able to recover the input values to within 1σ for all parameters for both the one-segment and two-segment IMF models.

After the best-fit cluster model is determined, we calculate revised cluster membership probabilities that take full advantage of the available kinematic and photometric information. The cluster model provides the distribution of cluster stars in CMD space, from which stars with proper motions similar to the cluster but with photometry similar to the field can be deweighted. First, we calculate the expected cluster CMD ${{\boldsymbol{k}}}_{{\rm{\Theta }},{cl}}$ and field star CMD ${{\boldsymbol{k}}}_{{\rm{f}}}$,

$$\begin{eqnarray}\begin{array}{rcl}{{\boldsymbol{k}}}_{{\rm{\Theta }},\mathrm{cl}} & = & \displaystyle \sum _{i=0}^{{N}_{\mathrm{obs}}}{P}_{\mathrm{pm},i}\ast p{\left({{\boldsymbol{k}}}_{\mathrm{int}}| {\rm{\Theta }}\right)}_{\mathrm{cl},\mathrm{obs}}\\ {{\boldsymbol{k}}}_{{\rm{f}}} & = & \displaystyle \sum _{i=0}^{{N}_{\mathrm{obs}}}{P}_{f,i}\ast p({{\boldsymbol{k}}}_{\mathrm{obs},f}),\end{array}\end{eqnarray} \tag{ 15 }$$

where $p{({{\boldsymbol{k}}}_{\mathrm{int}}| {\rm{\Theta }})}_{\mathrm{cl},\mathrm{obs}}$ and $p({{\boldsymbol{k}}}_{\mathrm{obs},f})$ are as defined in Equations (5) and (6). The revised membership probability for a given star then becomes

$$\begin{eqnarray}&&{P}_{\mathrm{clust},i}={\int }_{-\infty }^{\infty }\left(\displaystyle \frac{{{\boldsymbol{k}}}_{{\rm{\Theta }},\mathrm{cl}}}{({{\boldsymbol{k}}}_{{\rm{\Theta }},\mathrm{cl}}+{{\boldsymbol{k}}}_{{\rm{f}}})}\right)\ast \displaystyle \frac{1}{\sqrt{2\pi {\sigma }_{k,i}}}{e}^{\tfrac{-{\left({k}_{i}^{{\prime} }-{k}_{i}\right)}^{2}}{2{\sigma }_{k,i}^{2}}}{{dk}}_{i}^{{\prime} },\end{eqnarray} \tag{ 16 }$$

and P_clust is thus a combination of the proper-motion membership (which sets the scale of cluster and field CMD components) and the cluster and field CMDs themselves.

We also use the best-fit cluster model to infer the intrinsic properties (e.g., mass) for each star in the observed sample. These values are often estimated by tracing the star to a theoretical cluster isochrone along the reddening vector, but this approach is challenging near the pre-main-sequence turn-on where multiple intersections between the reddening vector and isochrone can occur. Instead, we calculate a probability distribution for the desired stellar property from ${{\boldsymbol{k}}}_{\mathrm{int}}$, based on the stars located at the observed star's location in the CMD. For example, the mass probability distribution within a given CMD bin k is

$$\begin{eqnarray}&&p{\left(m| {\rm{\Theta }}\right)}_{k}=\displaystyle \frac{{\displaystyle \sum }_{i}^{{N}_{i}}{m}_{i,b,k}}{{\displaystyle \sum }_{b}^{{N}_{b}}{\displaystyle \sum }_{i}^{{N}_{i}}{m}_{i,b,k}},\end{eqnarray} \tag{ 17 }$$

where m_i,b,k is the mass of the ith star in mass bin b in the CMD bin k, N_i is the number of stars in mass bin b, and N_b is the total number of mass bins. The mass bins are chosen to be 20 equal log-spaced bins between 0.8 and 70 M_⊙, which are the minimum and maximum masses in the cluster model.¹¹

For a given star, we calculate its mass probability distribution by multiplying $p{(m| {\rm{\Theta }})}_{k}$ by the position of the star in the CMD convolved with its photometric error:

$$\begin{eqnarray}&&\phi {\left(m\right)}_{i}={\int }_{-\infty }^{\infty }p{\left(m| {\rm{\Theta }}\right)}_{k}\ast \displaystyle \frac{1}{\sqrt{2\pi {\sigma }_{k,i}}}{e}^{\tfrac{-{\left({k}_{i}^{{\prime} }-{k}_{i}\right)}^{2}}{2{\sigma }_{k,i}^{2}}}{{dk}}_{i}^{{\prime} }.\end{eqnarray} \tag{ 18 }$$

We construct the observed IMF Φ_obs by summing the mass probability distributions over the sample, taking into account each star's revised cluster membership probability, observational completeness, and area completeness,

$$\begin{eqnarray}&&{{\rm{\Phi }}}_{\mathrm{obs}}=\displaystyle \sum _{i}^{{N}_{i}}\phi {\left(m\right)}_{i}\ast \displaystyle \frac{{P}_{\mathrm{clust},i}}{C(r)\ast a(r)},\end{eqnarray} \tag{ 19 }$$

where C(k, r) is the completeness as a function of CMD position and radius and a(r) is the area completeness. We reiterate that Φ_obs is dependent on the synthetic cluster and calculated after the best-fit model is found. It thus serves as a check that the IMF derived in the analysis is indeed a good match to the observations.

The best-fit cluster models for each of the different cases examined in this analysis (one-segment versus two-segment IMF, Pisa/Geneva versus MIST evolution models, with versus without multiplicity) are given in Table 4, and a breakdown of the corresponding likelihoods is given in Table 5. A detailed comparison of these cases is presented in Sections 5.2 and 5.3, but in summary: (1) the one-segment IMF model is slightly favored over the two-segment IMF model, but we cannot rule out the two-segment IMF model; (2) we cannot distinguish between the Pisa/Geneva and MIST evolution models in the one-segment IMF case, but the MIST models are favored in the two-segment IMF case; and (3) the fits without multiplicity are strongly disfavored. As a result, we adopt the one-segment IMF fit with the Pisa/Geneva models and multiplicity as the best-fit IMF for the Arches cluster, and we use the MIST model solution to estimate the systematic error. When discussing the two-segment IMF fit, we adopt the MIST model solution with multiplicity and use the Pisa/Geneva model solution to estimate the systematic error.

Table 4.  Best-fit Cluster Models

	One-segment IMF		Two-segment IMF
Parametera	Pisa/Genevab	MIST v1.0c	Pisa/Genevab	MIST v1.0c
α1	1.80 ± 0.05	1.68 ± 0.05	2.12 ± 0.11	${2.04}_{-0.19}^{+0.14}$
α2	⋯	⋯	0.95 ± 0.45	${1.10}_{-0.31}^{+0.39}$
mbreak	⋯	⋯	${5.4}_{-0.8}^{+2.4}$	${5.8}_{-1.2}^{+3.2}$
Mcl	${24400}_{-1600}^{+2000}$	${28600}_{-2800}^{+3000}$	${19600}_{-1600}^{+2000}$	${21000}_{-2800}^{+3400}$
log t	6.57 ± 0.02	6.56 ± 0.02	6.60 ± 0.05	${6.55}_{-0.04}^{+0.02}$
d	7900 ± 158	7900 ± 160	8030 ± 160	8100 ± 160
AKs	2.44 ± 0.01	2.44 ± 0.01	2.46 ± 0.02	2.45 ± 0.01
ΔAKs	0.15 ± 0.01	0.15 ± 0.01	0.13 ± 0.01	0.15 ± 0.01

Notes.

^aPriors and units are the same as described in Table 3. Note that only statistical uncertainties are reported in this table. Systematic uncertainties are estimated to be half the difference between the parameter values derived using different stellar evolution models. ^bPisa: Tognelli et al. (2011); Geneva: Ekström et al. (2012). ^cChoi (2016), Dotter (2016).

Download table as:  ASCIITypeset image

Table 5.  IMF Model Likelihoods

	One-segment IMF			Two-segment IMF
Component	Pisa/Genevaa	MIST v1.0b	No Multiples	Pisa/Genevaa	MIST v1.0b	No Multiples
CMD	−5058.6	−5060.2	−5067.0	−5055.8	−5051.9	−5057.4
Nstars	−4.48	−4.15	−4.23	−4.23	−4.22	−4.18
NWR	−3.24	−2.46	−4.23	−2.21	−3.45	−2.26
Spectroscopy	−19.48	−19.45	−21.08	−19.48	−19.54	−19.35
log(${ \mathcal L }$)	−5103.1	−5103.5	−5116.9	−5101.4	−5098.6	−5102.9
BIC	10247.5	10248.3	10275.1	10257.9	10252.3	10260.9

Notes.

^aPisa: Tognelli et al. (2011); Geneva: Ekström et al. (2012). ^bChoi (2016), Dotter (2016).

Download table as:  ASCIITypeset image

The posteriors for the IMF model parameters are provided in Appendix C. A comparison between the observed and model CMD is shown in Figure 9, and the subsequent F153M luminosity function shown in Figure 10. Good agreement is generally found between the observations and model, though perhaps with a slight excess of model stars at the bright end of the sample (F153M ≲ 16 mag). The agreement between the spectroscopic T_eff measurements and those predicted by the model is shown in a Hertzsprung–Russell (HR) diagram, where the (model-dependent) luminosity for each of the observed stars has been derived in the manner described in Section 4.2 (Figure 11). The luminosities $(\mathrm{log}L/{L}_{\odot }\sim 5.0\mbox{--}5.2)$ are noticeably smaller than what has been measured for O-type supergiants of similar spectral type $(\mathrm{log}L/{L}_{\odot }\sim 5.6\mbox{--}5.95$; Najarro et al. 2011; Bouret et al. 2012), though further work is required to determine if these stars are truly anomalous. The total number of cluster stars predicted by the model (618.9 ± 33) is in good agreement with the observed value $(\sum {P}_{\mathrm{pm}}=636.7)$, though the expected number of WR stars is ∼1.3σ higher than observed (18.4 ± 1.75, compared to N_WR = 13).

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

We obtain a high-mass power-law slope of α = 1.80 ± 0.05, which is either ∼1.6σ or ∼10σ lower than the local IMF value, depending on whether the uncertainty on the local IMF is considered. Perhaps a more informative comparison is that our result is ∼8.3σ lower than the measured IMF of young clusters in M31 ($\alpha ={2.45}_{-0.06}^{+0.03};$ Weisz et al. 2015). A comparison of these values is shown in Figure 12. This suggests that the Arches has a top-heavy IMF with an overabundance of high-mass stars relative to low-mass stars for M > ∼1.8 M_⊙. The α we derive does somewhat depend on which stellar evolution model we adopt, as the best-fit cluster with the MIST models has α = 1.68 ± 0.05. We thus add a systematic error term of 0.06 to our α measurement (the difference between the α values of the two fits divided by 2), so the final constraint becomes α = 1.80 ± 0.05 (stat) ± 0.06 (sys). Note that when the statistical and systematic errors are added in quadrature, the Arches result remains 6.6σ lower than the M31 result.

Zoom In Zoom Out Reset image size

The best-fit cluster age is log t = 6.57 ± 0.02 (∼3.7 ± 0.17 Myr), with negligible systematic error. This is generally older than previous ages reported in the literature. Past estimates come primarily from spectroscopic studies of the massive stars, with values of 2–2.5 Myr based on the observed nitrogen abundances and luminosities of WR stars (Najarro et al. 2004), 2–4 Myr based on the locations of WR + O stars on the HR diagram (Martins et al. 2008), 2–3.3 Myr based on the spectral types of candidate main-sequence stars (Clark et al. 2018), and ${2.6}_{-0.2}^{+0.4}\,\mathrm{Myr}$ based on the properties of an eclipsing binary in the cluster (Lohr et al. 2018). However, a cluster age of 3.7 ± 0.7 was obtained by Schneider et al. (2014) based on the shape of the PDMF relative to stellar population models with binary star evolution, which is more consistent with our result.

We infer a cluster mass of ${M}_{\mathrm{cl}}=({2.44}_{-0.16}^{+0.20}\pm 0.21)\,\times {10}^{4}\,{M}_{\odot }$, which represents the intrinsic mass between 0.8 and 150 M_⊙ out to a cluster radius of 3 pc. This assumes that the one-segment IMF model is valid over the entire mass range and that the radial profile is adequately modeled for r < 0.25 pc, which is outside the boundary of the observed sample (Section 6.4). However, the advantage of this result is that it is jointly constrained with the IMF, while previous photometric mass estimates of the cluster are needed to adopt an IMF and extrapolate it to achieve a similar depth (e.g., Serabyn et al. 1998; Figer et al. 1999; Espinoza et al. 2009).

As a consistency check, we compare the best-fit cluster model to dynamical mass estimates of the cluster made by Clarkson et al. (2012). Using the velocity dispersion of the cluster core region, they estimate the dynamical mass of the cluster to be ${0.9}_{-0.35}^{+0.40}\times {10}^{4}\,{M}_{\odot }$ for r_cl < 0.4 pc and ${1.5}_{-0.60}^{+0.74}\times {10}^{4}\,{M}_{\odot }$ for r_cl < 1.0 pc. Since our model only contains stellar masses down to 0.8 M_⊙, we would expect the enclosed mass at these radii to be lower than the dynamical estimate. This is indeed the case, with model enclosed masses of (0.74 ± 0.10) × 10⁴ and (1.2 ± 0.14) × 10⁴ M_⊙ for r_cl <0.4 and 1.0 pc, respectively.

We use the procedure outlined in Section 4.2 to calculate revised membership probabilities and Φ_obs. Figure 13 shows the P_pm and P_clust for the individual stars in the CMD. A comparison of the panels reveals the regions where P_pm > P_clust, suggesting that P_pm is overestimated due to field contamination, which is especially evident near the RC (the diagonal distribution of stars to the red of the cluster sequence at F153M ∼ 18 mag) and faint field star distribution (the stars to the blue of the cluster sequence at F153M ≥ 20 mag). The total number of cluster stars based on P_clust is 601.3, which is ∼6% smaller than what is calculated from P_pm. Thus, we estimate that P_pm (which was used in the IMF analysis) contains ∼6% field contamination.

Zoom In Zoom Out Reset image size

The observed IMF Φ_obs is shown in Figure 14. Also plotted is the Φ_obs we would obtain if we adopted a cluster model identical to the best fit but with the local IMF. The mass function obtained with the local IMF is significantly inconsistent with the observations, while the mass function obtained from the best-fit model is a good match to the observations.

Zoom In Zoom Out Reset image size

A catalog of the observed stars with membership probabilities and mass estimates is provided as a machine-readable table with this paper. A sample of the catalog is shown in Table 6.

Table 6.  Stellar Parameters

Namea	F127Mb	F153Mb	Xc	σx	Yc	σy	μαcosδ	${\sigma }_{{\mu }_{\alpha {\cos }\delta }}$	μδ	${\sigma }_{{\mu }_{\delta }}$	M	σM	AKs	Ppm	Pclust
	mag	mag	''	''	''	''	mas yr−1	mas yr−1	mas yr−1	mas yr−1	M⊙	M⊙	mag
Ae035_001	15.99	14.12	27.5032	9.50E-04	21.3643	1.15E-03	−0.03	0.025	−0.19	0.023	26.3	3.9	2.24	0.84	0.91
Aw061_002	16.38	14.47	−55.0359	1.63E-03	26.2556	1.46E-03	−0.04	0.054	−0.02	0.064	22.7	3.4	2.35	0.90	0.89
Aw048_004	16.39	14.49	−43.1954	1.12E-03	21.9127	1.22E-03	−0.22	0.013	−0.15	0.031	22.7	5.5	2.13	0.74	0.89
As017_001	16.95	14.71	7.6862	1.03E-03	−15.6326	1.15E-03	0.08	0.019	0.15	0.064	26.3	3.9	2.05	0.82	0.94
An022_002	17.27	14.72	−10.0457	1.03E-03	19.4646	1.15E-03	−0.06	0.087	−0.15	0.027	35.4	8.5	2.28	0.86	0.38
Aw006_001	17.28	14.86	−5.6225	1.03E-03	−0.0122	1.15E-03	−0.27	0.014	0.06	0.017	30.5	7.4	2.34	0.69	1.00
Ae010_001	17.04	14.92	7.3893	1.22E-03	6.8467	1.22E-03	0.15	0.040	−0.04	0.027	26.3	3.9	2.46	0.85	0.98

Notes.

^aNaming convention is as follows: The first letter is always "A", followed by "n/s/e/w" to designate whether the star is to the north, south, east or west quadrant relative to the central reference star, as determined using 45° and −45° line boundaries that intersect at the reference star position. The number immediately following gives the radius of the given star relative to the reference in arcseconds, while the second number (following the "_" separator) is a unique index for the star relative to all others at that same radius. ^bObserved magnitudes not corrected for differential extinction. ^cPositions are reported as offsets relative to a central reference star, chosen to be star 8 in the Clarkson et al. (2012) catalog. Note that this star is in the cluster core, which is excluded from this catalog. The x-direction is ${\rm{\Delta }}\alpha \cos {\rm{\Delta }}\delta$ while the y-direction is Δδ, where Δα and Δδ are the offsets in RA and DEC, respectively. Reported positional uncertainties are the uncertainties in the star position and the reference star position added in quadrature.

Only a portion of this table is shown here to demonstrate its form and content. A machine-readable version of the full table is available.

Download table as:  DataTypeset image

The best-fit two-segment cluster model is also significantly different than the local IMF, but in a different manner than the one-segment IMF model. While the high-mass IMF slope is perhaps slightly shallow $({\alpha }_{1}={2.04}_{-0.19}^{+0.14}\pm 0.04)$, the real discrepancy is in the detection of a significant m_break at ${5.8}_{-1.2}^{+3.2}\pm 0.02\,{M}_{\odot }$, which is an order of magnitude larger than the local IMF (m_break = 0.5 M_⊙). The power-law slope below m_break is ${\alpha }_{2}={1.10}_{-0.31}^{+0.39}\pm 0.08$, which is consistent with the local IMF values of 1.3 ± 0.3 for 0.08 M_⊙ ≤ M < 0.5 M_⊙ (Kroupa 2002). As a result of the high m_break, the two-segment IMF solution could be characterized as "bottom-light," with a deficit of low-mass stars relative to the local IMF. Figure 15 shows the two-segment model compared to the observed luminosity function and the derived Φ_obs.

Zoom In Zoom Out Reset image size

One of the advantages of the Bayesian framework is that we can distinguish between one- and two-segment IMF models by comparing the likelihoods of the best-fit solutions. We use the BIC (Schwarz 1978) for this comparison,

$$\begin{eqnarray}&&\mathrm{BIC}=\mathrm{ln}(n)\ast k-2\ast \mathrm{ln}({ \mathcal L }),\end{eqnarray} \tag{ 20 }$$

where n is the total number of stars in the sample (980 in this case), k is the number of free parameters in the model (i.e., six for the one-segment model and eight for the two-segment model), and ${ \mathcal L }$ is the best-fit likelihood of the model. When comparing the two models, the model with the lowest BIC is preferred, and the absolute value of the difference between the BIC values (ΔBIC) indicates if the preference is statistically significant. Table 5 contains the likelihoods and BIC values for the one- and two-segment IMF fits.

The one-segment IMF model is slightly preferred over the two-segment IMF model in both the Pisa/Geneva and MIST cases, with ΔBIC = 10.4 and 8.2, respectively. To assess the significance of these values, we generate artificial clusters with one- and two-segment IMFs as described in Appendix D (adopting the best-fit values in the Arches solutions), fit them in both the one- and two-segment cases, and calculate the corresponding ΔBIC values. In our simulations, the ΔBIC values between the one- and two-segment models are consistently a factor of 1.5–7.5 times greater the actual Arches data. Thus, the real data are substantially more agnostic on the distinction between one- and two-component IMF models than the synthetic data sets. We conclude that we cannot rule out the two-segment IMF model, though we adopt the one-segment IMF model as the overall best fit. In either case, our results show that the Arches cluster IMF is significantly different from the local IMF.

Table 4 reveals that the best-fit model parameters are only weakly dependent on the choice of stellar evolution model. Similar to Section 5.2, we use the BIC test to determine whether our analysis prefers one set of evolution models over the other. For the one-segment IMF model, the Pisa/Geneva model solution is slightly favored with ΔBIC = 0.8. However, artificial-cluster tests show that this ΔBIC is not significant, as the average difference between evolution model fits is ΔBIC = 2.1 ± 1.2. Thus, we conclude that we cannot distinguish between the two solutions, and we adopt the Pisa/Geneva solution as the result and use the MIST solution to estimate the systematic error. In the two-segment IMF case, the MIST solution is favored with a ΔBIC = 5.6. The simulations show that this discrepancy is indeed significant, with an average difference of ΔBIC = 5.39 ± 2.56 between two-segment IMF fits using different evolution models. As a result, we adopt the MIST solution for the two-segment IMF case and use the Pisa/Geneva solution to estimate the systematic error.

Whether stellar multiplicity is accounted for in the cluster model is found to significantly impact the quality of the fit. The BIC analysis strongly favors the models that include stellar multiplicity, with ΔBIC values of 27.6 and 8.6 for the one- and two-segment IMF model cases, respectively. As seen in Table 5, this difference is primarily driven by the CMD likelihood component. Artificial-cluster tests show that the observed ΔBIC values are significant; for artificial clusters that have intrinsic multiplicity, ΔBIC = 12.2 ± 0.5 in favor of the fit with multiplicity in the one-segment IMF case and ΔBIC = 8.8 ± 0.5 in the two-segment IMF case. Thus, we adopt the model fits with multiplicity included over those without.

Our result that the high-mass slope of the Arches IMF is significantly top-heavy differs from previous photometric studies of the cluster that have found the IMF to be largely consistent with the local IMF (Kim et al. 2006; Espinoza et al. 2009; Habibi et al. 2013; Shin & Kim 2015). However, a key advantage of this study is the use of proper motions to calculate cluster membership probabilities, which produces a significantly more accurate sample of cluster members than is possible through photometry alone. For example, Figure 16 shows a comparison between cluster samples obtained using proper motions versus a photometric color cut similar to Habibi et al. (2013). Even when limited to r < 1.5 pc and M > 10 M_⊙ (the range of PDMF was measured by Habibi et al. 2013), the photometric sample is systematically larger than the proper-motion selection due to field contamination. On the other hand, adopting stricter color cuts on photometric samples can be problematic as well, as Espinoza et al. (2009) noted that the color cuts they adopted forced them to eliminate stars that could be high-mass (M > 16 M_⊙) cluster members.

Zoom In Zoom Out Reset image size

An alternative approach is to statistically subtract the field from the cluster using the field population observed in nearby control fields (e.g., Kim et al. 2006; Shin & Kim 2015). However, differential extinction can alter both the average extinction and the distribution of extinction values between two fields (e.g., note the detailed extinction structures in Figure 5). As a result, it is challenging to obtain a sufficiently accurate model of the field stars in the cluster observations. In addition, care must be taken that the control fields are beyond the extent of the cluster, which H15 showed extends to a radius of at least 75'' (∼3 pc).

It is interesting to note that several previous studies have reported evidence of an enhancement in the PDMF at ∼6 M_⊙, whether it be evidence of a turnover (Stolte et al. 2005) or a localized "bump" in the mass function (Kim et al. 2006). The presence of such a feature may be driving the two-segment IMF model solution. Future studies are needed to extend the proper-motion-selected sample to lower masses in order to definitively distinguish between the one- and two-segment IMF models and determine whether an enhancement at 5–6 M_⊙ truly exists.

The top-heavy IMF we obtain for the Arches cluster (α = 1.80 ± 0.05 ± 0.06) is in good agreement with the YNC (α = 1.7 ± 0.2 for M > 10 M_⊙; Lu et al. 2013). This suggests that this unusual IMF extends beyond the central parsec of the Galaxy and into the CMZ, which spans a galactocentric radius of ∼200 pc (Morris & Serabyn 1996). Unfortunately, the exact birth location of the Arches is not well constrained due to the range of possible orbits allowed by the three-dimensional motion of the cluster (Stolte et al. 2008; Kruijssen et al. 2015). Further, the proper motion of the cluster in the galactocentric reference frame is not yet well determined, as current estimates are based on the relative proper motion between the cluster and a single-Gaussian kinematic model for the field (e.g., Clarkson et al. 2012). In reality, the field exhibits a more complex kinematic structure (see H15 and Appendix A), so the cluster motion may need to be revised. This is left to a future paper.

However, this result raises the question of whether the top-heavy IMF is truly due to the GC environment or if it is a general property of YMCs (see review by Portegies Zwart et al. 2010). In Figure 17, we compare IMF measurements of YMCs in the Milky Way disk to the YNC and Arches cluster at the GC. The YMC sample includes Westerlund 1 (Wd1; Gennaro et al. 2011; Lim et al. 2013; Andersen et al. 2017), Westerlund 2 (Wd2; Zeidler et al. 2017), NGC 3603 (Harayama et al. 2008; Pang et al. 2013), Trumpler 14 and 16 (Hur et al. 2012), and h and χ Persei (Slesnick et al. 2002).

Zoom In Zoom Out Reset image size

Figure 17 shows that the YMCs in the Galactic disk are generally consistent with the local IMF, though potential discrepancies exist. In particular, NGC 3603 has been found to be potentially top-heavy ($\alpha ={1.74}_{-0.47}^{+0.62}$ and 1.88 ± 0.15; Harayama et al. 2008 and Pang et al. 2013, respectively). However, these results may be biased due to mass segregation, which both studies find to be significant in the cluster. Indeed, the uncertainty in the Harayama et al. (2008) measurement is quite large in order to account for this (as well as other) systematic uncertainties, while Pang et al. (2013) acknowledged that their IMF measurement is restricted to the inner 60'' of the cluster. The IMF of Wd1 is potentially discrepant as well, with reported high-mass IMF slopes that are both near-standard ($\alpha ={2.44}_{-0.20}^{+0.08};$ Gennaro et al. 2011, via near-infrared photometry) and top-heavy (α = 1.8 ± 0.1; Lim et al. 2013, via optical photometry). The inconsistency between these studies makes it difficult to draw conclusions about the IMF of Wd1, though the low-mass stellar content of the cluster has been found to be consistent with the local IMF (Andersen et al. 2017). These cases highlight the difficulty of these measurements, as differences in cluster membership selection, stellar models, and methodology may significantly impact results.

Given the uncertainties surrounding the NGC 3603 and Wd1 measurements, the fact that other YMCs in the Galactic disk have been found to be consistent with the local IMF while the Arches and YNC are top-heavy provides tentative evidence that the top-heavy IMF is indeed caused by the extreme GC environment. We discuss the implications of a top-heavy IMF at the GC in Section 6.3 and the caveats of our Arches IMF measurement in Section 6.4.

The Quintuplet cluster is a third YMC in the CMZ that provides an additional probe of the IMF at the GC. A proper-motion-based analysis of the Quintuplet mass function was carried out by Hußmann et al. (2012), who found a top-heavy PDMF ($\alpha ={1.68}_{-0.09}^{+0.13}$) for the inner 0.5 pc of the cluster. However, it is uncertain whether this is due to mass segregation or a top-heavy IMF. A study of the Quintuplet IMF using a similar approach as this work is in progress.

At first, a top-heavy IMF at the GC appears to favor star formation models where the increased thermal Jeans mass leads to the formation of more high-mass stars (e.g., Larson 2005; Bonnell et al. 2006; Klessen et al. 2007; Bonnell & Rice 2008; Papadopoulos et al. 2011; Narayanan & Davé 2013). However, the main prediction of these models is that the break mass of the IMF should increase, leading to a deficit of low-mass stars, rather than an overabundance of high-mass stars and a shallow high-mass slope. This behavior is similar to the "bottom-light" two-segment IMF solution, but we do not yet have enough evidence to conclude that this is preferred over the top-heavy one-segment IMF solution.

However, our results are generally inconsistent with models where the IMF is set by the CMF (e.g., Padoan & Nordlund 2002; Hopkins 2012). Though the combination of turbulence and gravity naturally produces a CMF with a shape similar to the local IMF, these models predict a steeper mass slope and a bottom-heavy IMF near the GC (Hopkins 2013; Chabrier et al. 2014). This suggests that the CMF cannot be directly mapped to the IMF and that additional physical processes are involved. On the other hand, it has been suggested that the gravoturbulent fragmentation of a molecular cloud may lead to a top-heavy IMF (and possibly a bump around 5–6 M_⊙) if the coalescence of prestellar cores is highly efficient, as might be expected in Arches-like environments (Dib et al. 2007). In addition, recent simulations have shown that a top-heavy IMF can be produced in turbulence-dominated environments if the turbulent probability density distribution can be described as a power law at high densities (Lee & Hennebelle 2018).

Previous studies have shown that radiative feedback (e.g., Bate 2009; Offner et al. 2009; Krumholz 2011), protostellar outflows (e.g., Krumholz et al. 2012; Federrath et al. 2014), and magnetic fields (e.g., Hennebelle et al. 2011; Myers et al. 2013) can impact the IMF as well. Only recently have simulations begun to incorporate all of these processes simultaneously (Myers et al. 2014; Krumholz et al. 2016; Cunningham et al. 2018; Li et al. 2018). However, these simulations have been limited to molecular clouds with initial masses ≤1000 M_⊙ and are only applicable to low-mass stars in environmental conditions similar to local star-forming regions. Future simulations of higher-mass molecular clouds in starburst-like environments are needed in order to determine what physics is behind a shallow high-mass IMF slope in the GC.

A caveat of our IMF measurement is that we do not take the potential effects of tidal stripping into account. Tidal stripping might play a significant role in the evolution of the Arches cluster given the strong Galactic tidal field near the GC. Since tidal stripping preferentially removes low-mass stars from a cluster (e.g., Kruijssen 2009; Lamers et al. 2013), it could bias the mass function to appear top-heavy. However, it is unclear from current dynamical models of the Arches whether tidal stripping would significantly impact the mass range examined in this study (M ≳ 1.8 M_⊙). The N-body simulations by Habibi et al. (2014) predict the formation of tidal tails out to 20 pc from the cluster core that potentially contribute to the population of isolated massive stars observed near the GC. Additional simulations by Park et al. (2018) also predict the formation of tidal tails but find that ∼98% of the tidally stripped stars have masses less than 2.5 M_⊙ and that the impact on the mass function above this limit is minor. This is consistent with the observations of H15 that find no evidence of tidal tails down to ∼2.5 M_⊙ and a cluster radius of 3 pc. Thus, we assume that the effects of tidal stripping can be ignored for the mass range in our sample.

It should be noted that the dynamical models discussed above require assumptions regarding the initial conditions and orbit of the Arches cluster, both of which are quite uncertain. In addition, only stars are considered in the simulations, though the expulsion of excess gas after cluster formation is expected to have a significant impact on the dynamical evolution of the cluster as well (e.g., Bastian & Goodwin 2006; Goodwin & Bastian 2006; Farias et al. 2015).

Another caveat is that this analysis does not contain data for r < 0.25 pc, where the observational completeness is low due to stellar crowding. We adopt the radial profile of Espinoza et al. (2009) for this region when modeling the cluster (Section 4), but their profile was derived only for stars with M > 10 M_⊙. So, while magnitude-dependent radial profiles are used to account for mass segregation in the range 0.25 pc < r < 3.0 pc, the profile for all stars within the cluster core is assumed to be the same. Combining the HST data set from this study with higher-resolution ground-based observations of the cluster core would remove the need for this assumption.

Finally, we note that changing the extinction law does not have a significant impact on the IMF results. To demonstrate this, we repeat the analysis using the Nishiyama et al. (2009) and original H18 extinction laws, which are shallower (i.e., lower A_λ/A_Ks values) and steeper (i.e., higher A_λ/A_Ks values) than the law we ultimately adopt, respectively (Appendix B). In both cases, the only parameter that significantly changes is the overall extinction, which decreases for the H18 law (A_Ks = 2.07 ± 0.01 mag) and increases for the Nishiyama et al. (2009) law (A_Ks = 2.47 ± 0.01 mag). The high-mass IMF slope only changes by $| {\rm{\Delta }}\alpha | =0.01$ in the one-segment case and $| {\rm{\Delta }}\alpha | =0.04$ in the two-segment case, well within the 1σ uncertainties. Additionally, $| {\rm{\Delta }}{\alpha }_{2}| =0.03$ and $| {\rm{\Delta }}{m}_{\mathrm{break}}| =0.49\,{M}_{\odot }$ for the two-segment case, again well within the uncertainties. Therefore, the extinction law is not a significant source of uncertainty in this analysis.