WISE Discovery of Mid-infrared Variability in Massive Young Stellar Objects

Mizuho Uchiyama; Kohei Ichikawa

doi:10.3847/1538-4357/ab372e

1. Introduction

Massive stars have a significant role in star formation activities and metal enrichment in galaxies and, hence, the evolution of galaxies in the universe (e.g., Zinnecker & Yorke 2007; Tan et al. 2014). The powerful outflows of massive stars during the star-forming phase and subsequent expansion of H ii regions trigger next-generation star formation. At the end of their evolution, these stars supply heavy elements to circumstellar regions through stellar winds and supernova explosions, resulting in enriched chemical environments (e.g., Nozawa et al. 2007). However, the detailed mass growth processes have yet to be resolved due to observational difficulties. The birthplace of such massive stars, called massive young stellar objects (MYSOs), are deeply obscured by dense gas and dust, preventing observation in the optical and even sometimes near-infrared (NIR) regimes (Grave & Kumar 2009). In addition, because they are relatively rare and located at far distances (>1 kpc), it is difficult to spatially resolve the physical and geometrical structures within <10 au (Tan et al. 2014), where mass growth phenomena occur. Notably, these features are smaller than the currently achievable highest spatial resolution (20 mas at 230 GHz in the Atacama Large Millimeter/submillimeter Array extended configuration). Very long baseline interferometry observations, such as Motogi et al. (2016) in the H₂O 22.2 GHz band, have to date been the only means to resolve such compact structures, in which mostly non-thermal emissions, such as masers, have been detected. Physical parameters and structures of MYSOs, such as disk gas/dust density and temperature structure, stellar effective temperature and radius, and mass accretion processes, are still largely unknown. Thus, alternative observational methods are required to determine the properties of the spatially compact inner structures of MYSOs.

Variability studies may provide a way to address such difficulties. Observations in both the optical and NIR bands of low-mass young stellar objects (YSOs) have been conducted over the last 50 yr via the small extinction of dust with A_V < 10 (Ménard & Bertout 1999). Such optical and NIR variability studies involving the emission lines and continuum of lower-mass YSOs have proven to be powerful tools for deciphering the physics of star formation and pre-main-sequence stellar evolution, such as the existence of rotating cool/hot spots (Bouvier & Bertout 1989), infalling/rotating dusty objects (Cody et al. 2014), and accretion rate variation (Audard et al. 2014). Recent intensive, large, and deep mid-infrared (MIR; λ > 3 μm) variability surveys of low-mass YSOs have enabled statistical studies, such as the Young Stellar Object VARiability (YSOVAR) project with Spitzer (Morales-Calderón et al. 2011). This statistical study revealed that flux variability in low-mass YSOs is ubiquitous in the optical, NIR, and MIR regimes. They also categorized the causes of flux variability, such as stellar activity (e.g., cool spots), the density structure of disk, and mass accretion processes (e.g., hot spots, dust infalling, and episodic accretion), based on time and color variability characteristics (Günther et al. 2014).

Despite their importance, MYSOs have not been investigated extensively, especially in terms of IR thermal emission. However, recently, high-amplitude (ΔK_s > 1 mag) year-scale variability in the NIR K_s was reported for the first time, in 13 MYSOs using the Vista Variables in Via Lactea (VVV) survey data (Kumar et al. 2016). Following that study, lower-amplitude variables, approaching ΔK_s > 0.15 mag, have been detected in 190 out of 718 MYSO candidates in VVV survey data (Teixeira et al. 2018). These studies suggest that the variability is also common for MYSOs. Recent theoretical studies also indicated that variable accretion rate and luminosity are common during mass growth phase of MYSOs (e.g., Meyer et al. 2019) and support the above observational results. However, sample sets are still limited in these works, due to the survey area of the Galactic center, 520 degree², and the difficulty in detection, as most MYSOs are still highly obscured in the NIR (Grave & Kumar 2009). Thus, a flux variability survey that (1) covers a much more extensive area of the sky and (2) includes wavelengths longer than those of the NIR band is essential to the study of the properties of massive-star formation and mass growth processes.

In this paper, we report the first discovery of MIR-variable MYSOs using the the Red MSX Source (RMS) survey MYSO catalog, by cross-matching the MIR counterparts with the multi-epoch ALLWISE archives. We also discuss the physical properties of detected-variable MYSOs using ALLWISE and NEOWISE data and their possible origins, based on the obtained light curves and color–magnitude diagrams (CMDs).

