© ESO, 2017

1. Introduction

“The two Magellanic clouds, Nubecula major and Nebecula minor, are very remarkable objects... In no other portion of the heavens are so many nebulous and stellar masses thronged together in an equally small space”¹.

R 136 is a very massive young star cluster that lies at the center of the Tarantula nebula in the Large Magellanic Cloud (LMC). Hosting the most massive stars known in the Local Universe (Crowther et al. 2010, 2016), R 136 provides a unique opportunity to study the formation of massive stars and clusters in the early stages of evolution.

Our understanding of the true nature of R 136 has constantly improved with increasing telescope resolution. The fuzzy object in the core of the cluster was initially thought to be a super-massive star with a mass in excess of 1000 M_⊙ (Feitzinger et al. 1980; Cassinelli et al. 1981; Savage et al. 1983). Image-sharpening techniques such as speckle interferometry (Weigelt & Baier 1985) revealed that this object had many individual stars, however, which ere coined the R 136a Weigelt components. R 136a was clearly resolved into hundreds of sources by the Hubble Space Telescope (HST, Hunter et al. 1995; Campbell et al. 1992), which opened routes to better study each of Weigelt components at various wavelengths and resolutions. More recent multiconjugate adaptive optics instruments (AO, Campbell et al. 2010) on the VLT attempted to better resolve the core of R 136 with relative success.

The combination of photometry, ultraviolet spectrometry (STIS/MAMA, Crowther et al. 2016), visible (HST/FOS, Massey & Hunter 1998) and near-infrared (VLT/SINFONI, Schnurr et al. 2009) observations has resulted in more constraints on the R 136 stellar population and its most luminous stars. Still, questions persist on the true nature of these stars: their multiplicity, their mass and the age of the cluster as a whole are still unclear.

This paper presents the first observations of R 136 in the infrared by the second-generation ESO Very Large Telescope (VLT) instrument Spectro-Polarimetric High-contrast Exoplanet Research (SPHERE²; Beuzit et al. 2008) aiming at uncrowding the dense central region of R 136 in the near-infrared (near-IR). Thanks to SPHERE’s extreme AO system, we reached the same resolution as HST in the visible with the VLT 8.2 m Melipal telescope in the Ks band and have surpassed it in J band thanks to the angular resolution, namely 0.035–0.055 arcsec in near-IR and with a better pixel sampling (12.25 mas/pix).

2. Observations

We collected data in guaranteed time observation (GTO) runs to image R 136 using the classical imaging mode of IRDIS (Langlois et al. 2014). For our purpose we use the same spectral band split into two channels to correct for residual detector hot pixels and uncorrelated detector noise among other instrumental effects. With this method we maximized the image sharpness of the total exposure by discarding the single frames with poorer Strehl ratio and corrected the residual tip-tilt error on each short exposure through software. Observations were performed in September 2015, with high dynamic and high angular resolution imaging in J and Ks bands, over a field of view (FoV) of 10.9″ × 12.3″, centered on the core of the cluster (Fig. 1 bottom). The seeing was 0.63 ± 0.1″ during the observations. The night was rated as clear. Less than 10% of the sky (above 30 degree elevation) were covered in clouds, and the transparency variations were lower than 10%.

In order to qualify our data, we compared the reduced J- and Ks-band images with published images of R 136 from HST (WFPC2 and WFC3) in V band and the VLT/MAD imaging (Campbell et al. 2010) in K band (Fig. 1 top). This comparison confirms that our data present better spatial resolution and point-spread function (PSF) sampling more suitable for applying deconvolution techniques.

Fig. 1 Comparison of R 136 core images at different wavelengths with the highest available angular resolution telescopes. The FoV of all images is the same as IRDIS data (10.9″ × 12.3″). Top left: HST/WFPC2 in V band (λcen 526 nm), top right: VLT/MAD in Ks band (λcen 2200 nm), bottom right: SPHERE/IRDIS/Ks (λcen 2182 nm); bottom left: SPHERE/IRDIS/J (λcen 1245 nm).

Our data consist of 300 frames of 4.0 s exposures in the two IRDIS broad-band Ks and J filters (BB-Ks, BB-J). The Wolf-Rayet star R 136a1 was used to guide the AO loop of SPHERE, confirming its better-than-nominal performances that surpass NACO and MAD observations.

The range of airmass during these observations was 1.54 to 1.67. A log of observations is presented in Table 1. We used the SPHERE pipeline package³ to correct for dark (current), flat (fielding), (spatial) distortion, bad pixels and thermal background due to instrument and sky (in Ks). Furthermore and in order to reach the highest sensitivity and the largest number of detectable sources, additional corrections were carried out on the images. Based on a Gaussian fit using selected stars, we estimated and corrected the subpixel image drifts before combining the individual images. This allowed us to correct for residual tip-tilt errors with an accuracy of a few mass. Still, some uncorrected atmospheric leaks persist in our final images; they are due to the adaptive optics residual halo, which is stronger in J than in Ks. We accurately considered this to estimate correct error bars in addition to the Starfinder reduction tool, which provided photometric errors (see Sect. 3).

3. Photometry of the R 136 compact core

