Free Access
Issue
A&A
Volume 619, November 2018
Article Number A41
Number of page(s) 46
Section Stellar structure and evolution
DOI https://doi.org/10.1051/0004-6361/201833667
Published online 05 November 2018

© ESO 2018

1. Introduction

The paradigm of accretion in young stellar objects (YSO) has shifted from a model of constant mean accretion rate to that favouring short events of intense accretion (Vorobyov & Basu 2006, 2015; Zhu et al. 2009). This shift is largely to address the issue of the “protostellar luminosity problem” (Kenyon et al. 1990; Kenyon & Hartmann 1995; Dunham et al. 2014). A variety of models including turbulent or competitive accretion, accretion regulated by core, disk, and feedback, are invoked to understand the deviation from the idealized case of isothermal sphere (Kenyon et al. 1990; McKee & Offner 2010; Myers 2010; Vorobyov & Basu 2008; Dunham & Vorobyov 2012; Dunham et al. 2014 and references therein). However, most of these models share the variable accretion component, albeit differing at various mass regimes. The accumulated observational evidence appears to favour variable accretion instead of constant mean scenarios (Dunham et al. 2014). Photometric variability of YSOs can be related to their natal environment, accretion physics or a combination of both (Contreras Peña et al. 2017a; Kesseli et al. 2016; Meyer et al. 2017 and references therein). Some of the variability can be caused by cold and hot spots formed on the surface of the YSO by infalling material from the disc. Dust clumps in the stellar medium surrounding the YSO can cause variable extinction of star-light as it passes along the observers line of sight (e.g. Herbst & Shevchenko 1999; Eiroa et al. 2002 among others).

The FUors (FU Orionis) and EXors (EX Lupi) examples of high amplitude photometric variability result from variable accretion. Respectively, they last from a few years to a few months. These objects are known to be low-mass YSOs, although similar counterparts in the higher mass range have been found (Kumar et al. 2016; Caratti o Garatti et al. 2017). Kumar et al. (2016) uses highly variable light curves (LCs) of massive young stellar objects (MYSOs) candidates from the Vista Variables in the Via Lactea (VVV) survey (Minniti et al. 2010), arguing that they were signposts of ongoing episodic accretion. Photometric and spectroscopic variability in a 20 M MYSO was used by Caratti o Garatti et al. (2017) to conclude that disk-mediated accretion bursts are a common mechanism across stellar masses. ALMA observations were used by Hunter et al. (2017) as evidence that sudden accretion is responsible for the growth of a massive protostar. These findings suggest that episodic accretion maybe a common mechanism in star formation, independent of mass. Computational models predict luminous flares in MYSOs, which are morphologically similar to FUors and EXors (Meyer et al. 2017).

The findings in Kumar et al. (2016) raise the question of the overall nature of variability in massive YSOs. In this paper, we attempt to examine the variability phenomena in known extended green objects (EGOs; Cyganowski et al. 2008) and IR sources, deeply embedded in clumps identified by the APEX Telescope Large Area Survey of the Galaxy (ATLASGAL; Schuller et al. 2009). They represent unbiased large samples of point-like massive young stellar candidates, therefore, allowing us to use the point source photometry to examine variability. We surmise that the RMS and UCHII regions represent an important MYSO sample, however, it requires larger aperture photometry of extended objects to examine variability, which we postpone to a different study.

Employing point source photometry requires that the targets are point like in MIPS, have associated high mass star forming signposts, and finally they are also point-like in the K s band. The selection of such targets is described in Sect. 2. In Sect. 3 we describe the results obtained and discuss their implications in Sect. 4.

2. Target sample, data, and methods

Identification of point-like MYSO targets is based on the Spitzer GLIMPSE and MIPSGAL surveys (Carey et al. 2009), the ATLASGAL survey (Schuller et al. 2009), and the VVV survey (Minniti et al. 2010). These different surveys are highly complementary, covering much of the same area but at different wavelengths (from ∼1.2 to 870 μm). We searched for: (a) driving sources of EGOs (Cyganowski et al. 2008; Chen et al. 2013a,b) and (b) luminous MIPS 24 μm point sources embedded in ATLASGAL clumps. The two samples are expected to roughly represent two early evolutionary phases of massive stars; the EGOs, with an active phase of mass ejection, and non-EGOs which are likely yet to begin outflow activity.

2.1. MYSO sample

2.1.1. EGO sample

EGOs are objects with extended emission in the Spitzer 4.5 μm band (IRAC 2). This band is of particular significance since it contains both H2, and CO lines, which can be excited by shocks when outflows and jets interact with the interstellar medium (ISM). This is particularly the case when the extended emission in the 4.5 μm band is in excess with respect to emission in the other IRAC bands.

They were first catalogued by Cyganowski et al. (2008), and later the catalogue was extended by Chen et al. (2013a,b). EGOs are thought to represent the H2 flows driven by MYSOs (Cyganowski et al. 2008) or MYSO outflow cavities (Takami et al. 2012). A total of 270 unique EGO targets have been catalogued so far. By original classification (Cyganowski et al. 2008) these targets have a MIPS 24 μm detection, usually representing the driving source of the outflow. In order to find the near-infrared counterpart of these driving sources we searched for 2 μm sources in the VVV catalogue with a search radius of 1.0″ and 0.5″ from the known EGO positions. We find 187 and 153 driving sources. We allowed for sources classified as both point-like and extended to be selected, even though, 80% of the detected sources were point-like. Young stellar objects with disk and outflow activity are often surrounded by circumstellar nebulae in the near-infrared, leading to a classification as extended. These objects were kept in the sample list. Additionally, three colour composites, shown in Fig. 2, were used to visually examine whether the identified point sources are good representations of an outflow driving source. This examination led us to retain the 153 sources which clearly represent 2 μm counterparts of the 24 μm source, hence the putative driving source of the EGO target.