2. Sample

Here we describe the sample in this study and the selection process to find MIR-variable sources. The flow chart of sample selection is compiled in Figure 1.

**Figure 1.** Summary flow chart of the sample selection process. The detail of each process is discussed in Section 2.
Download figure:
Standard image High-resolution image

2.1. RMS Database

Our initial sample was from the RMS MYSO Database (Lumsden et al. 2013), which is a multi-wavelength Galactic-plane MYSO survey based on Midcourse Space Experiment (MSX) IR satellite data archived at 8.28, 12.13, 14.65, and 21.3 μm (Egan et al. 2003). It is the largest statistical catalog of young massive protostars with MIR detection, and the most suitable one for our purpose, i.e., searching for MIR-variable MYSOs in large samples. It covers almost the entire Galactic plane, with the exception of the extremely crowded region of 10° < l < 350°, with a latitude of $| b| \lt 5^\circ$ (Egan et al. 2003). In this catalog, 806 MYSOs are listed by multicolor selection using the MSX and 2MASS data archives with criteria of F₂₁ > 2F₈, F₂₁ > F₁₄, F₁₄ > F₈, F₈ > 5F_K, and F_K > 2F_J (Lumsden et al. 2002, 2013). These criteria effectively remove contaminants, such as evolved stars with dust shells. The RMS survey team also conducted follow-up radio observations; any sources detected at 5 GHz above 10 mJy were considered to be an H ii region or planetary nebulae (Urquhart et al. 2007a, 2009; Lumsden et al. 2013); hence, the contamination from the H ii region/planetary nebulae should be minimal in the data. There are some known MYSOs associated with jets and detection at radio frequencies (e.g., Guzmán et al. 2010; Lumsden et al. 2013). However, as all are well below 10 mJy, the catalog accounts for such MYSOs based on the radio detection criteria specified above.

2.2. Cross-matching of the RMS Catalog with ${WISE}$

We determined the MIR (band W1–W4; 3.4–22 μm) counterparts of the RMS MYSOs through positional matching with ALLWISE (for all four bands; Wright et al. 2010) and NEOWISE (only for W1 and W2; Mainzer et al. 2011, 2014). We used NEOWISE 2019 Data Release version⁴ for the following analysis. Hereafter, we will refer to both IR data sets as WISE.

In this study, we applied a cross-matching radius of 2'', based on the positional accuracy obtained with the 2MASS catalog (see also Ichikawa et al. 2012, 2017). When multiple WISE objects satisfied this criterion, we removed them from our sample due to unreliable point-spread function (PSF) photometry in WISE data reduction pipelines. After this matching, 660 MYSOs showed corresponding WISE objects. We only used sources with flux quality ph_qual=A, with a signal-to-noise ratio larger than 10.0, and applied the single-exposure images of qi_fact=1.0, which have the best image quality. To select data which was not affected from artificial source contamination, South Atlantic Anomaly, or moon light, qual_frame>0.0, saa_sep>5.0, and moon_masked=0 were applied. We also checked for sources of contamination and/or biased flux in the vicinity of an image artifact (e.g., diffraction spikes, scattered-light halos, and/or optical ghosts), using the contamination flag cc_flags. We removed sources containing a contamination flag with any letters, but retained sources flagged as cc_flags=0 that are known to be unaffected by the known artifacts. We then selected the sources with saturated pixels less than 5% within the photometric apertures, flagged as w1/2_sat<0.05 to remove unreliable photometry with heavy saturation. Saturation begins to occur for point sources brighter than ∼8 and 7 mag in W1 and W2 bands, respectively.⁵ We removed sources that failed background sky fitting or had a PSF profile fitting that was not good, flagged as w1/2_sky=NAN or w1/2_rchi2>150. These selection processes finally leave 331 sources, as shown in Figure 1.

2.3. Variability Selection

MIR variability signals were checked for the selected 331 sources having reliable multi-epoch WISE photometries. In this study, the only ALLWISE data were employed for MIR variability selections, due to the stability and accuracy of the photometries. In addition, we utilized only W1- and W2-band data, as multi-epoch data for these bands are readily available. Through this process, objects with clear MIR variability signals were chosen, and hence note that sources with possible variability might have been missed in our study.