We used the Starfinder package implemented in IDL (Diolaiti et al. 2000). Starfinder is designed for the analysis of AO images of crowded fields such as the Galactic Center (Pugliese et al. 2002). It determines the empirical local PSF from several isolated sources in the image and uses this PSF to extract other stellar sources across the FoV. Starfinder also estimates the formal error on the estimated photometry based on photon noise, variance due to the sky, and the PSF fitting procedure itself. This error is called “PSF fitting error” hereafter.

Our photometry analysis of R 136 was conducted in two steps: 1) stellar source detection using Starfinder and 2) the background analysis to obtain realistic error bars on the photometry of individual stars beyond the formal PSF fitting error.

In the first step, three well-isolated (within 0.47″ × 0.47″) stars were used to extract the PSF in J- and Ks-band data separately. The extracted PSFs from J- and Ks-band data were used as input for stellar source detection by Starfinder. The FWHM of the extracted PSFs are 54.71 and 65.16 mas in J- and Ks-band data, respectively. Consequently, 1110 and 1059 sources were detected in the J and Ks bands, respectively. We stopped source extraction after we attained a minimum correlation coefficient of 65% and 80% in Ks and J bands between the extracted star with the locally determined PSF according to Starfinder procedures. Stars with higher correlation coefficients, that is, with more similarity to the PSF, represent a higher reliability of their photometric measure. Figure 2 shows the map of the correlation coefficient across the field of J (left) and Ks (right) images of R 136. The PSF changes as a function of radial distance and azimuth along the field (Fig. 2). The AO correcting efficiency degrades as a function of distance from R 136a1, which is the reference star for the AO loop. At the borders of the FoV, the isoplanatic limits are approached. Overall, the PSF is not centro-symmetric at large distances from R 136a1. We took these distortions into account to estimate the local statistical errors, which become significant on the distant sources from the center of the image, typically >3′′.

In addition to the correlation coefficient criterion, we applied the limit of standard deviation from the sky brightness (σ_sky) to stop the extraction of sources by Starfinder, that is, the local PSF maximum value must exceed 2σ_sky over the sky. The faintest common detected stars between J and Ks band have a signal-to-noise ratio (S/N) better than 2.

The Strehl ratio (SR) in J and Ks is determined as 0.40 ± 0.05 and 0.75 ± 0.03, respectively, on average. These estimates are assessed from the SPARTA files recorded during simultaneous SPHERE runs with AO-corrected images of R 136. They are based on the slope measurements of the Shack-Hartmann wavefront sensor in SAXO by extrapolating the phase variance deduced from the reconstruction of SAXO open-loop data using deformable mirror, tip-tilt voltages, and wavefront sensor closed-loop measures (Fusco et al. 2004). The method has been proved to be quite robust in the past for differential imaging by NACO (Maire et al. 2014). We are therefore quite confident of our photometry corrected for the SR effect in J and Ks bands. To convert stellar fluxes into magnitudes, we set the zero-point magnitude to the instrumental zero-point (one ADU/s in J and Ks are 25.405 and 24.256 mag, respectively) in addition to the SR effect⁴.

Fig. 2 Map of the correlation coefficients between the input PSF and each detected stellar source, calculated by Starfinder across the FoV in J (top) and in K (bottom) images of R 136 by SPHERE.

In the second step, after extracting sources, the background image was used to estimate the residual errors in addition to the formal photometric PSF fitting errors of Starfinder. The background image contains 1) the AO halo from the atmospheric turbulence, 2) residuals from the photometric analysis in the first step and 3) undetected faint stars. We define the residual error as the fluctuation of the background caused by the remaining flux from the photometry using Starfinder and also from the undetected faint sources.

Since the core of R 136 is crowded with the brightest stars, the AO halo is brightest at the center of R 136. We removed this large halo in the Fourier space by applying a high pass filter (hat function with the diameter equal to the FWHM of the background halo), in order to estimate the fast variations of the background at the scale of the PSF. The residual error is the standard deviation of the background (after removing the stellar PSF and the AO halo) in the 30 × 30 pixel² area around the source. The final photometric error is set to the quadratic combination of the PSF fitting errors, residual errors, and the photometric calibration error (0.042 in K and 0.128 in J).

Figure 3 shows these errors separately for each detected source. The PSF fitting errors (red-pluses) are smaller in Ks band because the AO works better in longer wavelengths. The residual errors (blue crosses) are larger in Ks band, however, because the background fluctuations at longer wavelengths are higher. Figure A.1 shows the maps of these two errors in J and Ks across the field.

The error on the position of the detected stars is estimated by Starfinder and shown in Fig. 4. All stars brighter than 16.7 mag in Ks and 17 mag in J have position errors smaller than 0.1 mas (about 5 AU).

Fig. 3 PSF fitting errors (red pluses) and residual errors (blue crosses) in Ks (top) and J (bottom). The PSF fitting error is the outcome of the Starfinder. The residual errors are the outcome of the background analysis after removing the stellar sources from the image.

Fig. 4 Error in the position of the stars in Ks (top) and J (bottom) data. This error is the outcome of the Starfinder and decreases for the bright stars.

