High-contrast and resolution near-infrared photometry of the core of R136

We present the sharpest and deepest near-infrared photometric analysis of the core of R136, a newly formed massive star cluster at the centre of the 30 Doradus star-forming region in the Large Magellanic Cloud. We used the extreme adaptive optics of the SPHERE focal instrument implemented on the ESO Very Large Telescope and operated in its IRDIS imaging mode for the second time with longer exposure time in the H and K ﬁlters. Our aim was to (i) increase the number of resolved sources in the core of R136, and (ii) to compare with the ﬁrst epoch to classify the properties of the detected common sources between the two epochs. Within the ﬁeld of view (FOV) of 10.8 (cid:2)(cid:2) × 12.1 (cid:2)(cid:2) (2 . 7 pc × 3 . 0 pc), we detected 1499 sources in both H and K ﬁlters, for which 76 per cent of these sources have visual companions closer than 0.2 (cid:2)(cid:2) . The larger number of detected sources enabled us to better sample the mass function (MF). The MF slopes are estimated at ages of 1, 1.5, and 2 Myr, at different radii, and for different mass ranges. The MF slopes for the mass range of 10–300 M (cid:3) are about 0.3 dex steeper than the mass range of 3–300 M (cid:3) , for the whole FOV and different radii. Comparing the JHK colours of 790 sources common in between the two epochs, 67 per cent of detected sources in the outer region ( r > 3 (cid:2)(cid:2) ) are not consistent with evolutionary models at 1–2 Myr and with extinctions similar to the average cluster value, suggesting an origin from ongoing star formation within 30 Doradus, unrelated to R136.


INTRODUCTION
R136 is a very massive young star cluster that lies at the centre of the Tarantula nebula in the Large Magellanic Cloud (LMC).Hosting the most massive stars known in the Local Universe (Crowther et al. 2010;Crowther et al. 2016), R136 provides a unique opportunity to observationally study the formation of massive stars and clusters in the earliest stages of their evolution.Our under-2 Z. Khorrami et al. planet Research 1 (SPHERE, Beuzit et al. 2019) instrument of the Very Large Telescope (VLT).Thanks to SPHERE's high contrast and extreme adaptive optics (XAO), the sharpest images from the central region of R136 were recorded in J and K.For the first time, more than one thousand sources were detected in these bandpasses within the small field of view (FOV) of IRDIS (10.9" × 12.3") covering almost 2.7 × 3.1 pc of R136 core (Khorrami et al. 2017).The SPHERE data on R136 core were used to partially resolve and study the initial mass function (IMF) covering a mass range of 3 to 300 M , at ages of 1 and 1.5 Myr.The density in the core of R136 ( < 1.4 pc) was estimated and extrapolated in 3D for larger radii (up to 6pc).Even at high angular resolution in the NIR the stars in the core are still unresolved due to crowding and central concentration of bright sources.Using evolutionary models of stars more massive than 40 solar masses, Bestenlehner et al. (2020) finds that the initial mass function (IMF) of massive stars in R136 is topheavy with a power law exponent  = 2.0 ± 0.3.They also estimate the age of R136 between 1 and 2 Myr with a median age of around 1.6 Myr.
This paper presents the second epoch observation of R136 in the NIR by VLT/SPHERE.These observations, obtained with longer exposures in different NIR filters, aimed first to increase the number of resolved sources in the core of R136, and second to enable comparison with the first epoch to classify the property of the detected common sources between these two epochs.We divided this paper as follows.In Section 2 we explain in detail the data and the observing conditions.We then present the photometric analysis in Section 3 and the method we use to correct for completeness in Section 4. In Section 5 we describe the details of extinction estimation.On the basis of the results we investigate the stellar mass functions and the density in Section 6.In Section 7, we compare the photometric analysis between the two epochs and between the two sets of filters.We conclude by summarising our results in Sec. 8.