WISE has a 90 minute orbit and conducts ≈12 observations of a source over a ≈1 day period. A given location is observed every 6 months. In this study, we refer to this set of ≈12 observations on 1 day as a "1 epoch" observation. The average magnitude and dispersion (σ) of each epoch observation was derived to determine the significant variability in each MYSO (see also Figure 2). As ALLWISE covers observations between 2010 January and 2011 February (Modified Julian Date (MJD)–55,000 = 200–600), each source has observations for two epochs.

**Figure 2.** Infrared light curves of variable sources in W1 and W2 bands. The single exposures with error values are shown in dark-blue (W1) and dark-red (W2), respectively, and the averaged value in each epoch is shown in the light-blue (W1) and pink (W2) circles with error bars. The averaged values are shifted to 100 days after the real values for clarity.(The data used to create this figure are available.)
Download figure:
Standard image High-resolution image

Using the ALLWISE two-epoch observations, we searched the variable MYSOs to determine which ones satisfied the four criteria listed below:

1.
Difference in average magnitudes between the two epochs is larger than 3σ in each epoch observation.
2.
Difference in average magnitudes between epochs is larger than 0.3 mag.
3.
σ in each epoch is less than 0.3 mag.
4.
Number of secured photometric data in each epoch is larger than two points.

Out of 331 sources, 12 fulfilled the above criteria.

In addition to applying the criteria of RMS MYSO classification, we also removed the possible low-mass YSOs. The RMS project provides the integrated bolometric luminosity of each source using archive MIR/far-infrared (FIR) satellite data and the associated kinematic distances that they determined (Urquhart et al. 2007b, 2008; Lumsden et al. 2013). Objects with luminosity less than 500 L_⊙ were also removed, as they were likely lower-mass YSOs with corresponding masses of less than about 5 M_⊙ during the zero-age main sequence (Cox 2000; Siess et al. 2000). Because MYSOs increase their bolometric luminosity almost monotonically during the mass accretion phase (Hosokawa & Omukai 2009), the threshold was set to be slightly below the exact massive-star definition, 8 M_⊙, so as not to miss potential MYSOs. This criterion left the final 5 MIR-variable candidates out of 12 sources.

We checked whether there is a source of equal brightness in the W1/W2 band within about 15'' or a brighter one within 25'' and confirmed that such source does not exist for all variable candidates in the ALLWISE catalog. We also checked all variable candidates by visual inspection using the WISE image server,⁶ to investigate whether the possible artifacts or diffuse components contributed significantly. None of the candidates were rejected in these two criteria. After applying the criteria described above, ultimately, five objects remained as MIR-variable MYSOs.

3. Results and Discussion

Figure 2 shows the light curves for the five MIR-variable sources of ALLWISE, covering observations between 2010 January and 2011 February (MJD−55,000 = 200–600), and NEOWISE, covering observations between 2013 December 13 and 2018 December 13, UTC (MJD−55,000 = 1600–3500). The variability properties of each variable object are summarized in Table 1.

Table 1. Properties of Selected Variable Sources

(1)	(2)	(3)	(4)	(5)	(6)	(7)	(8)	(9)	(10)	(11)	(12)	(13)	(14)
Object Name	R.A.	Decl.	W1⁽¹⁾	W2⁽¹⁾	W3⁽¹⁾	W4⁽¹⁾	W1⁽²⁾	W2⁽²⁾	Δ W1	Δ W2	Variability	Color	L
			(mag)	(mag)	(mag)	(mag)	(mag)	(mag)	(mag)	(mag)			(L_⊙)

G076.1807+00.0619	306.012	37.610	8.57	6.85	3.65	0.67	⋯	6.50	⋯	−0.35	D Pl	BB	5.45 × 10²
G269.5205-01.2510	136.123	−48.823	8.63	6.36	2.96	−0.84	9.29	⋯	0.66	⋯	O	Unclear	1.0 × 10³
G323.7986+00.0173	232.738	−56.250	8.05	⋯	2.34	−0.44	7.70	⋯	−0.35	⋯	I Pl	⋯	4.58 × 10⁴
G335.9960-00.8532	248.795	−48.814	7.94	⋯	3.01	0.63	8.84	⋯	0.91	⋯	P	⋯	9.34 × 10²
G343.1880-00.0803	254.642	−42.832	8.28	⋯	3.84	1.66	9.47	⋯	1.17	⋯	Dip	⋯	7.58 × 10²