2.1.2. Non-EGO sample

Kumar et al. (2016) identified a sample of highly variable VVV objects and found MYSO counterparts in ATLASGAL clumps (Contreras et al. 2013). Here an inverse approach is used. Using ATLASGAL, Contreras et al. (2013) and Urquhart et al. (2014) built the Compact Source Catalogue (CSC) which identified ∼10 000 dense clumps. The mass, density, and distance to these clumps are provided by Urquhart et al. (2018) and they are believed to represent active sites of high-mass star formation. Assuming that ATLASGAL clumps host MYSOs, we searched for MIPSGAL point-like sources that matched with ATLASGAL CSC sources within a radius of 5″. This ensured that we matched red point-like sources in 24 μm band with the peak emission in the 870 μm observations of ATLASGAL. 873 point sources were found with this search. When there were multiple matches we chose the object with the closest centroid distance. The MIPS FWHM is equal to 6″, therefore, a further search of the 873 targets with a matching radius of 5″ was performed with the VVV catalogue, allowing us to find 574 Ks-band targets. These 574 targets display more than one Ks band source within the 5″ radius. In the next step, the point source closest to the MIPS peak was searched for, by constraining the search radius to 1.0″, only for the 574 targets. This retrieved a list of 2171 sources from the VVV source catalogue. Out of these, any source which had less than ten non-saturated epochs (over the full five-year period) was removed. This led us to find 367 single detections and 147 multiple detections in the centroid search with r ≤ 1″. The multiple detection targets were examined visually considering the source magnitude, colour and centroid distance, based on which 66 of the 147 sources were rejected, retaining 81 sources. These 448 (367 + 81) sources are, therefore, the Ks-band point sources representing the MYSO candidate at the peak of an ATLASGAL clump with a spectral energy distribution (SED) that can be assembled from at least 2 μm unto 870 μm.

The final MYSO sample we produced to study the variability is, therefore, composed of 153 EGO and 448 non-EGO sources, resulting in 601 targets. We note that 66 of the 153 EGO targets also lie within the ATLASGAL clumps, the non-EGO sources being exclusively those that coincide with the peak of ATLASGAL clumps.

2.2. VVV survey data

The VISTA Variables in the Via Lactea (VVV) survey has obtained photometric observations in the near-infrared (NIR) passbands (0.9–2.5 μm), covering multiple epochs, spread over five years (from 2010 to 2014) and covering a 520 deg2 area of the inner Galactic plane (see Fig. 1) (Minniti et al. 2010). The survey data is made publicly available through the Cambridge Survey Astronomical Unit (CASU), which is the photometry obtained on the final combined “tile” images. A tile image incorporates multiple “pawprint” that are single exposures on sky. The pawprint data (available to the VVV team and also made public at the VISTA Science Archive in Edinburgh), is the basic product of the observations which often holds better photometric and seeing information, as they are better calibrated and tend to have sharper image profiles than the tile data. In this work we have exploited this full potential by using the pawprint photometry.

thumbnail Fig. 1.

VVV survey area.

Open with DEXTER

2.3. Processing of the pawprint photometry

The pawprint photometry and photometric classification used are standard pipeline products from the Cambridge Astronomical Survey Unit (CASU), as detailed in Lewis et al. (2010). Matching and combining detections between multiple pawprints were made following the approach detailed by Smith et al. (2018). Sources are classified according to their morphology and flagged as 1, 0, –1, –2, –3, –7, –9, respectively as; a galaxy, noise, stellar, probably stellar, probable galaxy, bad pixel within 2″ aperture, and saturated.

The pawprint observing pattern of dithers and overlaps, implies that each source might have between two and six image frames for the same observing epoch, they can also be detected only once along certain edges of the survey. Since these observations are close in time, we chose to bin them together and compute the median magnitude of all observations in intervals of half a day. This binning prevents the detection of variations with timescales smaller than half a day, but reduces the level of scatter in short-periods. The gain in photometric sensitivity will thus be a factor of the number of observations binned, and scales according to where n is the number of binned observations. The typical error in the photometry will be K_err ≤ 0.05 mag which allows for the detection of low-level variability. Following the reasoning explained in Smith et al. (2018) we employ the aperMag2 (r ∼ 0.71″) as the Ks magnitude for all analysis.

For each source we assembled a database that contains: a unique identification, median co-ordinates in the ICRS, median magnitude (over half a day) in the K-band, the median absolute deviation (MAD), the standard deviation, the inter-quartile range (IQR), the number of pawprints in which the source was observed, the number of total observed epochs, the modal class, the number of epochs classified with each flag, the K-band magnitude, the quality classifier of the photometry, and the modified julian date (MJD) of the observation.

The median of the co-ordinates and magnitude, and the modal class were computed for all pawprint observations. The MAD and IQR were computed for each source as they are robust statistical indicators to measure the amplitude and dispersion of the variability (Hampel 1974; Upton & Cook 1996; Sokolovsky et al. 2017). These parameters are less sensitive to outliers than the standard deviation. A high value of the MAD or IQR can be a good indicator of the inherent variability of the source. The IQR measures the amplitude of the difference between the third and first quartiles (Q3 and Q1) of a distribution, in this case the distribution of magnitudes

(1)

The MAD is the median of the absolute differences between each data point and the median, as shown by the following equation:

(2)

in which, Ki is an observation and K represents all the observations.

2.4. Light curves and their reliability