In order to interpret the photometric distribution of the 1110 and 1059 sources in J and Ks bands, we conducted an incompleteness test to the core (r < 3″) of R 136 and outside its core (r > 3″) in both J- and Ks-band imaging data. Figure 5 shows the result of incompleteness tests, which were performed by adding artificial stars with known magnitudes, one by one, to the original images. These experiments were repeated 500 times for each flux value (magnitude).

The faintest star in the artificial star test has a flux of 83.14 [ADU/s] in J band (Jmag = 21.87) and 328.614 [ADU/s] in K-band (Kmag = 19.91). For this star we iterated the test 500 times in the inner region (r < 3″) and 500 times in the outer region (r > 3″). This experiment was made for 187 (in J) and 158 (in Ks) fluxes, 500 times for each flux value, with the Δ(flux) about 12.5 [ADU/s] on average. Each time, one artificial star was added to these images and Starfinder was used again to detect the artificial star. The core of the cluster is very crowded so that the incompleteness does not reach 100% even for the bright artificial stars. This effect is more important for J-band data where the core is fuzzy because of a lower AO correction.

Fig. 5 Incompleteness test in J (red) and Ks (blue) for two regions: the very core of R 136 (r < 3″), shown with pluses and crosses, and outside of the core (r > 3″), shown with circles and triangles. For each flux value (magnitude) we used 500 artificial stars in order to find the completeness value. Note that artificial stars are added one by one to the original image, and there is always one artificial star in each experiment.

Fig. 6 Top: CMD of 818 detected sources in J- and Ks-band images of SPHERE/IRDIS from the core of R 136. Solid black, green and blue lines show the PARSEC isochrones at the ages of 0.6, 1 and 1.5 Myr (corrected for distance modulus of 18.49 and central values of extinctions, AJ = 0.45 mag and AK = 0.2 mag). The horizontal dashed gray line depicts the incompleteness of 70%. The CMD is plotted for the whole FoV (818 sources), in the very core of the cluster (r < 3″) and outside (r > 3″), from left to right respectively. The error bar on each point is the combination of the PSF fitting errors and the residual errors from the background image after removing the stellar source signals from the images. The 54 spectroscopically studied stars by Crowther et al. (2016) are shown as black dots in the middle plot. Bottom: same as top CMD, but in 3D to show the number density of stars in the CMD. Solid red, black, pink, and yellow lines depict PARSEC isochrones at the age of 0.5, 1, 1.5, and 2 Myr, respectively.

4. Age and extinction

To estimate the stellar ages and the extinction of the core of R 136, we used the effective temperature (T_eff) and luminosity (log L) of 54 stars that were studied spectroscopically by Crowther et al. (2016). We also chose a grid of isochrones at different ages (from 0.1 up to 8 Myr) with the LMC metallicity (Z = 0.006) from the latest sets of PARSEC evolutionary model⁵ (Bressan et al. 2012), which is a complete theoretical library that includes the latest set of stellar phases from pre-main sequence to main sequence and covering stellar masses from 0.1 to 350 M_⊙. Figure 7 shows these selected 54 stars with their T_eff−log L (with their error bars) and the sets of isochrones covering them.

Fig. 7 Left panel: the IRDIS/Ks (FoV 6″ × 6″) on which the 54 spectroscopically known stars from Crowther et al. (2016) have been added as red circles. The right plot depicts the Teff, log L/L⊙, and corresponding error bars on these 54 sources taken from Crowther et al. (2016) in blue. The solid red lines indicate the PARSEC isochrones covering ages from 0.1 to 8 Myr.

By fitting the isochrones to each star, we estimated the age and intrinsic color of each star with the error bars. We adopted the LMC distance modulus (DM) of 18.49 mag estimated by Pietrzyński et al. (2013), which is consistent with the value suggested by Gibson (2000) for the LMC.

Table A.1 shows the estimated age, initial and actual mass and extinction of these 54 spectroscopically known stars. The values of T_eff and log L/L_⊙ in the second and third column are taken from Crowther et al. (2016). Figure 8 shows the generalized histogram of the age of these 54 sources. We note that the age of each star has a Gaussian distribution with a given σ (error) in the histogram, derived from the errors in T_eff and L_∗. We also note that the large errors on the age and extinction come from the large spectroscopic uncertainties (errors on T_eff and log L/L_⊙). We were also limited by the evolutionary tracks up to 350 M_⊙, which explains the upper mass limit of 348 M_⊙ for very massive stars like R 136a1.

Figure 9 shows the histogram of the extinction in J and K and their color-excess.

The age of $\mathrm{1.}{\mathrm{8}}_{-0.8}^{\mathrm{+}\mathrm{1.2}}$ Myr is the most probable age range for these stars. Comparing Fig. 8 with Fig. 11 of Crowther et al. (2016), the age estimate for this sample of stars is consistent within the error bars in both analyses.