Note. Columns: (1) Object name; (2) R.A.; (3) decl. (J2000); (4)–(7) profile-fitting flux densities from W1 to W4 obtained by ALLWISE for the first epoch, while only W1 and W2 are selected data as described in Section 2.2 ; (8)–(9) profile-fitting flux densities from W1 and W2 obtained by ALLWISE for the second epoch using selected data; (10)–(11) infrared (IR) color of the different epochs, defined as ΔW1 (ΔW2) = W1⁽¹⁾ (W2⁽¹⁾) − W1⁽²⁾ (W2⁽²⁾). (12) Variability type flag as discussed in Section 3.1.1. "P" represents "periodic variability," "dip" as "dipper," "I Pl" as "increase with plateau," "D Pl" as "decrease with plateau," and "O" as "other." (13) Color variability type as discussed in Section 3.1.2. "BB" represents "bluer when brighter." (14) Estimated bolometric luminosity with units of solar luminosity, obtained from Lumsden et al. (2013).

Download table as: ASCII Typeset image

3.1. Classification of Variability

To investigate the origins of the detected variability, each MYSO is categorized with respect to its time and color flux variation.

3.1.1. Time Variation

Figure 2 shows that there are distinct types of variabilities. To discuss each type of variability, the time variation of each MYSO is categorized as "periodic," "dipper," "increase/decrease with plateau," or "other."

In this work, the Lomb−Scargle method (Lomb 1976; Scargle 1982) was used to search for periodic variability with the AstroPy Lomb−Scargle periodogram package (Astropy Collaboration et al. 2013). This method is widely used to detect periodic variability in unevenly obtained data, which is seen frequently in astronomical observations (e.g., Morales-Calderón et al. 2011; Goedhart et al. 2014; Sugiyama et al. 2017). Considering that the cadence of the IR observations is roughly every half a year (∼180 days), the minimum detectable periods should be larger than twice the cadence, i.e., 360 days. On the other hand, our available data set spans up to nine years (∼3300 days), indicating that the upper-bound of the detectable period should be smaller than that of half of the obtained data set, i.e., 1650 days. Therefore, periods between 360 and 1650 days were searched in WISE data set in this work. When conducting a Lomb−Scargle periodogram analysis, power spectra versus frequency is derived as the calculation result. To evaluate the peaks that are significant in the power spectra, we used false alarm probability (FAP), which represents the probability of attaining a certain amplitude of power spectra without any periodic component, assuming that the obtained data error consists of only Gaussian noise (VanderPlas 2018). In this work, an FAP of less than 3 × 10⁻³, 3σ detection assuming Gaussian noise distribution, was adopted as the threshold of periodic variability detection.

The criteria for each variability class are given below.

1.
If the FAP of the Lomb−Scargle method periodogram is less than 3 × 10⁻³ in either band, the object is categorized as "periodic."
2.
Or if the flux drop between adjacent epoch magnitudes is larger than 1 mag in either data set, the object is categorized as "dipper."
3.
Or if the binning fluxes partially increase or decrease while the other epochs are overlapped within σ, the object is categorized as "increase/decrease with plateau."
4.
The object which is not categorized as one of the above criteria is labeled as "other."

The classification results are summarized in Table 1. There are two "plateau," one "dipper," one "periodic," and one "other" object.

3.1.2. Color Variation

A color change during the flux variation provides us important information on the possible origins of the variability. When the flux variation is caused by either a change in the dust temperature of emitting regions, possibly due to mass accretion-rate variation, or the level of dust extinction, its color should also vary simultaneously. In contrast, the occultation by binary stars or planets in the line of sight affects the normalization of the flux, but does not produce a color change between W1 and W2.

Figure 3 shows the two CMDs of W1−W2 and W1 for our MIR variability sources, which have selected more than four W1 and W2 data at the same observation time, G076.1807+00.0619 and G269.5205-01.2510. According to the MIR color variability analysis in YSOVAR project (Günther et al. 2014), we used Orthogonal Distance Regression (ODR) method to derive a slope of each CMD. When analyzed data of W1−W2 and W1 is noise limited, then an artificial slope of 45° is obtained (Günther et al. 2014). Therefore, we derive slopes of CMD to investigate whether obtained values could be 45° within the ODR fitting error.