The light-curves (LCs) of 448 non-EGO and 153 EGO targets were produced in the following way by querying the database assembled above. First, we queried a set of co-ordinates and search radius on the database. Secondly, we built a list with all observations that matched the query. Thirdly, we excluded all saturated observations (modal class = –9). Fourthly, we produced a LC for each target using the difference from median (KmedianKmagi).

In the top panel of Fig. 2 we show an example of a LC. To ascertain that the photometry of the target at a given time epoch is not affected by poor observing conditions, flat field errors, improper photometry, or poor seeing, we performed a few tests.

thumbnail Fig. 2.

LC of an eruptive event. Top panel: the LC of the source, the error bars represent MAD(ΔSimjd ). Bottom plot: the RGB image of the source using the Spitzer IRAC 3.6 μm, IRAC 4.0 μm, and the 24 μm MIPS band as blue, green and red, respectively. The VVV source is indicated by the blue circle and the green cross represents the MIPS co-ordinates. The contours of the RGB are in the interval of [Peak-5σ, Peak] from the ATLASGAL observation at 850 μm.

Open with DEXTER

2.4.1. Identifying the variable source

Stellar sources within two annuli defined by r = 1″ and r = 60″ from the target were selected. Typically 100–200 sources were found by this selection. For each such source Si, the magnitude deviation (ΔSimjd) from its median value () over all epochs was computed. The median value for all sources in the annulus, is a representation of the photometric deviation (if any) of the the individual epoch over the time-line. for each source, at each epoch, the offset was added to Simjd to produce the corrected light curve. The MAD of the deviations for all selected sources, MAD(ΔSimjd), is used as an approximation of the 1σ photometric error of a 1′ field around the target for a given epoch. And is shown as the error bars in Fig. 2.

2.4.2. Influence of magnitude on variability

Next we assessed the influence of using any and all, or, only sources with magnitudes comparable to that of the target inside the annuli for computing the 1σ error. For this purpose we filtered the sources within a magnitude range ±1 mag to that of the target, which decreased the number of sources by a factor of approximately ten. Figure 3 illustrates the results of this test. It can be seen that the difference in 1σ error by using the two comparison samples is Ks ∼ 0.0018 − 0.0031 mag, which is 1–2 orders of magnitude below the typical 1σ errors in the target fields.

thumbnail Fig. 3.

Systematic errors as a function of the magnitude of each target. The red points represents the case where we only consider the stellar sources with ±1 mag around our targets. The blue points represent the case in which we consider all stellar sources in the vicinity of our targets.

Open with DEXTER

2.4.3. Control field test

The targets of study are found in the midst of star forming regions, often deeply embedded in dark clouds, leading to reduced and non-uniform source distribution. Also YSOs (in general) are known to be variable objects, so, many sources in a given field may be variable. To address the influence of these effects, we used a control field region randomly selected to be 5′ away from the target field. For each control field the steps explained in the two subsections above were executed. We find that the control field variability was very similar to the MAD(ΔSimjd) computed above.

2.5. Periodograms, false alarm probability, and their aliases

Once the LCs are assembled, we computed the Lomb Scargle periodogram, identify the max power frequency component, and use it to produce a phase-folded LC. Scargle (1982) defined the false alarm probability (FAP) as a measurement of the probability of a signal without any periodic component to have a peak amplitude. The predictive power of the FAP decreases in the presence of correlated noise, non-Gaussian errors, and, highly non-sinusoidal variability. The 90%, 95%, and 99% FAP levels have been computed for periodograms of each target.

A given periodicity can, by a compound effect of binning, observational window, and noise, produce harmonics of itself, which appear in the periodograms as additional peaks, or aliases (VanderPlas 2018). In an effort to verify if the peaks determined were in fact real signals or their aliases, an additional verification step was added. The highest peak of the periodograms and the following 10 highest peaks were identified. Aliases were searched by examining: (a) multiples in the frequency range; (b) multiples in the period range; and (c) solving the following equation:

(3)

where fi is the frequency of the alias, ft is the true frequency, n is an integer, and fw is a frequency window, using the windows of 1 year (0.0027 day−1), one day (1 day−1), and a sidereal day (1.0027 day−1), as these are the most common aliases for Earth-based telescopes (VanderPlas 2018).

2.6. SED analysis

The target samples are generally considered to represent MYSOs based on signposts of high mass star formation and survey shallowness. To better understand the nature of the sources studied here in detail for variability, we have analysed their 1.2 μm–870 μm spectral energy distributions (SEDs).

The Python version of the SED fitting tool (Robitaille et al. 2007) was used to model SEDs of the target sources. The photometric bands, filter, and apertures used to construct the SEDs can be found on Table 1. The photometric data used to construct the SEDs was obtained from querying the public online archives of 2MASS, Spitzer, ATLASGAL, and Herschel (Huchra et al. 2012; Carey et al. 2009; Schuller et al. 2009; Pilbratt et al. 2010). Our SED fitting follows the method detailed in Grave & Kumar (2009). An uniform photometric error of 10% was assumed. Longer wavelength data were usually set as upper limits, because their large beams include emission from multiple sources, sometimes small clusters, even those that are well resolved at 24 μm and below. Data at wavelengths shorter than 24 μm are set as data points. However, for the EGO sample, the 4.5 μm IRAC band data was set as upper limit by default, as their main characteristic is to have excess emission in that band.

Table 1.

Filters and apertures used for building the SEDs.

We used a range of extinction of Av = 0 − 50 mag for all targets. Distances are available (Urquhart et al. 2018) for 105 targets, non-EGO and EGOs, they were used with an uncertainty of ±1 kpc while fitting. For the remaining 102 targets we allowed a full plausible range of d = 1 − 13 kpc.

