Time-domain Study of the Young Massive Cluster Westerlund 2 with the Hubble Space Telescope. I

We repeatedly observed the YMC Wd2 from 2016 to 2019 October for a total of 47 orbits of HST, using the UVIS channel of WFC3 (PI Sabbi, GOs-14087, 15362, and 15514). Each orbit consists of two 3 s (hereafter short exposures) and six 350 s (hereafter long exposures) observations in the filter F814W. For the long exposures, we adopted a six-point dither pattern with a step of 126 to guarantee an optimal sampling of the point-spread function (PSF), fill the gap between the UVIS detectors, and allow for hot pixel removal. The two short exposures were acquired at the beginning and at the end of each orbit, respectively, at the same position of the first and last long exposures. To identify and characterize the various types of Wd2 variable stars (periodic, semiperiodic, and nonperiodic), we organized the observations in groups of 5 or 10 visits using a nonuniform, logarithmic increasing spacing.

In this paper we limit our analysis to the long exposures, acquired during the first 26 visits, taken between 2016 October 30th and 2018 July 5th. Over this period we collected a total of 156 long exposures organized as follows:

• 1.  
  The first five epochs were separated by 51 days each, and were acquired between 2017 October and June.
• 2.  
  In 2017 August we observed Wd2 10 times with a 0.8 day cadence.
• 3.  
  And 1.6 days later we obtained the other five epochs, each separated by 1.6 days.
• 4.  
  The first epoch belonging to the 102 days cadence was observed on 2017 November 27th.
• 5.  
  Lastly, five epochs separated by 3.2 days were acquired between 2018 March 9th and July 5th.

In Zeidler et al. (2015) we showed that Wd2 consists of two nearly coeval clumps. To guarantee the best coverage of the entire cluster population, we centered the gap among the two UVIS CCDs between the two clumps, in a region almost devoid of stars. To allow for a more flexible scheduling of the observations, we did not request any specific orientation. The left panel of Figure 1 shows the stack image, obtained by combining all deep exposures, while the right panel shows the coverage map of the survey, and highlights how many images contributed to each pixel.

Zoom In Zoom Out Reset image size

We carried out the photometric analysis directly on the bias-subtracted, flat-fielded, and charge-transfer-efficiency corrected _flc.fits exposures produced by the standard calibration pipeline calwf3 v3.5.⁹ Compared to drizzled (_drc) images, _flc images have the advantage of not being resampled and therefore they are the most direct representation of the astronomical scene. However _flc images are still affected by geometric distortion, which we took into account by creating a distortion-free reference frame and relating the photometry and astrometry of each exposure to that frame.

We used the publicly available one-pass photometry routine img2xym (Anderson & King 2006) and a library of empirical PSFs to identify and measure all the high signal-to-noise ratio (S/N > 150), isolated (within 5 pixel radius), and nonsaturated sources. The PSF library properly accounts for the spatial variations caused by the optics of the telescope and the variable charge diffusion in the CCD. Seasonal and daily thermal changes, however, can cause variations in the optical path length of HST up to a few microns within the timescales of an orbit (Bély et al. 1993; Lallo et al. 2006). These changes affect the focus of the telescope and translate to small, but measurable differences in the PSF from one exposure to another. To take these variations into account, for each exposure we derived the correction necessary to minimize its average PSF-residuals. We then derived for each catalog a six-parameter transformation to align the UVIS observations to the Gaia Data Release 2 (DR2, Lindegren et al. 2018), and to create a distortion-free reference frame with a pixel scale of 40 mas/pixel, the UVIS native pixel scale.

To separate candidate members of Wd2 from the field of the MW, and to estimate masses and temperatures of the stars in the cluster, we combined our data with our previous WFC3/IR observations acquired in the filters F125W and F160W (PI Nota, GO-13038). We aligned the IR observations to Gaia using the same procedure followed for UVIS.

The final photometric analysis was performed using the software KS2 (J. Anderson 2020, in preparation; see Sabbi et al. 2013, 2016, for details), that can use all the exposures together to detect, measure, and subtract stars in subsequent passes, moving from the brightest to the fainter sources. KS2 can perform simultaneous detection in multiple filters. However, given the higher spatial resolution of the UVIS data set, we only used the long exposures acquired through the UVIS/F814W filter to create the reference catalogs. For each source KS2 provides the mean flux, the flux standard deviation, and an estimate of the quality of the fit of the PSF. The final catalog contains in total 9268 sources.

The fluxes of the sources found by KS2 in the UVIS/F814W combined observations were measured in the WFC3/IR exposures, and in each of the individual F814W visits. For each exposure we used its own "focus-specific" PSF library, derived as described above.