Our fitted slopes of G076.1807+00.0619 and G269.5205-01.2510, are over-plotted in Figure 3 and their equations are y = (3.2 ± 1.3)x + (2.8 ± 2.3) and y = (6.2 ± 5.0)x − (4.8 ± 11), respectively, where x and y are a color of W1−W2 and a magnitude of W1. In Figure 3, an extinction vector at R_V = 5.5, typical value at both interstellar medium and dense molecular clouds in 3–8 μm (Weingartner & Draine 2001; Wang et al. 2014), is also shown. G076.1807+00.0619 object shows a bluer (redder) color when the flux is increased (decreased) and the obtained slope is steeper than the 45° including fitting error. Therefore, the obtained slope shows a "bluer when brighter (BB)" color trend. On the other hand, the slope of G269.5205-01.2510 has a large error in the fitting. Therefore, we could not find a statistically significant color trend for G269.5205-01.2510.

3.2. Fraction of MIR-variable MYSOs and Comparison with Lower-mass YSOs

Our results show that the fraction of variable MYSOs is 5/331, or equal to 1.5%, which seems to be much smaller than that found by lower-mass variables such as 70% at YSOVAR (Morales-Calderón et al. 2011). However, year-scale variables are rare even in lower-mass YSOs, because most of the variability in lower-mass YSOs are related to rotation, and hence its periods are a few days or less (Morales-Calderón et al. 2011; Rebull et al. 2014). The fraction of year-scale variability (Δ [3.6] ≳ 0.05 mag) in YSOVAR reduces into roughly 20%–30% (Rebull et al. 2014). While the value is still lager than our result of 1.5% ( ${\rm{\Delta }}\left([3.4]/[4.6]\right)\gt 0.3$ mag), our detection-threshold requires higher amplitude with >0.3 mag, which might miss possible variabilities with smaller amplitude of $0.05\lt {\rm{\Delta }}\left([3.4]/[4.6]\right)\lt 0.3$ mag as seen in the sample of YSOVAR. Actually, amplitudes of year-scale flux variations are generally weaker than those of day-scale ones in lower-mass YSOs according to the Research Of Traces Of Rotation (ROTOR)-program in V band observations (Grankin et al. 2007) and if this is the case with MYSOs, the current selection method may miss many year-scale variables, although such ambiguous variables are out-of-scope in this paper. Thus, the obtained variability fraction of 1.5% in this study could be considered as a lower-limit.

3.3. Characteristics of Variability in MYSOs

The cross matrix of time variability (see Section 3.1.1) and color variability (see Section 3.1.2) is shown in Table 2. In this table, while about three objects were not able to conduct ODR analysis and one object does not show clear color trend in ODR analysis, one "plateau" object shows "BB" color tendency, suggesting that an origin of those MIR variabilities is from either the extinction and/or temperature change of the MIR emission region of the source, such as a change in accretion rate.

Table 2. Classification Based on the Time and Color Variations of MYSOs in this Study

	BB	Unclear	-	Total
Periodic	0	0	1	1
Dipper	0	0	1	1
Plateau	1	0	1	2
Other	0	1	0	1
Total	1	1	3	5

Note. "BB" represents "bluer when brighter" color variability ; "-" represents unavailable CMD data. Detailed classification criteria are described in Sections 3.1.1 and 3.1.2.

Download table as: ASCII Typeset image

G335.9960-00.8532 showed clear periodic variability in the W1 band, with an FAP of less than 3 × 10⁻³; however, no color variability was available because of the saturation of W2 band data. Detailed stellar parameters and possible origins of periodic variability are discussed in Section 3.4.4.

3.4. Possible Origins of Flux and MIR Color Variation

3.4.1. Plateau Class with BB Color

Long-term variability has been studied in lower-mass YSOs. ROTOR-program shows year-scale variability in the optical UBVR bands for classical T-Tauri stars (CTTs) from their 20 yr monitoring program (Grankin et al. 2007). Approximately 75% of the long-term variability observed in those CTTs shows moderate and almost monotonic flux variation with ΔV ≤ 1.6 mag, and color variations are consistent with the interstellar extinction vector. This suggests that the main cause of long-term variability in CTTs is from a difference in line-of-sight extinction or variation in the mass accretion rate (Grankin et al. 2007). The color variability trends which are consistent with the interstellar extinction vector or somewhat much more colorless ones are also generally found in MIR variability study of lower-mass YSOs (Günther et al. 2014). They also suggested the difference in line-of-sight extinction or variation in the mass accretion rate. Although a 90% quantile of obtained A_V corresponds to ΔV ≤ 10 mag and this value is larger than that in the ROTOR-program, this may be the result of much younger objects included in the work of Günther et al. (2014).