OBSERVATIONS
We obtained the data via ESO Time Observation run (0102.D-0271) to image R136 using the "classical imaging" mode of IRDIS Langlois et al. (2014) in H-and K-bandpass filters.For our purpose we used the same spectral band split into the two IRDIS channels to correct for residual detector artefacts such as hot pixels and uncorrelated detector noise among other instrumental effects.We maximized the image sharpness of the total exposure by discarding the single frames with poorer Strehl ratio and by correcting "a posteriori" the residual tip-tilt image motion on each short exposure before combining them.Observations were performed in October/November (K/H) 2018, achieving high dynamical range and high angular resolution imaging in both the H and K bandpass, over a FOV of 10.8" × 12.1" (2.7pc × 3.0pc), centered on the core of the cluster.The natural seeing was 0.69 ± 0.10 arcseconds in K and 0.58 ± 0.05 in H during the observations.The night was rated as "Clear" meaning that less than 10% of the sky (above 30 degrees elevation) was covered by clouds, and the transparency variations were less than 10% during the exposures.
Our data consists of 1088/544 frames of 2.0/4.0sexposures taken with the IRDIS broad-band H and K filters (BB-H, BB-K).and the filters information.Figure 1 shows the final reduced images in the H (left) and K (right).These images are the deepest and sharpest images taken from the core of R136 so far.The Wolf-Rayet star R136a1 was used for guiding the AO loop of SPHERE confirming the high level of performances even for faint guide stars which surpass both NACO and MAD performances.Figure 2 compares the highest resolution available images from the very central part of R136 taken by HST (WFPC2 by Hunter et al. 1995) in left (top in V, bottom in H), VLT/MAD (Campbell et al. 2010) in the middle (top in K, bottom in H), and IRDIS in right (top in K, bottom in H).
These images shows the effect of better contrast, pixel-sampling, sensitivity, and AO correction on resolving the stellar sources in the most crowded part of R136.
The range of airmass during these observations was 1.52 to 1.45 in K, and 1.70 to 1.55 in the H.A log of observations is provided in Table 2.For comparison we include the observing log for our first epoch data-sets taken in 2015.
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 K).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 mas.The final step is performed to combine the parallel images on the left and right side of the detector, after the anamorphism correction.Some uncorrected background leaks persist in our final K images due to thermal background fluctuations, which are stronger in K than in H.