Fluxes were converted into instrumental magnitudes and then calibrated in the Vega-mag photometric system following the same approach described in Sabbi et al. (2016), and using the zero-points listed on the STScI website.¹⁰ The aperture correction was derived using a 04 aperture photometry on the drizzled images.

The LCs of WTTS are thought to be dominated by the asymmetric distribution of cool or dark magnetic spots on the stellar surface, with stellar rotation being responsible for the modulation of the fluxes (Bouvier et al. 1993; Herbst et al. 1994). These sources show minimal to no evidence of ongoing accretion, possibly because the innermost part of their circumstellar disks have already been depleted.

Three percent of Wd2 stars fainter than m_F814W = 17.0 show periodic (with periods ranging between 1.6 and 46 days, and with an average value of 7.7 days), small-amplitude (0.035 × Δm_F814W ≤ 0.75 mag) fluctuations consistent with Class iii WTTSs. The position of the WTTS candidates in the m_F814 versus ${m}_{{\rm{F}}814{\rm{W}}}-{m}_{{\rm{F}}160{\rm{W}}}$ CMD is shown in panel (A) of Figure 7. The number of WTTS candidates considerably decreases below m_F814W ≃ 21.5, with the faintest source found at ${m}_{{\rm{F}}814{\rm{W}}}\,\lesssim 23.8$. This drop is likely caused by incompleteness effects. In a future paper, we will examine all the 47 orbits of HST time and perform artificial star tests to evaluate the impact of various selection biases, including temporal sampling, S/N, and crowding, on the five populations of variables. Panel (A) in Figure 8 shows that the WTTS candidates have relatively blue colors, in agreement with the hypothesis that these objects are likely quite evolved PMS stars.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

We used the Stefan–Boltzmann law to convert the luminosities and temperatures derived from Padova isochrones into stellar radii (1.5 ≤ R_WTTS ≤ 5.7R_⊙). In doing this we found that the rotational velocities v sin(i) of the WTTS candidates range between 2.3 and ∼168 km s⁻¹, with a peak at ∼12.5 km s⁻¹, in agreement with the v sin(i) measurements found in the literature (i.e., Clarke & Bouvier 2000; Herbst et al. 2002; Xing et al. 2006).

CTTSs are class ii young stellar objects that are still actively accreting material from their circumstellar disks. We identified 757 variable sources (∼15% of Wd2 low-mass stars) whose LCs are highly variable, but do not show evidence for periodicity on a short timescale. These objects cover the magnitude range $17.0\lesssim {m}_{{\rm{F}}814{\rm{W}}}\lesssim 24.2$ (panel (B) of Figure 7), and, compared to WTTS candidates, extend toward redder colors both in ${m}_{{\rm{F}}814{\rm{W}}}-{m}_{{\rm{F}}125{\rm{W}}}$ and in ${m}_{{\rm{F}}814{\rm{W}}}-{m}_{{\rm{F}}160{\rm{W}}}$ (Figure 8).

A special class of CTTS was discovered about 10 yr ago. These objects are characterized by largely flat LCs, interrupted by sharp and short dips, and they are therefore often referred to as "dippers" (Alencar et al. 2010; Morales-Calderón et al. 2011; Cody et al. 2014; Ansdell et al. 2016; Hedges et al. 2018). Because the dips are normally shallower in the IR than at optical wavelengths, they have been attributed to large and dusty structures (such planetesimals) passing along our line of sight (Bouvier et al. 2003; Alencar et al. 2010; McGinnis et al. 2015), although their precise location within the disk (inner versus outer disk, Bouvier et al. 1999; Morales-Calderón et al. 2011; Bodman et al. 2017) is still debated.

Seventy-two of the sources in Wd2 between $19.14\,\leqslant {m}_{{\rm{F}}814{\rm{W}}}\leqslant 23.79$ (0.3 ≤ M ≤ 2.4M_⊙) show relatively flat LCs, interrupted by sudden drops in luminosity as large as 2.78 mag in the m_F814W filter. Based on our temporal coverage, we conclude that the observed drops in luminosity are short events, that last a few tens of days at most. An example of this type of object is shown in Figure 9, where the luminosity of the star #5320 (R.A. $=\,{10}^{{\rm{h}}}{24}^{{\rm{m}}}05\buildrel{\rm{s}}\over{.} 463$, decl. $=\,-57^\circ 45^{\prime} 08\buildrel{\prime\prime}\over{.} 838$, J2000) in the filter F814W dropped by ∼2.0 mag two times in less than two years.

Zoom In Zoom Out Reset image size

In Wd2 about 8% of the stars fainter than ${m}_{{\rm{F}}814{\rm{W}}}\lesssim 16.62$ are characterized by eruptive behaviors, with the largest increase in magnitude we observed so far being Δm_F814W = 3.22. In some cases the variation lasts for only a few days, in others, once the star lights up, its luminosity remains elevated.