A similar trend is evident in G076.1807+00.0619, which shows "plateau" time-variability MYSO associated with the "BB" color variability; this suggests that these long-term variabilities may be caused by the transition of the line-of-sight extinction, such as the rotation of the accretion disk with a non-axisymmetric angular density distribution or change in mass accretion rate, possibly caused by binary interaction (Araya et al. 2010) due to high close binary rate in massive stars (Chini et al. 2012), or gravitational instability of accretion disk due to high accretion rate (Matsushita et al. 2017; Meyer et al. 2019). Assuming that the accretion disk of those MYSOs is Keplerian rotation with a central MYSO mass of ∼10M_⊙ and bolometric luminosity of ∼5000L_⊙, the dynamical timescale of the observed variation of ∼1000 days in this study corresponds to a radius of R_dyn ∼ 10 au, which is five-fold larger than the dust sublimation radius R_sub of ∼2 au (R_dyn ∼ 5R_sub). Further spectroscopic monitoring of hydrogen recombination lines related to the mass accretion rate (e.g., Alcalá et al. 2014) is critical to resolve the degeneracy of the two possible origins between the extinction and the change in accretion rate.

3.4.2. Dipper Class

One detected "Dipper" class object, G343.1880-00.0803, shows a sudden flux decrease with subsequent flux recovery. The above characteristics well resemble the "dipper" objects in lower-mass YSOs (Bibo & The 1990; Waters & Waelkens 1998; Cody et al. 2014; Günther et al. 2014), exhibiting a flux decrement due to the sharp line-of-sight extinction, while the CMD of G343.1880-00.0803 was unavailable.

3.4.3. Other Class

The object, G269.5205-01.2510, is categorized as "other" in terms of time variability. This would be mainly associated with the insufficient cadence or the duration of observations. While the origin of the WISE MIR color trend in this object is not clear yet, the K−W1 color referred from the RMS survey is reported to be very large, about 4–5 mag.⁷ Again, its origin is not clear at the current stage. Thus, further observations with higher cadence and/or longer duration in the NIR and MIR should reveal a hidden variability trend with a shorter timescale.

3.4.4. Periodic Variable MYSO Candidate G335.9960-00.8532

According to the periodic analysis described in Section 3.1.1, a clear periodic variability is detected in G335.9960-00.8532 (hereafter, G335). G335 has a period of 693.3 days in the W1 band, with an FAP of 5.6 × 10⁻⁴. Assuming its light curve is approximated as a sinusoidal curve, the best-fitting modeled light curve is overlaid in Figure 4. While the periodicity is clearly shown in Figure 4, the CMD and hence the MIR color trend was not available for G335 due to saturation in W2 band.

**Figure 4.** Best-fitted period of G335.9960-00.8532 in the W1 band. Periods of 693.3 days were derived.
Download figure:
Standard image High-resolution image

In further examination of the possible origins of periodic variability, we estimated stellar and accretion disk parameters, using a spectral energy distribution (SED) fitting based on the MYSO models, including various combinations of stellar parameters, such as extinction, total bolometric luminosity, stellar mass, effective temperature, and disk mass (Robitaille et al. 2007).

The SED of G335 consists of the archival data of MSX (RMS project), Spitzer, and Herschel (Hi-GAL program; Wright et al. 2010; Lumsden et al. 2013; Molinari et al. 2016) by cross-matching the coordinates obtained by the RMS catalog with the searching radius of 5''. SED-fitting results are shown in Figure 5, and the derived stellar parameters are summarized in Table 3. Because the distance to the object is obtained from the kinematic near-side distance (Urquhart et al. 2007b), the obtained bolometric luminosity and the other luminosity-related parameters correspond to the lower bound. The associated stellar parameters indicate that G335 contains a central MYSO of 10.5 M_⊙, with an associated small mass disk of ∼4.2 × 10⁻³ M_⊙. This indicates that the evaporation is due to the irradiation of the central MYSO or late-phase MYSO evolution. Although the model fitting shows that the temperature of the central MYSO is well above 10⁵ T, the H ii region cannot be detected based on high-spatial resolution radio observation (Urquhart et al. 2007a). One possible origin is that during the late-phase MYSO evolution, as the central engine begins to blow away the surrounding dust region, the central H ii region remains compact and sufficiently dense, such that it cannot be detected by radio observation.