PHOTOMETRY
To analyse the final images we applied the same method/tools as on first epoch data (Khorrami et al. 2017).For the photometry we used the Starfinder package implemented in IDL (Diolaiti et al. 2000).Starfinder is designed for the analysis of AO images of crowded fields, like the Galactic Center, for instance, as in Pugliese et al. (2002).This method determines the empirical local Point Spread Function (PSF) from several isolated sources in the image and uses this PSF to extract other stellar sources across the FOV.
For the present study, four and seven well isolated stars (within 0.47" × 0.47") were used to extract the PSF in the K and H, respec-   tively.The FOV in K is more crowded than in H (see Figure 1), so finding the bright isolated stars (within 0.47" × 0.47") in K is harder than in H.That is why we limit our selection of stars to be used for PSF extraction, to four in K. Figure A3 shows the separation with the closest neighbour for each star in H and K.The extracted PSFs, were used as an input for stellar sources detection by Starfinder.The FWHM of the extracted PSFs are 58.8 and 63.7 mas in H and K, with 1658 and 2528 sources detected, respectively.We stopped the source extraction after obtaining a minimum correlation coefficient of 65% and 80%, in K and H, between the extracted star with the locally determined PSF according to Starfinder procedures.In H, we put a higher threshold on minimum correlation coefficient than in K, because the AO correcting efficiency in H degrades faster as a function of distance from R136a1 -which is the reference star for the AO loop-than in K.The isoplanatic angle in H is smaller than in K, so at larger radii from R136a1, the PSF is not centrosymmetric.Indeed, stars with higher correlation coefficients, i.e. more similarity to the PSF, represent higher reliability on their photometry estimation.
In addition to the correlation coefficient criterion, we applied the limit of standard deviation from the sky brightness (   ) for stopping the extraction of sources by Starfinder, i.e. the local PSF maximum value must exceed 2   over the background.Figure 3 shows the signal to noise ratio (SNR) of 1658 and 2528 detected sources in H (purple circles) and K (green squares), respectively.Common sources (1499) between H and K are shown as filled circles/squares.The solid horizontal line shows the SNR=2.0where some of the detected sources have SNR less than this value.These stars are located further from the central region where the local    is smaller than the    of the whole FOV.76% of these sources have visual companions closer than 0.2" which is much higher than the value found in the first epoch in 2015 data.Figure 4 shows the histogram of closest neighbour (visual companion) separation from each source versus its distance from the centre of R136, in 2015 (top) and 2018 (bottom) data.Comparing the histogram of closest companions between two epochs, 85% of the new detected sources in 2018 have visual companions with separation less than 0.2".
The Strehl Ratio (SR) in H and K is determined as 0.71 ± 0.05 and 0.83 ± 0.03, respectively, on average within 5 arcseconds from the core.These estimates are assessed from the SPARTA files recorded during SPHERE runs simultaneously with the AO corrected images of R136.They are based on the slope measurements of the Shack-Hartmann wavefront sensor in Sphere AO for eXoplanet Observation (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 quite robust for FOV smaller than the anisoplanetism in the past for differential imaging by NACO (Maire et al. 2014).We are therefore quite confident on our photometry corrected for the SR effect in H and K. To convert stellar fluxes to magnitudes, we set the zeropoint magnitude to the instrumental zeropoint (one ADU/s in H and K, are 25.2 and 24.8 magnitudes, respectively).To double check our calibration, we compared H and K magnitudes of 26 sources (located  ≥ 3" from R136a1) common between our catalogue and the VLT/MAD (Campbell et al. 2010) in H and K2 .We cannot compare our zeropoints with the MAD/H and K magnitudes in the central 2.8" region because of the completeness issue and the very low AO correction (SR ∼ 15-30%) of their data.Figure 2 shows the completeness and contrast problem of the previous data, leading to an inability to distinguish certain sources.By way of example R136a6 (H19) and H263 -are very close in terms of brightness (flux ratio of about 90%) and location (Separation of about 70mas).These bright O-type stars are most easily detectable in our data compared with others (Figure 2 and Figure A2).Cases like these close bright stars, bring confusion even in the HST spectroscopic analysis (Bestenlehner et al. 2020).
The K magnitude of R136a1 and R136a2 (11.15 mag and 11.43 mag) increased about 0.08 and 0.11 mag, compare to their Kmag in the first epoch (11.07 mag and 11.32 mag).The brightness of R136a3 and R136b in K (11.45 mag and 11.67 mag) are consistent with the first epoch K data (11.44 mag and 11.66 mag).R136c in the K (11.31 mag) is 0.17 mag brighter than in 2015 (11.48 mag), making this source the second brightest in our K catalogue.In previous studies R136c shows indications of binarity based on its strong X-ray emission (Portegies Zwart, Pooley, & Lewin 2002;Townsley et al. 2006;Guerrero & Chu 2008) and possible lowamplitude RV variations (Schnurr et al. 2009).
Comparing the brightness of common detected sources within two epochs, the K magnitude of 56% of the sources have changed more than twice their photometric errors 4 .To make sure that this is not due to Starfinder photometry analysis (which is based on the reference input PSF), we have used DAOPHOT algorithm for this analysis as well.The variation of the K magnitude of these sources still remains 5 .The K mag of about 67% of the sources located outside of the core ( > 3"), has changed more than twice their photometric errors.Since the observations were made in different nights and the AO correction in the outer part of the FOV is not as good as its center, this inconsistency is probably an observational effect, otherwise it is originated from physical reasons like variables or multiple systems with periods less than 3 years.

COMPLETENESS
We produced 6.2 × 10 4 artificial stars (500 per magnitude) in H and K, in order to determine the completeness as a function of magnitude.These artificial stars are created using the same PSF as we used for source extraction process.They are added to the original image one by one (500 times per magnitude) and the same photometric tools and criteria to determine how often they can be recovered from the source extraction process are used.These tests were performed one star at a time to avoid the artificial stars from affecting one another.Figure 5 shows the completeness values as a function of magnitude in H (purple dots) and K (green squares).
The completeness is 80% at  = 20 mag and  = 20.5 mag.The completeness value in each magnitude is the average value for detection of 500 artificial stars which are distributed randomly in the FOV.But the completeness varies across the FOV depending on how crowded and bright that region is.Using these 6.2 × 10 4 artificial stars, we created completeness maps covering magnitude range of 17 to 22 (5 maps), both in H and in K. Figure 6 shows two examples of completeness map in the magnitude of 19 to 20 both in H (top) and K (bottom).The average value of completeness is about 90% both in H and K (Figure 5), but one can notice that the completeness in the very central part of the cluster is much lower than the outskirts.In the K we can also see the effect on the completeness of the thermal background noise from the instrument lying in the arcs located in the far east and west parts of the FOV as seen on Figure 6.The black pluses in Figure 6 shows the position of the observed sources (398 in H and 410 in K) within the magnitude range of 19 to 20.This plot clearly shows the effect of completeness on the detected sources which are used to create MF.To overcome this limitation we corrected the MF for all stars fainter than 17 magnitude.Each star contributes a mass distribution to the MF (see Section 6).When we added stars with a given mass to the histogram of mass, we corrected the contribution of each star to the mass function (fainter than 17 mag) for completeness using the completeness maps in every H and K magnitudes.We know that the completeness in a given magnitudes varies across the FOV in both filters, so depending on the position and magnitude of an observed star, we estimated its completeness using these maps (see Figure 6).Since the average value of completeness at 17 mag is above 95% (above 80% at the core), we stopped the completeness correction for sources brighter than K=17 mag. 4 see Figure A4 to compare the K mag of these sources within two epochs.

EXTINCTION
In order to fit the evolutionary models to our data and estimate the stellar masses, we need to measure the extinction first.We used two methods to estimate the extinction for the spectroscopicallyobserved massive stars in the core of R1366 .In both methods, we adopt the LMC distance modulus (DM) of 18.49 magnitude estimated by Pietrzyński et al. (2013) which is consistent with the value suggested by Gibson (2000) for LMC.
Method I: We used the effective temperature (T eff ) and luminosity (logL) of 49 stars by Bestenlehner et al. (2020) in the optical.The values of these temperatures and luminosities are provided in Table A1.We also chose a grid of isochrones at different ages (from 0.1 up to 10 Myr) with the LMC metallicity (Z=0.006),from the latest sets of PARSEC evolutionary model 7 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 .By fitting the PAR-SEC isochrones to the observed stellar parameters of each star (T eff , logL), we estimated the age and intrinsic colour for each star with an error.The extinction in H and K (A  and A  ) are estimated by comparing these intrinsic H and K magnitudes with the observed values from our catalogue.Figure 7 shows the histogram of the extinction values of those 49 stars.The mean extinction in H and K, and the mean colour excess are A  = 0.38 ± 0.55, A  = 0.22 ± 0.57, and E(H-K)= 0.11 ± 0.21.These values derived by fitting a Gaussian function on the extinction distributions, while these distributions are not exactly Gaussian and symmetric.That is why the mean colour excess is slightly (0.05 mag) higher than A  − A K .The large errors on these values results from the errors in the reported stellar parameters (T eff and log L) from the spectroscopic analysis in the optical (see Table A1).
Method II: Using the stellar parameters from Bestenlehner et al. (2020), we adopt the closest stellar atmosphere model of that star.2013) and with the one used in previous studies by Campbell et al. (2010).

MASS FUNCTION AND CORE DENSITY
Figure 9 shows the colour magnitude diagram (CMD) of 1499 detected sources common between the H and K data (top) for the full FOV in left (N=1499), for inner  < 3" region (N=627) in middle, and for outer  > 3" region (N=872) in right11 .Red, yellow and blue solid lines are PARSEC isochrones at DM=18.49 with   = 0.35 and  ( − ) = 0.1, at the ages of 1, 1.5, and 2 Myr, respectively.We note that there is a large scatter in the CMD, especially for the lower part of the main sequence and pre-main sequence sources which are located (detected) in the outer region.A significant scatter in CMD and colour-colour diagrams of 30 Doradus has previously been reported in the visible and near infrared (Hunter et al. 1995;Andersen et al. 2009;De Marchi et al. 2011;Cignoni et al. 2015).The scatter is likely due to a combination of observational confusion (affected mainly by the visual multiple systems or variables), photometric errors, differential extinction, a possible age spread, colour excess due to warm circumstellar matter and ongoing star formation within 30 Doradus.Pre-main sequence stars are often associated with circumstellar disks and outflows which will introduce additional extinction for the clusters low-mass content.Brandl et al. (1996) found that the extinction varies significantly from star to star within the cluster, with the range of 1-2 mag.The HST observations also reveal the presence of considerable differential extinction across the 30 Doradus region.De Marchi et al. (2011) quantified the total extinction toward massive main sequence stars younger than 3 Myr to be in the range 1.3 <   < 2.1.Figure 9 shows the evolutionary models (isochrones) can not be fitted for these scattered sources.We prefer to consider a (larger) error of 0.5 mag on the extinction, to estimate the stellar mass range for each star at a given age.We estimated the stellar masses just for common sources between H and K data (1499 sources) using both their H and K magnitudes fitted to PARSEC isochrones at three different ages: 1, 1.5, and 2 Myr.
The mass function (MF), is plotted in Figure 10 considering a Gaussian distribution for each stellar mass.We convert the photometric/extinction uncertainty into a mass (gaussian) distribution.This Gaussian uncertainty in the mass of each star is accounted for, when constructing the MF.Each star fainter than K=17 (about 13.4 M at 2 Myr) is corrected for completeness, according to its brightness and its location in the FOV (see Section 4 for more information).Figure 10 shows the MF at three different ages (1, 1.5, and 2 Myr), for the whole FOV (top), central  ≤ 3" region (middle), and outer r > 3" region (bottom).The MF slope values in these regions, at different ages and mass ranges are provided in Table 3, both for the completeness corrected (Γ  ) and not corrected MF (Γ   ).The MF slopes are calculated for the minimum mass of 3 M (about K=19.8 at 2 Myr) where the average value of completeness is above 85% (Figure 5) but the difference between Γ  and Γ   is not negligible for the inner-region of the cluster where the completeness reaches to 30% (see Figure 6) and this affects the MF slope for the whole FOV (top values in Table 3).The slopes for the mass range of 3−300 M for the inner region is flatter than the outer region and consequently for the whole FOV.This can be explained by completeness since the Γ  and Γ   are not compatible.But for the mass range of 10 − 300 M where the data are complete (see Figure 5 and Figure A1) and the MFs are not affected by completeness (compare Γ  and Γ   in table 3), the slope values in the inner region is slightly flatter than the outer region, and consequently than the whole FOV.Comparing the stellar population in the inner region with the outer in CMD (Figure 9), detected sources are Table 3. Γ CC is the slope of the completeness corrected MF and Γ NC is the not-corrected one.The stellar masses are estimated using parsec isochrones with DM=18.49,A  = 0.35 and E(H-K)=0.1,at different ages (first column) and mass ranges (second column), for the whole FOV (top), and central r < 3" region (middle), and outer r > 3" region (bottom).N is the number of detected sources and N  is the number of stars fainter than K=17 mag that are corrected for completeness according to their position on the image and their brightness (see Section 4).To check the presence of mass-segregation, we used the minimum spanning tree (MST) algorithm by Allison et al. (2009) comparing the length of the MST connecting massive stars to that connecting randomly selected stars.Mass-segregation ratio (Λ   , Eq. 1 in Allison et al. (2009)), defined as the average random path length of the MST and that of the massive stars.   is the number of selected massive stars and randomly selected stars.Figure 11 shows Λ   calculated for different    subsets.Λ   reaches to a constant value when    is larger than the total  2020) and the black squares represents the simulated magnitudes for these 49 sources (explained in Section 6 Method II).Bottom plots are same as tops but for 818 sources detected in J and K in our first epoch observation (Khorrami et al. 2017).
number of massive stars (m > 10 M ) in the FOV which is about 339 using 1.5 Myr isochrones.The bottom plot in Figure 11 shows the Λ   calculated within different radii of the cluster, centered at R136a1.The numbers in parentheses denote the total number of stars more massive than 10 M which are chosen for    .The colour shows the ratio of number of massive stars ( > 10 ) to the total number of stars (down to 1 M ) within each radius.This plot shows the variation of Λ   locally for different sample of    .Λ   is less than 1.2 within the whole FOV and for    > 350, and it is larger than 1.0 in all the conditions.The values of Λ are point to a moderate but significant (above 2 sigma) degree of mass segregation.This must be taken with caution, considering the spatial distribution of completeness (Figure A1).The fact that, globally, we detect less low mass stars towards the centre of the cluster can have an effect on the mass segregation ratio.If more low mass stars are injected in the centre of the cluster, the size of the N most massive stars would still be the same, but the average size of a group of N stars within the complete sample is expected to decrease.The values of Λ   would also decrease in that case, but quantifying this is out of the scope of this work.
Still, these MF slopes in our study are limited to the resolution of the instrument (55 mas in K) and in future, using higher angular resolution data, we may resolve binaries and low-mass stars which affects the shape of MF.
Using the stellar masses estimated at the age of 1, 1.5, and 2 The colour shows the ratio of number of massive stars ( > 10 ) to the total number of stars within each radius.
Myr, which are corrected12 for the completeness (for 1226 sources), we plot the 2D (projected) mass density at a given radius, i.e. the mass between  and  +  divided by the corresponding area (, Figure 12-top).We used an Elson-Fall-Freeman (EFF) profile Elson et al. (1987) to fit  in the core of R136 (Eq. 1) up to the maximum radius of   = 1.32pc from R136a1.
The total mass of the cluster (  ) and the cluster's surface density (Σ) which is the projected mass density within a given radius, i.e. all the mass between 0 and r divided by the area of the corresponding circle, are shown in Figure 12.Σ and   become consistent at different ages for larger radii (starting at 0.3pc up to 1.3pc), where lots of pre-main sequence and low-mass stars are located/detected.The surface density of R136 at  = 0.3pc (and   = 1.3pc) reaches to 1.4 +0.7 0.4 ×10 4 [M /pc 2 ] (and 2.7 +1.1 0.6 ×10 3 10 Z. Khorrami et al. [M /pc 2 ]), and the total mass down to 1.0±0.1 M is 4.0 +1.7 −1.2 ×10 3 M (and 1.5 +0.5 −0.4 × 10 4 M ), respectively.The densities and total estimated masses at 1 Myr is slightly higher than the values reported in the first epoch (compare Figure 12 to Figure 13. in Khorrami et al. (2017)).
Fitted  and  parameters (in Eq. 1), vary from 0.8 to 0.85 and 0.14 to 0.19, respectively, depending on the selected age.
Comparing these values to the previous ones estimated for 818 stars in the 2015 data in J-K data (Khorrami et al. 2017),  0 and  are consistent within the given errors to their previous values of  0 = (3.89+1.60  −1.14 ) × 10 4 [ / 2 ] and  = 0.17 ± 0.05, but  is decreased by about 0.4.

JHK COLOURS
Among 818 sources detected in J and K in 2015 data, and 1499 sources detected in H and K in 2018, we could detect 790 sources common in total (J, H and K).30% of the newly detected sources in 2018 had enough SNR to be detected in the K data in 2015 (they are brighter than the faintest star in the 2015 K catalogue) but these sources were not listed in the J since the AO correction for J data is not as good as in K.These sources were either located in the central part of the cluster where AO halo (uncorrected part of the stellar signal with the FWHM of seeing) was large leading to decreasing their SNR, or in the outer part of the cluster were the PSF were distorted (see photometric selection criteria in Khorrami et al. (2017)) due to anisoplanatism.
For these 790 common sources between two epochs we have plotted colour-colour diagram in (J -K) vs (H -K) for 413 stars in the central region  < 3" (Figure 13 top) and 377 sources in the outer region  > 3" (Figure 13-bottom).The black circles in Figure 13 shows 49 stars studied by Bestenlehner et al. (2020) in the centre of R136, and the black solid line is the PARSEC isochrone at 2 Myr.These plots show that the detected sources in the inner region of R136 are more consistent with the evolutionary models, than the sources in the outer region of the cluster.Comparing the CMD in H-K and J-K (Figure 9) for inner and outer regions, the dispersion in the colour of the detected sources in the outer region is higher than the inner region.In the first epoch, 48% of the detected sources in J and K were in the outer region, but in the second epoch this number increases to 58%.Only about 30% (70%) of the new detected sources in the second epoch (H-K) are located in the inner (outer) region.This can be explained by the effect of incompleteness in the core, so that even with the longer exposure time and better observing condition (e.g.SR > 70%) central region remains too bright for low-mass MS stars to be detected.The completeness in H and K, for the central region of R136, is less than 50% (20%) for the magnitude range of 19-20 (20-21), in our completeness maps (Figure A1).

CONCLUSIONS
We presented a new photometric analysis of the core of R136 using the second epoch data from VLT/SPHERE instrument in the near-IR.We observed R136 in H and K filters in better atmospheric observing conditions and longer exposure time than the first epoch.This enabled us to detect twice as many sources (in H and K) as we have detected in the first epoch (in J and K), in the FOV of IRDIS (10.8" × 12.1") covering almost 2.7 × 3.0 pc of R136 core.Among 1658 and 2528 detected sources in H and K, respectively, we found Top: projected mass density () profile of R136 in IRDIS FoV centered on R136a1, using 1499 stars in H-K.The stellar masses are estimated at the age of 1 (purple squares), 1.5 (green circles), and 2 Myr (red triangles) with extinction values of   = (0.35 ± 0.5) ,  ( −  ) = (0.1 ± 0.1)  in H and K. Eq. 1 is used to fit the purple, green, and red solid lines to the data at 1, 1.5, and 2 Myr, respectively.Center: surface densities (Σ) within a given radius.Bottom: total stellar masses (  ) within a given radius.Solid blue lines shows the average value of surface density (Σ = 2.7×10 3  /  2 ) and total mass (  = 1.5×10 4  ) up to   = 1.3  1499 common sources between these two sets of data, where 76% of these sources have visual companion closer than 0.2", which is higher than the value found in the first epoch in 2015 data (Khorrami et al. 2017).About 71% of the newly detected sources are located in the outer region ( > 3") of the cluster.
Using the stellar parameters (    , log g, and logL) of 49 stars studied spectroscopically by Bestenlehner et al. (2020) in the optical, PARSEC isochrones, Tlusty SEDs (Hubeny & Lanz 1995;Lanz & Hubeny 2003), and synthetic extinction curves from Draine (2003a,b,c); Li & Draine (2001); Weingartner & Draine (2001) we estimated the extinction at DM of 18.49 magnitude, using two methods (Section 6).We adopt the extinctions of A  = 0.45 and A  = 0.35 magnitude, which are consistent with the two methods (within the error-bars) and De Marchi & Panagia (2014).The colour excess of  ( − ) = 0.1 is also consistent with Tatton et al. (2013) and with the one used in previous studies by Campbell et al. (2010).Consequently, the stellar masses were calculated at the ages of 1.0, 1.5 and 2.0 Myr.The MF slope for 1, 1.5, and 2 Myr isochrone at the inner ( < 3") and outer region ( > 3") of the cluster, are estimated and shown in Table 3.One can check the effect of incompleteness, by comparing the MF slopes before and after the completeness correction, shown as Γ   and Γ  , respectively in Table 3.The completeness-corrected MF slopes for the whole FOV (Γ 1  = −0.93 ± 0.08, Γ 1.5  = −0.98 ± 0.09), are consistent the values derived from the photometric analysis of the first epoch data (Khorrami et al. 2017) for the mass range of (3-300) M (Γ 1  = −0.90± 0.13, Γ 1.5  = −0.98 ± 0.18), and are closer to the Salpeter value (Salpeter 1955) for the high mass range 10-300 M (Γ 1.5  = −1.26± 0.01, Γ 2  = −1.34± 0.03).The MF slopes for the mass range of 10-300 M are about 0.3 dex steeper than the mass range of 3-300 M , for the whole FOV and for different radii.The MF slopes in the inner region are shallower than the outer region for different mass ranges.In Figure 11, Λ   with a value (about 2 sigma) above 1, shows a degree of mass segregation.Considering the spatial distribution of completeness (Figure A1), both MF slope in the core and Λ   would decrease, if number of low-mass stars increases in the centre where the completeness is very low.
We corrected the MF for completeness for sources fainter than 17mag (1226 stars in the whole FOV and 444/782 stars in inner/outer region).Still these values are low limits to the steepness due to incompleteness and central concentration of bright stars.
Comparing data with the first epoch ones, we could detect 790 sources common in total (J, H and K) and the majority (67%) of detected sources in the outer region ( > 3") are not consistent with the evolutionary models at 1 − 2 Myr and with extinction similar to the average cluster value, suggesting an ongoing star formation within 30 Doradus.A significant scatter in the CMD (Figure 9) and colour-colour diagram (Figure 13) is originated mainly from the lower part of the main sequence and pre-main sequence sources, located (detected) at the outer region.The observed scatter in CMD and colour-colour diagrams of 30 Doradus has previously been reported in the visible and NIR (Hunter et al. 1995;Andersen et al. 2009;De Marchi et al. 2011;Cignoni et al. 2015).The scatter is likely due to a combination of observational confusion (affected mainly by the visual multiple systems or variables), photometric errors, differential extinction, and a possible age spread.In addition, pre-main sequence stars are often associated with circumstellar disks and outflows which will introduce additional extinction for the clusters low-mass content.Brandl et al. (1996) found that the extinction varies significantly from star to star within the cluster, with the range of 1-2 mag.The HST observations also reveal the presence of considerable differential extinction across the 30 Doradus region.De Marchi et al. (2011) quantified the total extinction toward massive main sequence stars younger than 3 Myr to be in the range 1.3 <   < 2.1.This motivates us to observe this cluster again in J and K in future with even longer exposure time, so the number of common sources in three filters (J, H, and K) increases.

DATA AVAILABILITY
All data are incorporated into the article and its online supplementary material.A2 for their IRIDS (J, H, K) and HST (U, V, I) magnitudes.Left: whole FOV, Right: same as left but zoomed in the core to avoid confusion.Background image is the SPHERE/IRDIS/K in 2018.
Table A1: List of brightest spectroscopically known stars in the core of R136 studied in optical.T eff , log L/L , and log g, and spectral types are taken from Bestenlehner et al. 2020 Crowther & Dessart (1998), Crowther et al. (2010), bestenlehner et al. (2014), and Crowther et al. (2016), respectively.HSH95 WB85, and ID  are the stars identification in Hunter et al. (1995), Weigelt & Baier (1985), and SPHERE/IRDIS/K catalogue in 2018 (this work), respectively.K and H are the SPHERE IRDIS magnitudes and J 15 is the SPHERE IRDIS magnitude in our first epoch data (Khorrami et al. 2017).r is the distance from r136a1 in our K data.The stellar initial masses (M  ) and ages are estimated by fitting grids of PARSEC evolutionary models (0.1 to 10 Myr) to the T eff and log L/L .

HSH95
ID   (1985), followed by their U, V, and I magnitudes in the HST/WFPC2 filters.Last column is V-K (F555W-K( 2018)).See Figure A2 for their positions in the image.This paper has been typeset from a T E X/L A T E X file prepared by the author.

Figure 2 .
Figure 2. Comparison of R136 central images at different wavelengths with the highest available angular resolution telescopes.Top, left to right: HST/WFPC2 in V , VLT/MAD in K, and VLT/SPHERE in K. Bottom, left to right: HST/NICMOS/NIC1, VLT/MAD, and VLT/SPHERE in the H.The identification of bright sources is shown in Figure A2.

Figure 3 .
Figure 3. Signal to noise ratio of detected sources in H (purple circles) and K (green squares) versus their magnitudes.Filled dots are common sources between H and K data.The solid horizontal line shows the SNR=2.0.

Figure 4 .
Figure 4. Histogram of the separation of the close detected sources versus their distance from the core of R136.Top: 818 common sources between J and K data in 2015.Bottom: 1499 common sources between H and K data in 2018.

Figure 5 .
Figure 5. Average value of completeness values for the whole FOV as a function of magnitude using 62,000 artificial stars in H (purple dots) and K (green squares).

Figure 6 .Figure 7 .Figure 8 .
Figure 6.Completeness map for magnitude range of 19 -20 in H (top) and in K (bottom).Black pluses shows the position of photometric detected sources within the given magnitude range.colour shows completeness in percents.

Age
60 ± 0.05 −1.59 ± 0.04 compatible with the evolutionary models in the inner region (Figure9-top middle), but for the outer region detected sources differ from the isochrones.About 50% of the high mass stars are located in the right side of the main-sequence, far from the evolutionary models.If, these population truly belong to R136 then they have K bandpass excess emission due to hot dust (K-excess sources).The MF slopes in the outer region are complete for the both mass ranges and the slope values are steeper than the inner region, but might be the effect of the existence of these K-excess sources.The MF slopes for the massive part (M > 10 M ) in the inner region is in line withBestenlehner et al. (2020) ( ∼ 2.0 ± 0.3) and for the whole FOV is consistent withSchneider et al. (2018) ( ∼ 1.90+0.37−0.26 ) for the 30 Doradus region, considering the large error bars from these studies.Massey & Hunter (1998) found the MF slope of −1.4 < Γ < −1.3 for the mass range of 15-120 M which is consistent with the slope for the whole FOV at 2 Myr.

Figure 9 .
Figure9.Top: CMD of detected sources common between the H and K data for the full FOV in left (N=1499), for central  < 3" region (N=627) in middle, and for outer  > 3" region (N=872) in right.Red, yellow and blue solid lines are PARSEC isochrones at DM=18.49 with   = 0.35 and  ( −  ) = 0.1, at the ages of 1, 1.5, and 2 Myr, respectively.The black circles are the 49 spectroscopically-known stars fromBestenlehner et al. (2020) and the black squares represents the simulated magnitudes for these 49 sources (explained in Section 6 Method II).Bottom plots are same as tops but for 818 sources detected in J and K in our first epoch observation(Khorrami et al. 2017).

Figure 13 .
Figure 13.colour-colour diagram for the common sources between J, H and K, within two epochs J-K (2015) and H-K (2018), for the central region  < 3" (top) and outer region  > 3" (bottom).The colour indicates the K magnitude of detected sources and the black circles represent 49 stars studied by Bestenlehner et al 2020 listed in TableA1.The black solid line is the PARSEC isochrone at 2 Myr with  ( −  ) = 0.1 and  (  −  ) = 0.25.The right hand side curve represents pre-main sequence part of isochrone followed by the main sequence locus on the straight line on the centre and the post-main sequence part on the top.

Figure A1 .Figure A2 .
Figure A1.Completeness maps in K (top) and H (bottom), used for correcting stellar masses for plotting MF. colour range [0:100].Stars shows the location of the observed detected sources within the given magnitude range for each map.If a star with a given magnitude is located in the yellow area, it will be detected 100% and it has difficulty to be detected if it is located in the dark-purple.

Table 1 .
Table 1 shows the exposure time log of the two epoch observations Exposure time log of VLT/SPHERE observations on R136.

Table A2 :
Hunter et al. (1995)rces in the centre of R136 brighter than 16 mag in K. ID  is the stars identification in our K catalogue.K(2018) and H(2018) are the K and H magnitude of stars in the second epoch data in 2018.K(2015) and J(2015) are the K and J magnitude of stars in the first epoch data in 2015.HSH95 and WB85 are the stars identification inHunter et al. (1995)andWeigelt & Baier

Table A3 :
Sample catalogue of detected sources in H and K data in 2018.Full table is available online.ID  is the stars identification in our K-band catalogue.K and H are the K and H magnitude of stars in 2018.Correlation  and Correlation  are the correlation between the PSF of detected source with the reference PSF in H and K images. is the photometric error estimated for the positions and magnitudes.