The extinction in J and Ks is 0.45 ± 0.5 and 0.2 ± 0.5 mag, respectively. The estimated extinctions are consistent with the values derived in previous studies by Campbell et al. (2010) and De Marchi & Panagia (2014). The color excess of E(J−K) = 0.25 ± 0.1 is also consistent with Tatton et al. (2013). Figure 6 shows the color–magnitude diagram (CMD) of detected sources in J- and Ks-band IRDIS data with their error bars. The CMD is plotted for the whole FoV (818 sources), in the very core of the cluster (r < 3″), and outside (r > 3″), from left to right. The error bars on each point are the combination of the PSF fitting errors and the residual errors from the background image after removing the stellar source signals from the images. The PARSEC isochrones at three different ages (0.6, 1, and 1.5 Myr) also are plotted in this figure using DM = 18.49 and central extinction values in J (0.45 mag) and K (0.2 mag).

We note that based on the CMD isochrone fitting, the age of 1 and 1.5 Myr for the inner (r < 3″) and outside core (r > 3″) fits best to our data. We can also see a young population with an age of about 0.6 Myr in both regions. In the histogram of the estimated age for the 54 spectroscopically known stars there is a peak at 0.6 Myr that shows a probable age for some stars. We can clearly see the main sequence and pre-main sequence conjunction, which is very sensitive to age. The age and extinction can be estimated using the precise photometric analysis for the 818 stellar sources in the core of R 136. Interestingly, the age of stars in the very core (r < 3″) is younger than in its outer region (r > 3″). This result appears to be consistent with the older population found beyond the core of R 136, like the northeast clump observed by Sabbi et al. (2012) and at a greater distance (3′), the old cluster, Hodge 301, studied by Grebel & Chu (2000). Understanding this apparent age trend in R 136 and also 30 Dor region is possible in the future from the formation point of view. Moreover, the question is how the younger population in the center of the cluster can be explained by star cluster formation scenarios (Wünsch et al. 2017)? We note that this age difference in the two regions can also be explained by an observational bias because the central region of R 136 is very compact and bright so that the incompleteness level is very low. Figure A.2 shows the incompleteness map at J = 18 (top) and J = 19 (bottom). The incompleteness is low in the very central regions, which means that we did not detect the true number of stars. This might affect the CMD in this magnitude range, which is sensitive to age.

Considering the errors on the extinction, we can estimate the stellar mass range for each star at a given age. We estimated the stellar masses only for common sources between J and Ks data (818 sources shown in Fig. 6) using their J and Ks magnitudes fitted to PARSEC isochrones at three different ages: 0.6, 1, and 1.5 My. The histogram of mass, which is the mass function (MF), is plotted considering a Gaussian distribution for each stellar mass. The Gaussian uncertainty in the mass of each star was accounted for when constructing the MF. Figure 10 shows the generalized histogram of the mass (MF) at three different ages (0.6, 1, and 1.5 Myr). The MF slope for 1 and 1.5 Myr isochrone is Γ_{1 Myr} = −0.90 ± 0.13 and Γ_{1.5 Myr} = −0.98 ± 0.18, respectively, for the mass range of (3−300) M_⊙. These values are low limits to steepness because of incompleteness and central concentration of bright stars. These slopes⁶ are flatter than the (Γ = −1.3, Kroupa 2001) and value (Γ = −1.35, Salpeter 1955). We note that the incompleteness correction is applied for the low-mass bin considering the number of stars in the inner (r < 3″) and outer (r > 3″) region of the cluster (see Fig. 5). The derived MF is limited to the resolution of the instrument and also to the detection limit of the observation. In future, using higher angular resolution data, we may resolve binaries and low-mass stars, which affect the shape of MF.

Fig. 8 Generalized histogram of the age of 54 known stars from Crowther et al. (2016, shown in Fig. 7). We used PARSEC isochrones at different ages to estimate the age-range for each spectroscopically known star. This result can be compared with the age distribution of this sample that was observed in the UV by Crowther et al. (2016, their Fig. 11).

Fig. 9 Generalized histogram of the extinction of 54 spectroscopically known stars from Crowther et al. (2016, shown in Fig. 7). We used PARSEC models to estimate the extinction for each of these stars according to its age range.

Fig. 10 Generalized histogram of the stellar masses (MF) at 0.6, 1 and 1.5 Myr. PARSEC models used to estimate the stellar-mass range for each source using extinction-range.

5. Visual companions

For each star detected in both J and Ks, we determined the distance between the star and its closest neighbor. Figure 11 shows the number of closeby detected stars in Ks (red) and J (blue) versus their separation in arcseconds. More than 250 star pairs have a closest neighbor at a separation smaller than 0.2′′. Over 90% of massive objects (brighter than 16.7 mag in Ks and 15 mag in J) have a closest neighbor with a separation smaller than 0.2′′. Figure 11 shows the separation between the visual close stars versus their distance from R 136a1 in the core. This figure indicates that even the sources at larger radii have close visual companions, so that the large number of close visual companions in not just an effect of 2D projection on the sky across the FoV. For the sake of simplicity, regardless of whether the stars are physically bound, we call these closeby stars “visual companions” hereafter.