**Figure 5.** Obtained spectral energy distribution (SED) for G335.9960-00.8532 and the best-fitting result. The black points with the error bars correspond to the obtained infrared (IR) photometries from 2MASS, *MSX*, *Spitzer*, and *Herschel*/PACS. As all 2MASS bands are non-detections, we show them as upper bounds, with the triangles of the *J, H,* and K bands. The best-fitting SED is indicated by a solid black line.
Download figure:
Standard image High-resolution image

Table 3. The Main Parameters for the Best-fitting Models of G335.9960-00.8532

Columns	Parameter	G335
General

(1)	A_V	69.7
(2)	Evol. Stage	Disk
(3)	L_tot/L_⊙	7.6 × 10³

Central Source

(4)	M_⋆/M_⊙	10.6
(5)	R_⋆/R_⊙	4.13
(6)	T_⋆/K	2.65 × 10⁴

Disk Component

(7)	M_disk/M_⊙	4.2 × 10⁻³
(8)	R_disk,in/au	61.4
(9)	R_disk,out/au	6.39 × 10³

Note. Columns: (1) foreground dust extinction. (2) Evolutional stage of the YSOs. "Disk" refers to sources that are surrounded only by a circumstellar disk. (3) Total bolometric luminosity of G335. (4) Central stellar mass. (5) Central stellar radius. (6) Central stellar temperature. (7) Disk mass. (8) Disk inner radius. (9) Disk outer radius.

Download table as: ASCII Typeset image

Some origins may provoke periodic MIR flux variation, as observed in G335. The stellar rotation with cool/hot spots is a well-observed phenomenon in lower-mass YSOs, as periodic day-scale variability (Morales-Calderón et al. 2011). However, because the observed period for G335 is much longer at >100 days, an origin of stellar rotation is unlikely (Rosen et al. 2012). The other interesting idea is a close binary origin. Because main-sequence massive stars frequently have binary stars, a close binary system could be easily realized for G335 (Chini et al. 2012). However, it is difficult for occultation of the binary system to achieve a sinusoidal-like light curve, as observed in Figure 4.

Close interaction of a MYSO binary system could produce the periodic change in mass accretion rate, and thus a resulting sinusoidal-like light curve as suggested in Araya et al. (2010). Since the change in mass accretion rate causes "BB" color variability in the MIR (Günther et al. 2014), multicolor monitoring observations are necessary to evaluate this hypothesis. Similarly, recent theoretical studies indicate that a change in accretion rate caused by either a disk-fragmented companion in disk or gravitational disk instability can occasionally make sinusoidal-like light curves (Matsushita et al. 2017; Meyer et al. 2019), although the current spatial resolution of those simulations is not high enough to resolve an inner 10 au region. Therefore, our results would motivate such simulation studies to explore more finer spatial resolutions by resolving the year-scale flux variation.

Another possible scenario is the periodic/aperiodic change in extinction in the line of sight, such as the non-axisymmetric dust density distribution in the rotating accretion disk. The CMD diagram gives us to judge whether it is also the case for G335, while the CMD is not available because of the saturation of W2 band.

Recent theories have proposed that stellar pulsation occurs with late-phase accreting MYSOs (e.g., Inayoshi et al. 2013); the obtained SED-fitting result for G335 is consistent with the scenario proposed above. In the prediction, a ∼100 day scale periodic variation that follows a period–luminosity relation was assumed; thus, the observed period for G335 should be the luminosity of 2.8 × 10⁵L_⊙, which is one order of magnitude larger than the bolometric luminosity of 7.6 × 10³L_⊙ obtained by the SED fitting. When the far-side distance, 12.1 kpc, was employed, the bolometric luminosity became 1.0 × 10⁵L_⊙, consistent with the period–luminosity relation within a factor of three. However, it would be difficult to observe the MYSO located at 12.1 kpc beyond the Galactic center region, even in the MIR and FIR regions. Because the current model of Inayoshi et al. (2013) assumes a constant accretion rate, it may affect the predicted relation. From a more dynamic perspective, a time-dependent accretion model for stellar pulsation may be necessary to fill in the gap between theory and observation for G335.