For each target all the models which have a were used, and the parameters of the source were computed by performing a weighted mean, weighted by the inverse χ2 as described in Grave & Kumar (2009). The observational data used to construct the SEDs is listed in the Table A.1.

3. Results

The LCs of 601 (448 non-EGO + 153 EGO) were visually examined and compared with the source IQR, while considering the deliberation made in Sect. 2.4. We consistently find that an IQR > 0.05 is associated with visually > 20% of the data-points in the light-curve that are above the 1σ error of the field, as shown by the error bars for each source.

This selection criteria resulted in 51 (of the 448) non-EGO and 139 (of the 153) EGO targets to be classified as variable sources. They are listed in Table 2, along with the LC classification. In Fig. 4 we display some of the clear LCs of both periodic and aperiodic nature. For each source (see Fig. 5), the LC, periodogram, phase-folded LC, and a three colour composite image of the target is made available.

Table 2.

Source co-ordinates, photometric data and variability.

thumbnail Fig. 4.

Some of the clearer LCs, periodic (left column) and aperiodic (right column). Each figure shows the LC of the source, error bars represent MAD(ΔSimjd). The vertical axis represents the variability from the median normalized by max(|Ki−median(K)|).

Open with DEXTER

thumbnail Fig. 5.

Prototypical LPV-yso sources: Top panel: the LC of the source, error bars represent MAD(ΔSimjd). Left middle panel: the corresponding periodogram in logarithmic scale (also plotted are the 99%, 95%, and 90% false probability levels, respectively: the green dot-dashed line, the cyan full line, and the red dashed line). Bottom left panel: the phase-folded light curve of the source using the best period fitted (also shows the corresponding value in days). Bottom right plot: the RGB image of the source using the Spitzer IRAC 3.6 μm, IRAC 4.0 μm, and the 24 μm MIPS band as blue, green and red, respectively. The VVV source is indicated by the blue circle and the green cross represents the MIPS co-ordinates. The contours of the RGB are in the interval of [Peak-5σ, Peak] from the ATLASGAL observation at 850 μm.

Open with DEXTER

3.1. Light curve classification

Light curves can be classified based on their behaviour and, often, such classification represents a close connection with certain physical processes. A classification scheme similar to the one used in Contreras Peña et al. (2017a) was followed here, and LCs were divided into: (a) long period variables (LPV-yso); (b) short timescale variables (STV); (c) dippers and faders; (d) eruptive. In defining periodic variables, Contreras Peña et al. (2017a) only included periods of the highest power, while we include all significant periods. Four LCs for which we have only a short time coverage were considered to be unclassified.

Long period variables (LPV-yso) are defined in Contreras Peña et al. (2017a) as sources with periodic photometric variability and periods larger than P > 100 days. LPV-ysos have periods larger than the stellar rotation or inner disc orbits of young stellar objects, which are typically P < 15 days. Figure 5 shows the example of two LPV-ysos, source G309.91+0.32 and G343.50-0.47, which have periods of ∼545 days, and ∼1156 days. The RGB image of source G309.91+0.32 reveals distinct extended green emission, a signpost of the presence of an outflow, its periodogram shows a prominent signal well above the 99% FAP level. The source has a median brightness of Ks = 13.65 mag in the VVV and the amplitude between the brightest and dimmest point of its LC is ∼0.81 mag. The other prototypical LPV-yso selected, G343.50-0.47, and is part of a complex of three MIPS bright sources. It is a source with Ks = 15.38 mag, the amplitude of its variability is close to ∼0.86 mag and, the periodogram of the source shows a distinct peak well above the 99% FAP level. There are no aliases in the periodograms of either of these sources.

Short timescale variables are objects with short timescales of periodic variability (P < 100 days), or without an apparent periodicity. Periods larger than the stellar rotation or inner disc orbits, 15 < P < 100 days can be explained by phenomena such as obscuration from circumbinary disc or by variable accretion (Contreras Peña et al. 2017a; Bouvier et al. 2003). Sources MG352.2452-00.0636 and G351.78-0.54 are typical examples of STVs, shown in Fig. 6, with typical periods of approximately 29 and 18 days, respectively. Both sources match well (r < 2″) with the 870 μm emission peak, additionally MG352.2452-00.0636 coincides with an IRDC filament, and G351.78-0.54 is close (r < 2″) to the VLA1a source studied by Zapata et al. (2008) as part of a compact cluster of MYSOs. G351.78-0.54 also coincides with an highly variable maser studied by Goedhart et al. (2014).

thumbnail Fig. 6.

Prototypical STV sources: Top panel: shows the LC of the source, error bars represent MAD(ΔSimjd). Left middle panel: the corresponding periodogram in logarithmic scale (also plotted are the 99%, 95%, and 90% false probability levels, respectively: the green dot-dashed line, the cyan full line, and the red dashed line). Bottom left panel: the phase-folded light curve of the source using the best period fitted (also shows the corresponding value in days). Bottom right plot: the RGB image of the source using the Spitzer IRAC 3.6 μm, IRAC 4.0 μm, and the 24 μm MIPS band as blue, green and red, respectively. The VVV source is indicated by the blue circle and the green cross represents the MIPS co-ordinates. The contours of the RGB are in the interval of [Peak-5 * σ, Peak] from the ATLASGAL observation at 850μm.

Open with DEXTER