The most massive stars R 136a1, R 136a3, and R 136c have visual companions. R 136a3 is resolved as two stars with the PSF fitting. Both stars have a high correlation coefficient (above 70%) with the input PSF. The separation between the R 136a3 primary and secondary is about 58.9 ± 2.14 mas, which is larger than the FWHM of the PSF. We note that even the closest visual companions (like R 136a3) are physically far from each other (0.059′′ is 2890 AU). This visual separation produces a period over P = 10⁴ yr, so that these sources are probably not gravitationally bound to each other. Table 2 shows the list of the 20 stars detected in both J and Ks data that have visual companions closer than 0.08′′. The flux ratio between two companions in Ks and J band are given in the third and fourth column, respectively. Their separation [in mas] is also given in the last column.

Fig. 11 Top: histogram of the separation of the close detected sources. For each star that is detected in both J and K data, we determined a distance between the star and its closest neighbor. Bottom: separation of the visual close detected sources versus their distance from the core of R 136.

6. Density and surface brightness profile

The unexpected number of detected sources in a small FoV and of resolved companions in R 136 indicates that this compact cluster is more crowded than thought before. The error bars on the stellar masses and also on the age and extinction of the cluster itself are large enough to make it difficult to study the density profile of R 136. Instead, we can scrutinize the surface brightness profile (SBP) of this cluster, which is less affected by the confusion and crowding. We obtained the SBP by measuring radial profiles for J and Ks images centered on R 136a1. The SBP informs us on the average magnitude per pixel at different radii. Figure 12 depicts the SBP of the core of R 136 in J and K. On average, the SBP in Ks is brighter than J, which can be caused by the extinction or brighter stars in Ks. Several bumps on the SBP are visible at approximately 0.08, 0.15, 0.46, 1.17 pc radii. These radial distances are the locations of known WR stars. The position of five WRs in the FoV is shown in the SBP plot. These stars have extensive emissions in the Ks band because of their wind and mass loss.

Fig. 12 SBP (mag/pixel) of R 136 in IRDIS FoV in IRDIS FoV centered on R 136a1. The radial positions of the five WR stars are shown with the solid vertical purple lines.

Fig. 13 Projected mass density [M⊙/pc2] profile of R 136 in IRDIS FoV centered on R 136a1. The stellar masses are estimated at the age of 1 Myr with extinction values of AJ = (0.45 ± 0.5) mag and AK = (0.2 ± 0.5) mag in J and Ks band. Equation (1)is used to fit the blue solid line to the data.

Using the stellar masses estimated at the age of 1 Myr and extinction values in J and K, A_J = (0.45 ± 0.5) mag and A_K = (0.2 ± 0.5) mag, we plot the 2D (projected) density profile (Fig. 13). We used an Elson-Fall-Freeman (EFF) profile (Elson et al. 1987) to fit the projected mass density in the core of R 136 (Eq. (1)), $\mathit{\rho }\mathrm{[}{\mathit{M}}_{\mathrm{\odot }}\mathit{/}{\mathrm{pc}}^{\mathrm{2}}\mathrm{]}\mathrm{=}\frac{{\mathit{\rho }}_{\mathrm{0}}}{{\mathrm{(}\mathrm{1}\mathrm{+}{\frac{{\mathit{r}}^{\mathrm{2}}}{{\mathit{a}}^{\mathrm{2}}}}^{\mathrm{)}}}^{\frac{\mathit{\gamma }\mathrm{+}\mathrm{1}}{\mathrm{2}}}}\mathrm{·}$(1)We estimated the central mass density of ${\mathit{\rho }}_{\mathrm{0}}\mathrm{=}\mathrm{\left(}\mathrm{3.8}{\mathrm{9}}_{-1.14}^{\mathrm{+}\mathrm{1.60}}\mathrm{\right)}\mathrm{\times }{\mathrm{10}}^{\mathrm{4}}\mathrm{[}{\mathit{M}}_{\mathrm{\odot }}\mathit{/}{\mathrm{pc}}^{\mathrm{2}}\mathrm{]}$ and the parameters γ = 1.23 ± 0.24 and a = 0.17 ± 0.06. The total observed mass of the clusters for r< 1.4 pc is ${\mathit{M}}_{\mathrm{obs}}\mathrm{=}\mathrm{\left(}\mathrm{1.3}{\mathrm{3}}_{-0.22}^{\mathrm{+}\mathrm{0.30}}\mathrm{\right)}\mathrm{\times }{\mathrm{10}}^{\mathrm{4}}{\mathit{M}}_{\mathrm{\odot }}$, down to 2 M_⊙. The total mass of the cluster depends on the shape of the MF and the lowest mass limit. The real stellar masses remain open because we are limited to the angular resolution of SPHERE/IRDIS, so that the estimated central projected density is a lower limit to the real central density.

We note that the estimated density is projected in 2D. In order to estimate the 3D density approximately, we considered R 136 to be spherically symmetric and have a radius of R_cluster. The density profiles were estimated for different R_cluster values from 2 pc to 6 pc. In this way, the 3D central densities, γ, and a were computed by fitting the EFF profile (Eq. (1)). Table 3 shows the fitting (3D) parameters for R 136 considering different values of R_cluster. The total mass of the cluster can be estimated by extrapolating to the considered R_cluster. The ratio of the observed total mass within r< 1.4 pc to the total mass estimated of the cluster within a given radius (R_cluster) also is given in the last column of Table 3.