Eruptive low-mass PMS stars are common and, based on the length of their burst, have been classified as either EX Orionis-type (short burst Herbig 2007) or FU Orionis-type (few year long burst, Hartmann & Kenyon 1996; Hartmann et al. 1998) stars. It is not clear yet if EXor- and FUor-type stars belong to the same class of objects, or if all the young low-mass stars undergo a bursting period. Nevertheless for both EXor- and FUor-type stars, the large variation in brightness is ascribed to unstable pile-up of gas in the inner disk, which then releases a cascade of material onto the star. It is believed that most of the stellar mass is accumulated via these dramatic outbursts (e.g., Bell & Lin 1994; Zhu et al. 2009; D'Angelo & Spruit 2010), during which a star can accrete up to 0.01M_⊙ of gas.

WISE observations of "burster" stars, identified during the K2 mission, show clear IR excess, typical of large inner circumstellar disks, and the strength of the outburst appears to correlate to the IR excess (Cody et al. 2017). These bursts are likely "cooking" the dust in the circumstellar disks, and thus modifying their chemical composition (Green et al. 2006; Quanz et al. 2007; Cieza et al. 2016). Understanding the properties and the frequency of these bursting systems is therefore important for the models of stellar mass accretion, protoplanetary gas-rich disk evolution, and, possibly, planet formation.

The phase-folded LCs of nearly 2% of the objects between $0.1\lesssim M\lesssim 5.8$ M_⊙ (panel (E) in Figure 7) show two minima, characteristic of eclipsing binaries. For each of these objects we derived physical and orbital properties by iteratively fitting their phase-folded LCs with the publicly available software Nightfall.¹¹

We assumed the orbital period and the ephemeris zero-point derived by the Lomb–Scargle statistics, and the total mass inferred from Padova isochrones. Nightfall then compared the observed LC with simulations obtained for different mass ratios, stellar temperatures, inclination angles, and Roche lobe filling factors using a chi-square function to determine the best-fit solution. The derived values depend on the assumptions made in modeling the stellar atmospheres in the Nightfall software and in the Padova isochrones, used to infer the total masses of the binaries. Therefore they should be used only for comparisons and not as absolute values. Because of the sparse cadence and the relatively short temporal coverage of our observations, in this preliminary study, we assumed only circular orbits, and, in the case of interacting binaries, in the LC model we did not include the contribution of a mass-transferring disk.

Because of the limited temporal sampling, our analysis is biased toward the shortest periods. Therefore we were able to find only eclipsing binaries with orbital periods ranging between 1.9 ≤ P ≤ 86 days, and separations between 6.8 ≤ R ≤ 108 R_⊙ (0.03-0.5 au). The upper panel of Figure 10 shows the distribution of periods covered in four different bins of system masses. The lower panel shows the distribution of mass ratios (Q = M₂/M₁) for the same bins of mass. It is interesting to note that in all four bins the distribution of the Q parameter is bimodal, with only ∼one-third of the binaries consisting of nearly equal-mass objects, similar to what observed by Raghavan et al. (2010) for binary systems, with primary star having mass M₁ ∼ 1 M_⊙.

Zoom In Zoom Out Reset image size

In 50% of the systems the LCs are consistent with those of semidetached binaries where one of the two components (the primary in 40% of the cases and the secondary in the remaining 10%) has filled its Roche lobe. In 2% of the cases the LCs are consistent with those of contact binaries, with both stars having filled the Roche lobe. Because of the considerable fraction of interacting binaries, with ongoing mass transfer, it is likely that the mass ratio distribution will change over time, increasing the number of equal-mass systems. The distribution of the filling factors for both the primary (in red) and the secondary stars (in orange) for four different bins of masses is shown in Figure 11.

Zoom In Zoom Out Reset image size

Figure 12 shows the mass to radius relation for the primary (blue circles) and the secondary (orange triangle) star as derived from the fitting of the EB LCs. The radii of the primary components range between 0.1 ≤ R_p ≤ 33.8R_⊙ and for the secondary between 0.1 ≤ R_s ≤ 13.6. In many cases these values are considerably larger than those derived from Padova isochrones using the Stephan–Boltzmann law. Phenomena such as tidal heating, ohmic dissipation and thermal tides have been considered to explain the inflated radii observed in close binaries and hot Jupiters; however, in this case the fact that the library of models of atmosphere in Nightfall is not optimized to describe PMS stars and that the objects are quite far from the approximation of blackbodies is likely responsible for at least part of the discrepancies.

Zoom In Zoom Out Reset image size

Figure 13 shows the spatial distribution of the five populations of variable PMS stars projected on a image of Wd2 taken through the filter F814W. Isodensity contours for the stars belonging to Wd2 are shown for reference in all five panels.

Zoom In Zoom Out Reset image size

A visual inspection of Figure 13 suggests that WTTS candidates (panel (A)) are more concentrated around the two main clumps, while CTTS candidates (panel (B)) and bursters (panel (D)) extend toward larger radii. Given that WTTS candidates are relatively massive objects, this is likely only a consequence of the fact that Wd2 is mass segregated (Zeidler et al. 2015).

In relatively low-density star-forming regions, like NGC 2264, dippers represent about 20%–30% of the CTTS population (Alencar et al. 2010; Cody et al. 2017; Cody & Hillenbrand 2018) and from the analysis of K2 mission data Hedges et al. (2018) found that dippers and bursters should have similar spatial distributions. However, panel (C) shows that in Wd2 dippers clearly avoid the two regions of highest density, where the majority of the high-mass stars (up to 80 M_⊙ Bonanos et al. 2004) are concentrated. In particular, based on the number of bursters and CTTS found within 150 pixels (∼1.12 pc) from the main clump of Wd2, in the same region we should find at least five dippers if they had the same spatial distribution of the other variable stars.

We used Ripley's K function (Ripley 1977, 1979) to further investigate the clustering properties of the variable stars in Wd2. According to Ripley's statistics, a population is considered clustered within a certain scale r if K(r) is above the response of the K function for a Poisson distribution. On the contrary, distributions whose K function at the distance r are below the Poisson distribution are considered dispersed. Figure 14 compares the K functions for WTTS candidates (in blue), CTTS candidates (in orange), dippers (in green), bursters (in purple), and EBs (in red) to the K function of a Poisson distribution (black dashed line) over the entire range of magnitudes considered (panel (A)), for stars between 17 < m_F814W ≤ 20.0 (corresponding to the mass range between 1.75 and 3.3 M_⊙ when using Padova isochrones, panel (B)), between 20 < m_F814W ≤ 22.0 (corresponding to the mass range between 0.8 and 1.8 M_⊙, panel (C)), and between 22 < m_F814W < 24.5 (corresponding to the mass range between 0.2 and 0.8 M_⊙, panel (D)).

Zoom In Zoom Out Reset image size

Over the entire range of magnitudes (17 < m_F814W < 24.5, Figure 14 panel (A)) EBS, CTTS candidates, and bursters have the same clustering scale r = 1.41 ± 0.14 pc and their K functions are indistinguishable from each other. On the same range of magnitudes, the clustering properties of WTTS candidates are less significant than those of the previous three populations, but their clustering scale is comparable (r_WTTS = 1.32 ± 0.14 pc).

The clustering properties of the dippers between 17 < m_F814W < 24.5 are only marginally significant, and in the innermost 0.08 ± 0.01 pc they are dispersed. Their clustering scale is also smaller than the other populations, being only r_Dip = 1.05 ± 0.20 pc.

Because the cluster is mass segregated (Zeidler et al. 2017), we checked if Ripley's statistics could give different results over different ranges of magnitude/masses. If we consider only the variable PMS stars brighter than m_F814W ≤ 20.0, we find that EBs are significantly more clustered than CTTS and WTTS candidates. In this range of magnitudes CTTS and WTTS candidates have identical distributions, while bursters are almost indistinguishable from a Poisson distribution. Only four dippers are brighter than m_F814W ≤ 20.0, and therefore we do not considered their distribution in this range of masses.

Between 20.0 < m_F814W ≤ 22.0, bursters have clustering properties and significance similar to CTTS candidates, while the properties of WTTS candidates resemble those of the EBs. On scales smaller than r ≃ 800 pixels (∼0.7 pc) dippers' properties are similar to those of EBs and WTTS candidates, but at larger scales they become more similar to bursters and CTTS candidates. Finally for magnitudes fainter than m_F814W > 22.0 only CTTS candidates and bursters are significantly clustered.

The lack of dippers in the regions of highest stellar density suggests that the properties of stellar disks vary with the position in the cluster, and that in less than 2 Myr phenomena such as dynamical truncation (Vincke et al. 2015; Portegies Zwart 2016; Vincke & Pfalzner 2016), stellar winds (Pelupessy & Portegies Zwart 2012), dynamical interaction (Olczak et al. 2008; Reche et al. 2009; de Juan Ovelar et al. 2012), and photoevaporation caused by the UV radiation coming from nearby bright OB stars (e.g., Haworth et al. 2017), can significantly affect disk growth and evolution.

Similarly the different behavior at different masses of the various types of variables can provide information on the life span and size of the circumstellar disks as a function of the mass of the star, although these trends will have to be further investigated to properly evaluate the effect of incompleteness.