Follow-up observations in the NIR and MIR regimes of G335 are planned for upcoming ground telescope observations through The University of Tokyo Atacama Observatory (TAO)/Simultaneous-color Wide-field Infrared Multi-object Spectrograph and Mid-Infrared Multi-mode Imager for gaZing at the UnKnown Universe (Kamizuka et al. 2016; Motohara et al. 2016). To evaluate variations in the mass accretion rate and dust extinction during flux variability, the NIR spectroscopic monitoring of hydrogen recombination lines and the MIR narrow-band imaging of dust silicate features are critical. Monitoring observations in the 10 μm band is currently difficult, due to a lack of access to the MIR instrument. The new ground telescope, TAO, will enable the acquisition of >1 yr scale monitoring data.

4. Summary and Conclusions

We conducted a systematic search for long-term (>1 yr scale) MIR (λ > 3 μm) variability of MYSOs utilizing all-sky ALLWISE and NEOWISE data. This provided a periodic MIR light curve at the bands W1 (3.4 μm) and W2 (4.6 μm), with a cadence of every 6 months during the 9 yr long coverage period. We cross-matched between the WISE catalogs and 806 MYSOs obtained by the RMS survey, the largest MYSO catalog in the MIR band with wide wavelength coverage, including those of the radio band. After the cross-matching process and the variability selection, we found five MIR-variable sources.

We categorized the time variability classes based on the obtained 9 yr long MIR light curves into periodic, dipper, plateau, and others. G335.9960-00.8532 (G335) shows 693.3 day periodic variability in the W1 band, two sources show a plateau-like light curve, and one shows "eruption"-like, and one is categorized as "other."

For the "plateau" class, one out of the two objects showed "BB" MIR color trends, suggesting that the MIR variability originates from a change in the accretion rate and/or dust extinction, which is also the case with many lower-mass MIR-variable YSOs. The "dipper" object G343.1880-00.0803 showed a drastic flux change, similar to "dipper" objects in lower-mass YSOs, whereas its CMD is not available and possible origin is not currently clear.

For G335, the only "periodic" class, we applied SED model fitting to estimate the stellar and accretion disk parameters. Our results indicated that G335 is in the phase of a relatively evolved MYSO; thus, we may be witnessing the very early stages of a hyper- or ultra-compact H ii region.

The remaining objects, G269.5205-01.2510, categorized as "other" class and showing no clear color trend because of the large scatter, are good candidates for future follow-up observations with various cadence and multi-band observations, including spectroscopy to determine the time and color MIR variabilities to resolve their origins.

We thank Kohei Inayoshi and Kei Tanaka for discussions and feedback, which improved our paper. We appreciate the support of Sizuka Nakajima for collection of important data. This study benefited from financial support from the Japan Society for the Promotion of Science (JSPS) KAKENHI program (18K13584; KI), and the Japan Science and Technology Agency (JST) grant "Building of Consortia for the Development of Human Resources in Science and Technology" (KI). This publication makes use of data products from the Near-Earth Object Wide-field Infrared Survey Explorer (NEOWISE), which is a project of the Jet Propulsion Laboratory (JPL)/California Institute of Technology. NEOWISE is funded by the National Aeronautics and Space Administration (NASA). This research made use of data products from the Midcourse Space Experiment. Processing of the data was funded by the Ballistic Missile Defense Organization with additional support from NASA Office of Space Science. This research has also made use of the NASA/ IPAC Infrared Science Archive, which is operated by the Jet Propulsion Laboratory, California Institute of Technology, under contract with the National Aeronautics and Space Administration. Finally, this paper made use of information from the Red MSX Source survey database at http://rms.leeds.ac.uk/cgi-bin/public/RMS_DATABASE.cgi which was constructed with support from the Science and Technology Facilities Council of the UK.

Facilities: 2MASS - , WISE - , Spitzer - , Herschel - , MSX. -

Software: astropy (Astropy Collaboration et al. 2013), Matplotlib (Hunter 2007), Pandas (McKinney 2010).

WISE Discovery of Mid-infrared Variability in Massive Young Stellar Objects

Article metrics

Permissions

Author e-mails

Author affiliations

ORCID iDs

Dates

Abstract

1. Introduction