The estimated value of γ and a in 2D and 3D are consistent and the shape of the densities is flatter than the Plummer model (close to the King model). All the 3D central densities γ and a are lower than the previous values given by Mackey & Gilmore (2003) and Selman & Melnick (2013).

7. Discussion and conclusion

We presented a photometric analysis of the core of R 136 using the VLT/SPHERE instrument in the near-IR. The high quality and resolution of these data open a new perspective on our understanding of R 136. For the first time, more than one thousand sources have been detected in Ks- and J-band data in the small FoV of IRDIS (10.9″ × 12.3″) covering almost 2.7 × 3.1 pc of R 136 core. For the ground-based telescopes, the best data come from VLT/MAD, whose the AO quality (Strehl ratios of 15–30% in Ks) is poorer than that of SPHERE (Strehl ratios of 75% in Ks). The confusion, especially in the core, therefore remains large enough for the sources that are undetectable in its core.

In SPHERE/IRDIS data, more than 60% of the stars have companions closer than 0.2′′(0.05 pc). Ninety percent of the very massive bright stars that have been studied spectroscopically by Crowther et al. (2016) have visual companions. The large error bars on the spectroscopic parameters (T_eff and log L) prevent us from estimating the age, extinction, and stellar masses accurately. Using the stellar parameters by Crowther et al. (2016) from the UV spectroscopically analysis, the most probable age of the core is about $\mathrm{1.}{\mathrm{8}}_{-0.8}^{\mathrm{+}\mathrm{1.2}}$ Myr (Fig. 8), which is consistent with the Crowther et al. (2016) values. We find the best age estimate to be about 1 and 1.5 Myr for the very core (r< 3″) and outside (r> 3″) of R 136, however, from fitting isochrones with CMD. We clearly see the conjunction of pre-main sequence and main sequence, which is very sensitive to the age. The extinction in J and Ks are A_J = (0.45 ± 0.5) mag and A_K = (0.2 ± 0.5) mag, respectively. The extinction was estimated from stellar parameters of 54 spectroscopically studied stars by Crowther et al. (2016). The large errors on the extinction values therefore come from the large errors on the spectroscopical analysis on log L and T_eff.

Considering the photometric errors, the stellar masses are estimated at different ages with a broad extinction range. The MF slope for the 1 and 1.5 Myr isochrone is Γ_1Myr = −0.90 ± 0.13 and Γ_{1.5 Myr} = −0.98 ± 0.18, respectively, for the mass range of (3–300) M_⊙.

As the core becomes better resolved, more stars are detected. The MF slope is flatter than the Kroupa (Γ = −1.3) and Salpeter (Γ = −1.35) values. Thanks to the SPHERE high-contrast data, we estimated the MF from 3.0 M_⊙ up to 300 M_⊙ at different ages. The estimated MF slopes are consistent with the values estimated by Malumuth & Heap (1994), Hunter et al. (1996), and Selman et al. (1999) for the core of R 136. These values are flatter than the values estimated by Massey & Hunter (1998) and Brandl et al. (1996). Table 4 shows the MF slopes estimated for R 136 in previous works for different mass range and regions. The direct comparison is not logical since 1) none of the previous observations, within the IRDIS FoV, covers the mass range of this study and 2) our MF slopes are estimated considering the errors on each stellar mass (see Sect. 4), while the previous studies did not apply the stellar mass errors in the MF.

The derived MF is limited to the resolution of the instrument and also to the detection limit of the observation. Higher angular resolution data may resolve binaries and low-mass stars that affect the shape of the MF. Figure 13 shows the density profile of the R 136 core at 1 Myr (0.1 pc <r< 1.4 pc). The lower limit of the central density of R 136 is ${\mathit{\rho }}_{\mathrm{0}}\mathrm{=}\mathrm{\left(}\mathrm{3.8}{\mathrm{9}}_{-1.14}^{\mathrm{+}\mathrm{1.60}}\mathrm{\right)}\mathrm{\times }{\mathrm{10}}^{\mathrm{4}}\mathrm{[}{\mathit{M}}_{\mathrm{\odot }}\mathit{/}{\mathrm{pc}}^{\mathrm{2}}\mathrm{]}$ at 1 Myr, which is about ρ₀ = 288 [stars/pc²]. As the mass estimate in the very core of R 136 is biased by crowding and also for massive WR stars, we excluded the first point from fitting the density profile. The observed total mass of R 136 for r< 1.4 pc is ${\mathit{M}}_{\mathrm{obs}}\mathrm{=}\mathrm{\left(}\mathrm{1.3}{\mathrm{3}}_{-0.22}^{\mathrm{+}\mathrm{0.30}}\mathrm{\right)}\mathrm{\times }{\mathrm{10}}^{\mathrm{4}}{\mathit{M}}_{\mathrm{\odot }}$ which covers the mass range of (2–300) M_⊙. This mass is about 51% of the whole mass range of (0.01–300) M_⊙ using the Chabrier MF (Chabrier 2003, 2005).