Faders and dippers are two classes of LCs with aperiodic photometric variability. Dippers are characterized by long-lasting (lasting months to years) dimming events followed by a return to normal brightness, while the same terminology is found in works of the YSOVAR team (Morales-Calderón et al. 2011) they use it to classify phenomena on shorter timescales (hours to days). Faders show light-curves that slowly decline over time, or a period of continuous brightness followed by a sudden decrease sustained over a year. Dippers are often associated with increased extinction from surrounding material. Faders can be caused by either a return to a quiescent accreting phase or a long lasting increase in extinction. It should be noted that both dippers and faders share common LC morphologies and can be easily mistaken for one another. A snapshot of a dipper event not returning to normal brightness can likely be mistaken as a fader. Source MG328.0494-00.0487 (Fig. 7) is the prototypical example of a fader event. It matches the peak emission of the ATLASGAL observations, and is an extended object in the 8 μm Spitzer band, it has a close-by companion. There is no clear peak in its periodogram and the LC shows some periodic variability until around MJD 56500, at which point there is a drop in brightness of close to ΔK ∼ 1.4 mag. The morphology of Dipper events is typified by source MG354.4384+00.4185 which is plotted in Fig. 8. There is no clear peak in its periodogram, and its light curve could be considered as non-variable, except for an abrupt drop in brightness of ΔK ∼ 0.8 until the target recovers more than half its brightness about 750 days later.

thumbnail Fig. 7.

Typical Fader event. Colours and symbols are the same as in Fig. 5.

Open with DEXTER

thumbnail Fig. 8.

LC of a dipper event. Colours and symbols are the same as in Fig. 5.

Open with DEXTER

Eruptive LCs are also aperiodic, but they are characterized by outbursts and increases in brightness, typically over periods of months or years but, in some cases, lasting a few weeks. Objects with increases in their luminosity, likely from ongoing accretion events, will produce such LCs. FUors or EXors are classic examples showing eruptive morphologies. Similarly to Medina et al. (2018) we employ a subdivision of the eruptive class, “low amplitude eruptives” for sources with ΔK < 1.0 mag. This distinction is made to emphasize that such variations are much less extreme than those in FUors and EXors in the optical wavelengths, and more similar to common short term variability in their amplitude. Nevertheless, for certain disk geometries and high extinction it is possible for a FUor or EXor-like eruption to appear as low-amplitude variability in the NIR. A low-amplitude eruptive LC can either correspond to a low-amplitude variable source or to a high-amplitude source with a geometry + extinction combination such that it appears as low-amplitude in the Ks band. Overall, there are 26 low-amplitude eruptives and 15 normal eruptives in our samples. While this identification is presented in Table 2, for the analysis they were not considered as separate classes.

One source that features an ongoing eruptive event is MG303.9304-00.6879, plotted in Fig. 2, showing multiple stages of increased brightness over the entire time-line, with two large amplitude brightness changes over years. Figure 9 shows an example of low-amplitude eruptive LC. This source, G335.59-0.29, displays C-shaped green emission, characteristic of jet emission. The main feature of this LC is its sustained increase in brightness over time, with a total amplitude of ΔK ∼ 0.68 mag. Overall, the variable sources can be split into the periodic category, composed of LPV-yso and STVs, or the aperiodic category, including faders, dippers, and eruptive sources. Their detailed distributions represented by the different MYSO samples can be found in Table 3. The corresponding plots for all sources can be found in Fig. A.1. Analysis of periodogram aliases (see Sect. 2.4) indicates that 1 and 15 members of the non-EGO + EGO samples, respectively, could be classified differently (see also Sect. 4).

thumbnail Fig. 9.

LC of the example of a low-amplitude eruptive event. Colours and symbols are the same as in Fig. 5.

Open with DEXTER

Table 3.

Observed parameters of LC classes, for both EGO and non-EGO samples.

The sources MG300.3241-00.1985, MG322.4833+00.6447, MG342.3189+00.5876 have also been studied as highly variable objects (Contreras Peña et al. 2017a; Kumar et al. 2016). Of these, MG300.3241-00.1985 was studied spectroscopically by Contreras Peña et al. (2017b) and classified as an eruptive MNor, an object with a mixture of characteristics from FUors and EXors. We note that other ΔK > 1 mag sources listed here were not found in Contreras Peña et al. (2017a) because they were not highly variable in the 2010–2012 period.

3.2. Variable source SEDs

The goal of SED fitting was to test if the variable targets indeed represent MYSOs. The SEDs of the variable sources were fitted by YSO models (Robitaille et al. 2006) (see Sect. 2.6) allowing us to constrain the properties of these objects. The results of this fitting procedure can be found in Table A.2. This table contains the full sample of variables with 190 entries. However, as mentioned in Sect. 2.6 only 105 targets have known distances, where the SEDs can be reasonably constrained. We note that in those cases where distances are not available, fitting with the full range of 1–13 kpc has resulted in some model fits that outputs sub-stellar masses. This result is likely to be a consequence of unknown distance rather than the true nature of the source because the indicators of high-mass star formation used in the original selection are more reliable. In Fig. 10, the data and model fits can be visualized for the example targets with different LC classes mentioned in the previous section. Figure A.2 shows the SED fits for all sources. The masses of these example targets range from 1.84 to 10.30 M, with luminosities between 57 and 6918 L, representing evolutionary ages between 104 and 106 yr. Table 4 summarizes the SED results by listing various properties of the sources grouped in mass ranges roughly separating the low, intermediate, and high-mass sources. It can be seen that about ∼35% of the targets are modelled in the 4–8 M range and only 6% representing ≥8 M objects. A large fraction (∼60%) are fitted with YSO models representing sources with M < 4M.

thumbnail Fig. 10.

Grid of SEDs for our prototypical sources. The dark line corresponds to the best fit model. The grey lines correspond to other models.