Considering that R 136 is a spherically symmetric cluster with radius R_cluster (Table 3), we estimated the 3D density profile. The 3D central densities are lower than the values estimated in previous studies given by Mackey & Gilmore (2003) and Selman & Melnick (2013). The computed values of γ and a (Eq. (1)) are consistent in 2D and 3D considering different R_cluster. All density profiles are flatter than the Plummer model (γ = 4.0).

Very massive stars in R 136 have similar characteristics as the galactic WR stars in the core of NGC 3603. NGC 3603 is almost eight times closer than R 136. The effect of confusion is visible in Fig. A.3, in which we visualize NGC 3603 at the distance of R 136 (Fig. A.3 middle). Using Starfinder we detected 408 and 288 stars in J and Ks images, respectively (Khorrami et al. 2016). Using the same criteria for Starfinder, we only detected 109 and 52 sources in J and Ks images of NGC 3603 at the distance of R 136, which means that more than 70% of the stars cannot be detected. This implies that about 1000 detected stars in the R 136 core (r< 6″) are possibly only 30% of the real number. The average density in this region (r< 6″) would increase from 71 [star/pc²] to 230 [star/pc²]. The lack of resolution prevents us from accurately estimating stellar masses, core density, and density profile, while SBP is less affected. Figure A.4 shows the SBP of NGC3603 from IRDIS data in J and Ks band, directly and simulated as it would be located at a distance of R 136. The general trend is not affected.

Using SPHERE data, we have gone one step further and partially resolved and understood the core of R 136, but this is certainly not the final step. R 136 needs to be observed in the future with higher resolution (E-ELT) and/or a more stable PSF (JWST), therefore deeper field imaging. The cluster would then be better characterized for its age, individual and multiple stars, and ultimately, its kinematics on a long enough temporal baseline of observation.

Acknowledgments

Z.K. is supported by the Erasmus Mundus Joint Doctorate Program by Grant Number 2012-1710 from the EACEA of the European Commission and the University of Nice foundation. The authors are indebted to Dr. A. Stolte for her comments and suggestions for improving the original manuscript. SPHERE is an instrument designed and built by a consortium consisting of IPAG (Grenoble, France), MPIA (Heidelberg, Germany), LAM (Marseille, France), LESIA (Paris, France), Laboratoire Lagrange (Nice, France), INAF–Osservatorio di Padova (Italy), Observatoire de Genève (Switzerland), ETH Zurich (Switzerland), NOVA (The Netherlands), ONERA (France) and ASTRON (The Netherlands), in collaboration with ESO. SPHERE was funded by ESO, with additional contributions from CNRS (France), MPIA (Germany), INAF (Italy), FINES (Switzerland) and NOVA (The Netherlands). SPHERE also received funding from the European Commission Sixth and Seventh Framework Programmes as part of the Optical Infrared Coordination Network for Astronomy (OPTICON) under grant number RII3-Ct-2004-001566 for FP6 (2004-2008), grant number 226604 for FP7 (2009-2012) and grant number 312430 for FP7 (2013-2016). We acknowledge financial support from the Programme National de Planétologie (PNP) and the Programme National de Physique Stellaire (PNPS) of CNRS-INSU. This work has also been supported by a grant from the French Labex OSUG@2020 (Investissements d’avenir – ANR10 LABX56). We thank the anonymous referee for constructive comments, which improved the paper.

References

Appendix A: Additional table and figures

Fig. A.1 Map of the photometric, residual, and total relative error bars. Top: map of the PSF fitting errors (outcome pf the Starfinder photometry) along the IRDIS FoV. Middle: map of the residual errors, outcome of the background analysis after removing the stellar source signals from the image. Bottom: map of the total relative error, combination of the PSF fitting errors and the residual background errors. Left images are in the J band. Right: images are in the Ks band.

Fig. A.2 Surface incompleteness map for magnitude 18 (left) and 19 (right) in J band. The incompleteness increases toward the central crowded core of R 136 and also close to the bright stars. The black pluses indicate the position of the stars that are almost as bright as the artificial star test.

Fig. A.3 Comparison of NGC 3603 and R 136 images from VLT/SPHERE. Left: core of R 136 (1.56″ × 1.56″) at its real distance. Right: core of NGC3603 (12.5″ × 12.5″) at its real distance. Middle: NGC 3603 as it would appear at the same distance as R 136. Upper and bottom panels are IRDIS Ks- and J-band images, respectively.

Fig. A.4 SBP of NGC 3603 in SPHERE/IRDIS J- and Ks-band data. Left: NGC 3603 real images. Right: simulated NGC 3603 at the distance of R 136.

All Tables

All Figures

	Fig. 1 Comparison of R 136 core images at different wavelengths with the highest available angular resolution telescopes. The FoV of all images is the same as IRDIS data (10.9″ × 12.3″). Top left: HST/WFPC2 in V band (λcen 526 nm), top right: VLT/MAD in Ks band (λcen 2200 nm), bottom right: SPHERE/IRDIS/Ks (λcen 2182 nm); bottom left: SPHERE/IRDIS/J (λcen 1245 nm).
In the text

	Fig. 2 Map of the correlation coefficients between the input PSF and each detected stellar source, calculated by Starfinder across the FoV in J (top) and in K (bottom) images of R 136 by SPHERE.