Open with DEXTER

Table 4.

SED results by mass bin.

The 4–8 M sources display evn ∼ 10−4 M yr−1, Disk ∼ 10−6 M yr−1 and a few hundred solar luminosities. The number of EGO and non-EGO sources fitted as low, intermediate, and high-mass stars are 87, 45, 10 and 25, 21, 1 respectively. It is worth noting that all but one of the sources fitted by models ≥8 M are EGO objects. These SEDs are well-fitted by MYSO models similar to those represented in Grave & Kumar (2009). Four of the 11 objects (≥8 M) are included in the 6.7 GHz class II methanol maser surveys and they show emission. These four also show class I methanol maser emission.

4. Discussion

The results show that 139 of 156 (91%) EGO sample presents variability in contrast to 51 of the 433 (12%) non-EGO targets, implying that variability is strongly correlated with the outflow activity in MYSOs. Table 3 summarizes the variability statistics. More than half (64%) of the variable EGOs are classified as periodic contrasting more than half (59%) of the non-EGO sample that are classified as aperiodic. Table 5 allows us to discern the differences between EGO and non-EGO samples, and sources classified as periodic or otherwise. It can be seen from Fig. 11 and Table 5, that the amplitude range of variation in non-EGOs is roughly twice as much as that of EGOs. Of the modelled parameters, the circumstellar extinction (AVcircum) for non-EGO targets clearly stand out as twice the median value for EGOs. Also, it appears that non-EGO variable sources may simply be more luminous objects located in slightly farther away targets. Together, the ΔKs, Av and L comparison indicate that the non-EGO variable sources are relatively more embedded objects when compared to EGOs.

Table 5.

Summary of the median fit parameters, for both EGO and non-EGO samples divided by periodicity.

thumbnail Fig. 11.

Histogram of ΔK divided by sample and periodicity. EGO and non-EGO sources are shown, respectively, on the right and left plots.

Open with DEXTER

The results of the search for aliases among the ten frequencies with greater power, found aliases for the highest peak of the periodogram in five non-EGO targets (∼9%), and 22 (∼15%) EGO targets. Of these, only 1 (∼2%) of the non-EGO targets would change their classification from LPV-yso to STV, while, for the EGO sample 15 (∼10%) of the targets could change from LPV-YSO to STV or vice versa. Therefore, these aliases would not change any periodic to aperiodic source, as period length is not the only condition defining a periodic source (LC morphology is also one of the main factors).

It can be seen from Fig. 12 that the envelope accretion rate (see also Table 5) for non-EGO sources is an order of magnitude smaller than that for the EGO sources. The same effect can be noticed for aperiodic sources (bottom panel). We note that the non-EGOs are dominated by aperiodic sources, that should be indicative of the differences observed in the top and bottom panels of Fig. 12. Aperiodic LCs, represented by eruptive, dippers and fader classes, are thought to trace objects with low level of quiescent accretion, that will undergo short periods of intense accretion. The lower level of accretion found in these objects can therefore be explained by this behaviour.

thumbnail Fig. 12.

Mass versus envelope-accretion rate for the fitted SEDs of EGO and non-EGO sources, in logarithmic scale. EGO, non-EGO, periodic, and aperiodic, are plotted at the top left, top right, bottom left, and bottom right, respectively.

Open with DEXTER

We compared the SED fitted model properties with the amplitude of variation and did not find any correlations to understand the variability as a function of mass, accretion rate, luminosity or temperature. Figure 13 shows an HR diagram, by plotting the luminosity versus temperature of all variable sources derived from the SED fitting. The zero-age-main-sequence (ZAMS) curve (Siess et al. 2000) is shown by a solid curve. The seven dashed curves display the pre-main-sequence tracks (also from Siess et al. 2000, for solar metallicity) for objects of 1–7 M in steps of 1 M. EGOs are concentrated closer to the putative birth-line position of the massive stars and are also largely lower mass objects (< 4 M). The precursor to a high mass star is considered to be a lower mass object which continues to accrete material for more than half of its life until it contracts on the main sequence. In view of that conjecture, it is not surprising that a majority of the EGO driving sources are modelled by young low to intermediate mass stars. Furthermore, the HR diagram, with the associated PMS tracks, validates the fitted masses.

thumbnail Fig. 13.

HR diagram for our sources. Symbol size corresponds to M < 4, 4 ≤ M < 6, 6 ≤ M < 8, M ≥ 8M, from smaller to larger, respectively. The dashed lines(from bottom to top) are the PMS tracks for 1, 2 , 3, 4, 5, 6, and 7 M, the filled line is the ZAMS. Blue and red symbols are, respectively, EGOs and non-EGOs.

Open with DEXTER

Most parameters resulting from SED fitting are model dependent, with known correlations within the model grid between the age, mass and accretion rates (Robitaille et al. 2006). However, the observed data is scaled to match the luminosity and temperature of the selected models from the grid, therefore, these are relatively more reliable fitted parameters. Unlike the lower mass stars, the luminosity and temperature differences prominently distinguish the sparsely populated massive stellar models in the grid. These two parameters are used for comparison in this analysis, to ensure that the inferences made are reasonably free of biases in the grid models.

There is an apparent concentration of non-EGO objects, probably with slightly more higher mass objects on the ZAMS. It was previously noted that the non-EGO targets are significantly more embedded objects displaying larger ΔKs compared to EGOs. The objects located closer to the ZAMS may therefore be candidate sources to test the hypothesis of bloated and pulsating young massive stars. Hosokawa et al. (2010) argue that high-mass stars are bloated objects. Such objects are also thought to be pulsationally unstable, or at least go through a period of significant pulsations as they settle down on the ZAMS (Inayoshi et al. 2013). Contreras Peña et al. (2017a) indicates that eruptive variable behaviour is more common or recurs more frequently at earlier stages of stellar PMS. The analysis in this work show that ∼70% eruptive variables are concentrated on the birthline and the ZAMS in nearly half proportion. Protostellar envelopes are a prominent feature of objects located on the birthline, therefore suggesting that most eruptive MYSOs are indeed the result of envelope accretion. High-mass protostellar objects ingesting a burst of accreting matter enter a “bloated phase” before re-adjustment and contraction (Hosokawa et al. 2010). This could be the case for those eruptive sources located on the ZAMS.

In Table 6, all the variable sources (32 targets) with known 6.7 GHz class II methanol maser detection are listed, along with, the simultaneous detection of class I methanol maser. Those sources with only class I methanol maser detections are not listed. The detection of class II methanol maser is considered as a strong sign-post of high-mass star formation, especially massive outflow activity (de Villiers et al. 2015). Of the 32 sources, only two are non-EGOs, therefore, reinstating the association of class II methanol masers with MYSO outflow activity. Goedhart et al. (2014) have studied variability of methanol masers, and two of the infrared variable sources presented here, G351.78-0.54 and G298.26+0.74, were analysed in that study and G351.78-0.54 is considered as an highly variable maser, while G298.26+0.74 does not present maser variability above instrumental noise.

Table 6.

EGO and non-EGO MYSO candidates with nearby methanol masers.

Table 7.

Best fit parameters for the SEDs of example targets.

Our selection criteria for non-EGO sources, 24 μm MIPS sources matching ATLASGAL CSC objects (r < 5″) might lead us to miss some of the most important sources in the clumps. Since the most luminous source inside each clump can be offset by more than r < 5″ we have missed many of the MYSOs in these regions. The criteria used ensures that the targets are good MYSO candidates but, the most luminous FIR sources and their counterparts will be examined in a future work.

5. Summary

This study has investigated the nature of near-infrared variability in MYSOs, focusing on the driving sources of EGOs and luminous 24 μm point sources coinciding within 5″ of the massive star forming clumps mapped at 870 μm by ATLASGAL. The search led us to examine the Ks-band light-curves of 718 point sources.

  • 190 sources (139 EGOs and 51 non-EGOs) were found to be variable with an IQR > 0.05 and ΔKs > 0.15. 111 and 79 of these objects are classified as periodic + aperiodic, respectively.

  • The 2 μm–870 μm spectral energy distribution of the variable point sources were assembled and fitted with YSO models. 47 and 6 sources were modelled as ≥4 M and ≥8 M, respectively.

  • On an HR diagram, most lower mass EGO sources concentrate along a putative birth-line.

  • A high rate of detectable variability in EGO targets (139 out of 153 searched) implies that near-infrared variability in MYSOs is closely linked to the accretion phenomenon and outflow activity.

Further to the discovery of a dozen high-amplitude variable MYSOs (Kumar et al. 2016), this is the first large scale systematic study of near-infrared variability in MYSOs. The variable sources identified in this work are excellent targets with which to undertake follow-up studies to understand the circumstellar environment of MYSOs in detail.

Acknowledgments

G.D.C.T. is supported by an FCT/Portugal PhD grant PD/BD/113478/2015. MSNK acknowledges the support from Fundação para a Ciência e Tecnologia (FCT) through Investigador FCT contracts IF/00956/2015/CP1273/CT0002, and the H2020 Marie-Curie Intra-European Fellowship project GESTATE (661249). Support for JB is provided by the Ministry of Economy, Development, and Tourism’s Millennium Science Initiative through grant IC120009, awarded to The Millennium Institute of Astrophysics, MAS. A.C.G. has received funding from the European Research Council (ERC) under the European Unions Horizon 2020 research and innovation programme (grant agreement No. 743029). PWL acknowledges the support of consolidated grants (ST/R000905/1 and ST/M001008/1) funded by the UK Science and Technology Facilities Research Council. CCP acknowledges support from the Leverhulme Trust. JFG is supported by Fundação para a Ciência e a Tecnologia (FCT) through national funds (UID/FIS/04434/2013) and by FEDER through COMPETE2020 (POCI-01-0145- FEDER-007672).

References

Appendix A: Additional figures

thumbnail Fig. A.1.

LC of the source, with error bars representing MAD(ΔSimjd), periodograms (also plotted are the 99%, 95%, and 90% false probability levels, respectively: the green dot-dashed line, the cyan full line, and the red dashed line), the phase-folded LC using the best period fitted, the RGB image of the source using the Spitzer IRAC 3.6 μm, IRAC 4.0 μm, and the 24 μm MIPS band as blue, green and red, respectively. The VVV source is marked by the blue circle and the green cross represents the MIPS co-ordinates. The contours of the RGB are in the interval of [Peak-5σ, Peak] from the ATLASGAL observation at 850 μm.

Open with DEXTER

thumbnail Fig. A.2.

The fitted SEDs with the best fit model in the black line, the grey lines are other models, data points, upper, and lower limits are, respectively, circles, inverted triangles, and triangles.

Open with DEXTER

All Tables

Table 1.

Filters and apertures used for building the SEDs.

Table 2.

Source co-ordinates, photometric data and variability.

Table 3.

Observed parameters of LC classes, for both EGO and non-EGO samples.

Table 4.

SED results by mass bin.

Table 5.

Summary of the median fit parameters, for both EGO and non-EGO samples divided by periodicity.

Table 6.

EGO and non-EGO MYSO candidates with nearby methanol masers.

Table 7.