In the text

	Fig. 3 PSF fitting errors (red pluses) and residual errors (blue crosses) in Ks (top) and J (bottom). The PSF fitting error is the outcome of the Starfinder. The residual errors are the outcome of the background analysis after removing the stellar sources from the image.
In the text

	Fig. 4 Error in the position of the stars in Ks (top) and J (bottom) data. This error is the outcome of the Starfinder and decreases for the bright stars.
In the text

Fig. 5 Incompleteness test in J (red) and Ks (blue) for two regions: the very core of R 136 (r < 3″), shown with pluses and crosses, and outside of the core (r > 3″), shown with circles and triangles. For each flux value (magnitude) we used 500 artificial stars in order to find the completeness value. Note that artificial stars are added one by one to the original image, and there is always one artificial star in each experiment.

In the text

Fig. 6 Top: CMD of 818 detected sources in J- and Ks-band images of SPHERE/IRDIS from the core of R 136. Solid black, green and blue lines show the PARSEC isochrones at the ages of 0.6, 1 and 1.5 Myr (corrected for distance modulus of 18.49 and central values of extinctions, AJ = 0.45 mag and AK = 0.2 mag). The horizontal dashed gray line depicts the incompleteness of 70%. The CMD is plotted for the whole FoV (818 sources), in the very core of the cluster (r < 3″) and outside (r > 3″), from left to right respectively. The error bar on each point is the combination of the PSF fitting errors and the residual errors from the background image after removing the stellar source signals from the images. The 54 spectroscopically studied stars by Crowther et al. (2016) are shown as black dots in the middle plot. Bottom: same as top CMD, but in 3D to show the number density of stars in the CMD. Solid red, black, pink, and yellow lines depict PARSEC isochrones at the age of 0.5, 1, 1.5, and 2 Myr, respectively.

In the text

	Fig. 7 Left panel: the IRDIS/Ks (FoV 6″ × 6″) on which the 54 spectroscopically known stars from Crowther et al. (2016) have been added as red circles. The right plot depicts the Teff, log L/L⊙, and corresponding error bars on these 54 sources taken from Crowther et al. (2016) in blue. The solid red lines indicate the PARSEC isochrones covering ages from 0.1 to 8 Myr.
In the text

	Fig. 8 Generalized histogram of the age of 54 known stars from Crowther et al. (2016, shown in Fig. 7). We used PARSEC isochrones at different ages to estimate the age-range for each spectroscopically known star. This result can be compared with the age distribution of this sample that was observed in the UV by Crowther et al. (2016, their Fig. 11).
In the text

	Fig. 9 Generalized histogram of the extinction of 54 spectroscopically known stars from Crowther et al. (2016, shown in Fig. 7). We used PARSEC models to estimate the extinction for each of these stars according to its age range.
In the text

	Fig. 10 Generalized histogram of the stellar masses (MF) at 0.6, 1 and 1.5 Myr. PARSEC models used to estimate the stellar-mass range for each source using extinction-range.
In the text

	Fig. 11 Top: histogram of the separation of the close detected sources. For each star that is detected in both J and K data, we determined a distance between the star and its closest neighbor. Bottom: separation of the visual close detected sources versus their distance from the core of R 136.
In the text

	Fig. 12 SBP (mag/pixel) of R 136 in IRDIS FoV in IRDIS FoV centered on R 136a1. The radial positions of the five WR stars are shown with the solid vertical purple lines.
In the text

	Fig. 13 Projected mass density [M⊙/pc2] profile of R 136 in IRDIS FoV centered on R 136a1. The stellar masses are estimated at the age of 1 Myr with extinction values of AJ = (0.45 ± 0.5) mag and AK = (0.2 ± 0.5) mag in J and Ks band. Equation (1)is used to fit the blue solid line to the data.
In the text

Fig. A.1 Map of the photometric, residual, and total relative error bars. Top: map of the PSF fitting errors (outcome pf the Starfinder photometry) along the IRDIS FoV. Middle: map of the residual errors, outcome of the background analysis after removing the stellar source signals from the image. Bottom: map of the total relative error, combination of the PSF fitting errors and the residual background errors. Left images are in the J band. Right: images are in the Ks band.

In the text

	Fig. A.2 Surface incompleteness map for magnitude 18 (left) and 19 (right) in J band. The incompleteness increases toward the central crowded core of R 136 and also close to the bright stars. The black pluses indicate the position of the stars that are almost as bright as the artificial star test.
In the text

	Fig. A.3 Comparison of NGC 3603 and R 136 images from VLT/SPHERE. Left: core of R 136 (1.56″ × 1.56″) at its real distance. Right: core of NGC3603 (12.5″ × 12.5″) at its real distance. Middle: NGC 3603 as it would appear at the same distance as R 136. Upper and bottom panels are IRDIS Ks- and J-band images, respectively.
In the text

	Fig. A.4 SBP of NGC 3603 in SPHERE/IRDIS J- and Ks-band data. Left: NGC 3603 real images. Right: simulated NGC 3603 at the distance of R 136.
In the text