Best fit parameters for the SEDs of example targets.

All Figures

thumbnail Fig. 2.

LC of an eruptive event. Top panel: the LC of the source, the error bars represent MAD(ΔSimjd ). Bottom plot: the RGB image of the source using the Spitzer IRAC 3.6 μm, IRAC 4.0 μm, and the 24 μm MIPS band as blue, green and red, respectively. The VVV source is indicated by the blue circle and the green cross represents the MIPS co-ordinates. The contours of the RGB are in the interval of [Peak-5σ, Peak] from the ATLASGAL observation at 850 μm.

Open with DEXTER
In the text
thumbnail Fig. 3.

Systematic errors as a function of the magnitude of each target. The red points represents the case where we only consider the stellar sources with ±1 mag around our targets. The blue points represent the case in which we consider all stellar sources in the vicinity of our targets.

Open with DEXTER
In the text
thumbnail Fig. 4.

Some of the clearer LCs, periodic (left column) and aperiodic (right column). Each figure shows the LC of the source, error bars represent MAD(ΔSimjd). The vertical axis represents the variability from the median normalized by max(|Ki−median(K)|).

Open with DEXTER
In the text
thumbnail Fig. 5.

Prototypical LPV-yso sources: Top panel: the LC of the source, error bars represent MAD(ΔSimjd). Left middle panel: the corresponding periodogram in logarithmic scale (also plotted are the 99%, 95%, and 90% false probability levels, respectively: the green dot-dashed line, the cyan full line, and the red dashed line). Bottom left panel: the phase-folded light curve of the source using the best period fitted (also shows the corresponding value in days). Bottom right plot: the RGB image of the source using the Spitzer IRAC 3.6 μm, IRAC 4.0 μm, and the 24 μm MIPS band as blue, green and red, respectively. The VVV source is indicated by the blue circle and the green cross represents the MIPS co-ordinates. The contours of the RGB are in the interval of [Peak-5σ, Peak] from the ATLASGAL observation at 850 μm.

Open with DEXTER
In the text
thumbnail Fig. 6.

Prototypical STV sources: Top panel: shows the LC of the source, error bars represent MAD(ΔSimjd). Left middle panel: the corresponding periodogram in logarithmic scale (also plotted are the 99%, 95%, and 90% false probability levels, respectively: the green dot-dashed line, the cyan full line, and the red dashed line). Bottom left panel: the phase-folded light curve of the source using the best period fitted (also shows the corresponding value in days). Bottom right plot: the RGB image of the source using the Spitzer IRAC 3.6 μm, IRAC 4.0 μm, and the 24 μm MIPS band as blue, green and red, respectively. The VVV source is indicated by the blue circle and the green cross represents the MIPS co-ordinates. The contours of the RGB are in the interval of [Peak-5 * σ, Peak] from the ATLASGAL observation at 850μm.

Open with DEXTER
In the text
thumbnail Fig. 7.

Typical Fader event. Colours and symbols are the same as in Fig. 5.

Open with DEXTER
In the text
thumbnail Fig. 8.

LC of a dipper event. Colours and symbols are the same as in Fig. 5.

Open with DEXTER
In the text
thumbnail Fig. 9.

LC of the example of a low-amplitude eruptive event. Colours and symbols are the same as in Fig. 5.

Open with DEXTER
In the text
thumbnail Fig. 10.

Grid of SEDs for our prototypical sources. The dark line corresponds to the best fit model. The grey lines correspond to other models.

Open with DEXTER
In the text
thumbnail Fig. 11.

Histogram of ΔK divided by sample and periodicity. EGO and non-EGO sources are shown, respectively, on the right and left plots.

Open with DEXTER
In the text
thumbnail Fig. 12.

Mass versus envelope-accretion rate for the fitted SEDs of EGO and non-EGO sources, in logarithmic scale. EGO, non-EGO, periodic, and aperiodic, are plotted at the top left, top right, bottom left, and bottom right, respectively.

Open with DEXTER
In the text
thumbnail Fig. 13.

HR diagram for our sources. Symbol size corresponds to M < 4, 4 ≤ M < 6, 6 ≤ M < 8, M ≥ 8M, from smaller to larger, respectively. The dashed lines(from bottom to top) are the PMS tracks for 1, 2 , 3, 4, 5, 6, and 7 M, the filled line is the ZAMS. Blue and red symbols are, respectively, EGOs and non-EGOs.

Open with DEXTER
In the text
thumbnail Fig. A.1.

LC of the source, with error bars representing MAD(ΔSimjd), periodograms (also plotted are the 99%, 95%, and 90% false probability levels, respectively: the green dot-dashed line, the cyan full line, and the red dashed line), the phase-folded LC using the best period fitted, the RGB image of the source using the Spitzer IRAC 3.6 μm, IRAC 4.0 μm, and the 24 μm MIPS band as blue, green and red, respectively. The VVV source is marked by the blue circle and the green cross represents the MIPS co-ordinates. The contours of the RGB are in the interval of [Peak-5σ, Peak] from the ATLASGAL observation at 850 μm.

Open with DEXTER
In the text
thumbnail Fig. A.2.

The fitted SEDs with the best fit model in the black line, the grey lines are other models, data points, upper, and lower limits are, respectively, circles, inverted triangles, and triangles.

Open with DEXTER
In the text

Current usage metrics show cumulative count of Article Views (full-text article views including HTML views, PDF and ePub downloads, according to the available data) and Abstracts Views on Vision4Press platform.

Data correspond to usage on the plateform after 2015. The current usage metrics is available 48-96 hours after online publication and is updated daily on week days.

Initial download of the metrics may take a while.