A highly variable methanol maser in G111.256-0.770
MMNRAS , 1–7 (2019) Preprint 19 February 2019 Compiled using MNRAS L A TEX style file v3.0
A highly variable methanol maser in G111.256 − M. Durjasz (cid:63) , M. Szymczak and M. Olech
Centre for Astronomy, Faculty of Physics, Astronomy and Informatics, Nicolaus Copernicus University,Grudziadzka 5, PL-87100 Torun, Poland
Accepted XXX. Received YYY; in original form ZZZ
ABSTRACT
G111.256 − − ∼ > Key words: masers – stars: formation – stars: massive – radio lines: ISM – stars:individual: G111.256 − Methanol masers are one of the earliest observable signa-tures of high-mass young stellar objects (HMYSOs), lasting ∼ yr (van der Walt 2005) before or during the formationof ultra-compact H II regions (Walsh et al. 1998; Ellingsen2006; Urquhart et al. 2013; de Villiers et al. 2015). It isthought that they predominantly arise in the circumstellardiscs or outflows (e.g. Moscadelli et al. 2016; Sanna et al.2017) and one of the strongest transition at 6.7 GHz probesthe molecular gas of density 10 − cm − and temperaturebelow 150 K (Cragg, Sobolev & Godfrey 2005).Recent studies have revealed the maser lines as sensitiveindicators of sudden changes in the pumping conditions inthe environments of two HMYSOs likely triggered by accre-tion bursts (Caratti o Garatti et al. 2017; Moscadelli et al.2017; Hunter et al. 2018). Changes in the stellar luminositydue to accretion episodes appear to affect all maser featuresand the overall maser structure at 6.7 GHz in the HMYSOS255IR-IRS3 (Moscadelli et al. 2017) while turbulence in themolecular clouds may produce short-lived (a week to a fewmonths) random flares not correlated across different fea-tures (Goedhart, Gaylard & van der Walt 2004). When themaser beaming effect in an assembly of clouds of arbitrarygeometry is considered then quasi-periodic variations of themaser intensity may result from the rotation of maser cloudsacross the line of sight (Gray, Mason & Etoka 2018). Long-lived (several months to a few years) bursts may be causedby outflows or shock waves passing through the masing re-gion (Goedhart et al. 2004; Dodson, Ojha & Ellingsen 2004).In the paper we report the results a long term monitoring of (cid:63) E-mail: [email protected]
G111.256 − . × − L (cid:12) (Wu et al. 2010)spatially coincides with the 22 GHz H O masers and the ra-dio continuum source of spectral index 0.64 between 1.3 and3.6 cm wavelength is consistent with free-free emission froma thermal jet or with a partially thick H II region (Trinidadet al. 2006). Observations of outflow tracers revealed bothredshifted and blueshifted components that largely overlapon a scale lower than 10 − (cid:48)(cid:48) suggesting an outflow almostcompletely along the line of sight (Wu et al. 2010; S´anchez-Monge et al. 2013) while the H line emission knots and theradio continuum elongation delineate a collimated jet of size5 − (cid:48)(cid:48) at a position angle − ◦ (Massi et al. 2018). Most ofthe water maser components are blueshifted with respect tothe systemic velocity of − − (S´anchez-Monge et al.2013) and come from a region of size up to 800 − c (cid:13) a r X i v : . [ a s t r o - ph . S R ] F e b M. Durjasz et al. F l u x d e n s i t y ( J y )
43 42 41 40 39 38 37 36 35Velocity (kms )0.00.51.01.5 V a r i a b ili t y m e a s u r e red VIFI 0.0153045 r Figure 1.
Top:
Bottom:
Plots of variability index (VI), fluctuation index (FI) andreduced χ . For the analysis presented in the paper we used published(Szymczak et al. 2018) and unpublished archival 6.7 GHzmaser data and new observations, all of which were obtainedwith the Torun 32 m radio telescope. The new observationswere carried out from February 2013 to October 2018. Thebeam full width at half maximum of the antenna at 6.7 GHzwas 5.8 arcmin and rms pointing error was ∼
25 arcsec but itwas reduced to 10 arcsec from mid-2016 (Lew 2018). The sys-tem temperature was between 25 and 40 K. The data weredual polarization taken in frequency switching mode using a1 MHz shift. We used an autocorrelation spectrometer to ac-quire spectra with a resolution of 0.09 km s − and the typicalrms noise level of 0.35 Jy before 2015 May and 0.25 Jy af-terwards. The system was calibrated continuously against anoise diode of known constant temperature and this calibra-tion was daily checked by observing the non-variable masersource G32.744 − The average spectrum for the time-span from MJD 54583 to58418 (685 observations) together with the spectra at high-and low-emission levels and plots of three variability mea-sures are shown in Figure 1. The variability index (
V I ) andfluctuation index (
F I ) quantify the flux amplitude changein variable sources (Aller, Aller & Hughes 2003) and are de-scribed in Appendix A. All four persistent features displayedhigh variability over the whole observing period (Table 1).There are significant differences between
V I and
F I profiles which result from their properties;
F I is sensitive to changes in the flux density with low signal-to-noise ratioand its value is highest in the wings of the features, whereas
V I better depicts high amplitude variations.The light curve of the velocity-integrated flux densityis presented in Figure 2 where the new and archival ob-servations are complemented with the discovery detection(Szymczak et al. 2000) and a single VLA observation (Huet al. 2016). The emission peaked around MJD 55296 and56340 and the second outburst has a FHWM of ∼
256 d andthe rising phase a factor of two shorter than the decliningphase. In periods of low activity the emission linearly de-creased with a rate of − ± − yr − . The to-tal flux density ranged from 0.9 to 14.5 Jy km s − and theaverage value was 4.4 ± − . For the adopted dis-tance of 3.34 kpc the isotropic maser luminosity reached apeak of 1.0 ± × − L (cid:12) around MJD 56340 and declinedto 8.9 ± × − L (cid:12) at the end of the monitoring period.The median luminosity equals to 2.8 × − L (cid:12) is a factorof 3.2 lower than that reported for the sample observed ina high sensitivity, untargeted survey with the Arecibo tele-scope (Pandian et al. 2009). This confirms that the sourcebelongs to a population of weak methanol sources (Wu et al.2010).The dynamic spectrum in Figure 3 visualizes the bulkvariability and the bursting variability of 6.7 GHz maseremission. It was created using linear interpolation betweenconsecutive 32 m dish spectra and the emission above 3 σ level is shown. There are four spectral features in the spec-trum and all of them display complex and high variability(Table 1). In the following we deal with specific aspects of variabilityof the source.The Gaussian analysis of profiles revealed a strongblending effect in some velocity ranges at different epochs(Fig. 4). For instance the emission of middle velocity fea-tures from − − − is composed of two Gaus-sian components at two time intervals of MJD 54680 − − − − − and at − − − , respec-tively. Changes in the flux density and line width of thebursting features near − − − generallyfollowed that for the persistent emission with a mean veloc-ity of − − . From MJD 55418 to 57011 this persis-tent emission showed a drift in velocity of 0.071 km s − yr − .Figure 6 illustrates the behaviour of the emission near − − . The feature showed a long lasting ( ∼
700 d)burst with a peak of ∼ <
200 d) less visible bursts and it dropped belowour sensitivity level at the end of monitoring. No signifi-cant variations in the FWHM were seen. A velocity drift of − − yr − occurred before MJD 57039 and changedto − − yr − afterwards.Figure 7 depicts the variability of two features duringthe bursts. For the feature − − the rising phase ofa burst is clearly seen; the flux density increased from ∼ ∼ ∼
306 d, the peak velocity remained stablebut the FWHM decreased from 0.50 to 0.24 km s − . Thismay be attributed to unsaturated amplification despite the MNRAS , 1–7 (2019) ethanol variability in G111.256 − Table 1.
Variability measures of the 6.7 GHz features. Variability index (
V I ), fluctuation index (
F I ) and χ are for the mean velocityof feature while the average values < V I > , < F I > and < χ > and their standard errors are calculated for the given velocity range(∆ V ). V V I F I χ ∆ V < V I > < F I > < χ > (km s − ) (km s − ) − − − − − − − − − − − − I n t e g r a t e d f l u x d e n s i t y ( J y k m s ) MJD 51360
Figure 2.
Time series for the velocity-integrated flux density of the 6.7 GHz methanol maser line. The black squares and magentacircles denote the archival Torun 32 m data and from Szymczak et al. (2018), respectively. The blue triangle refers to the first detection(Szymczak et al. 2000) and the black star marks the VLA observation at MJD 55990 (Hu et al. 2016).
Figure 3.
False-colour image of the 6.7 GHz maser flux density versus velocity and time. Velocities are measured in the frame of thelocal standard of rest. Individual observation dates are indicated by the vertical bars below the top horizontal axis. The two longestintervals with no data are blanked.MNRAS000
False-colour image of the 6.7 GHz maser flux density versus velocity and time. Velocities are measured in the frame of thelocal standard of rest. Individual observation dates are indicated by the vertical bars below the top horizontal axis. The two longestintervals with no data are blanked.MNRAS000 , 1–7 (2019)
M. Durjasz et al. F l u x d e n s i t y ( J y )
43 40 37Velocity (kms )0.50.00.5 43 40 37 Figure 4.
Examples of Gaussian profile fits to the spectra takenat MJD 55203 (left) and MJD 56276 (right). V e l o c i t y ( k m s ) F W H M ( k m s ) F l u x d e n s i t y ( J y ) Figure 5.
Peak velocity, line width at half maximum (FWHM)and peak flux density for Gaussian components of the emissionnear − − . The red circles and green squares denotemaser components appeared at slightly different velocities. fact that the canonical relationship between the intensityand linewidth (Goldreich & Kwan 1974) is not fulfilled. Onthe other hand the declining phase is poorly seen and thereis no evidence for profile broadening after the maximum.The emission near − − showed a short (49 d)burst with rising and declining phases of 11 and 38 d, re-spectively. It is striking that the line width was nearly con-stant while the peak velocity drifted by 0.086 km s − dur-ing a ∼
31 d long burst. This phenomenon can appear whena factor triggering the burst propagates to nearby layers ofsimilar gas temperature but slightly different radial velocity.As the maser features are redshifted relative to the systemicvelocity the drift may trace inflow motion.We conclude, the variability of 6.7 GHz maser in thesource is characterised by short duration (2–5 months)bursts at slightly different velocities superimposed onlong lasting (1.5 − − V e l o c i t y ( k m s ) F W H M ( k m s ) F l u x d e n s i t y ( J y ) Figure 6.
Same as in Fig. 5 but for the emission near − − . The green horizontal lines refer to the mean valueof FWHM, the dashed violet line denotes the detection level of0.8 Jy. The blue (dashed) and red (dotted) lines mark 1 σ and 3 σ levels, respectively. V e l o c i t y ( k m s ) F W H M ( k m s ) F l u x d e n s i t y ( J y ) Figure 7.
Same as in Fig. 5 but for two exemplary bursts. Leftand right panels refer to − − − features, re-spectively. The green vertical lines denote the start time of burstsand the red vertical marks the end of burst. The grey symbolsrepresent the measurements with the flux density below 3 σ level. changes in the flux density showing velocity drifts (up to ∼ − yr − ) for some features. In order to determine a variability timescale we used thediscrete structure function and minimum-maximum method(see Appendix A, Fuhrmann et al. 2008). The structure func-
MNRAS , 1–7 (2019) ethanol variability in G111.256 − (d)10 S t r u c t u r e f un c t i o n -41.3 kms -38.7 kms -37.8 kms -37.0 kms Figure 8.
Structure functions of the persistent maser features.
Table 2.
Variability timescales for the spectral features. t var and t sf are the timescales calculated with the minimum-maximummethod and the structure function, respectively.Feature ( km s − ) t var (d) t sf (d) − ± +60 . − . − ± +58 . − . − ± +90 . − . − ± +198 . − . tions for the four persistent features are presented in Figure8 and values of variability timescales obtained with bothmethods are given in Table 2. The timescales for the features − − − are 1.5 to 1.7 yr. The uncertain-ties of these estimates are usually lower than 26 per centwith the exception of the − − feature that showsonly a monotonic decrease and t sf could not be accuratelydetermined.The values of the variability timescales for features at − − − obtained by the two methods dif-fer significantly (Table 2). This is probably due to the factthat t var strongly relies on the individual choice of the min-imum and maximum measurements that make it irrelevantfor describing the long lasting burst variability. We concludethat the variability timescales of the methanol masers inG111 range from 0.7 to 1.8 yr.In order to search for time lags we computed the dis-crete cross-correlation function (DCF, Edelson & Krolik1988) between the three spectral features. The maximumvalue was obtained with the centroid τ c of the DCF givenby τ c = (cid:80) i τ i DCF i / (cid:80) i DCF i . Following the method pre-sented in detail by Peterson et al. (1998) and Fuhrmannet al. (2008) we performed Monte Carlo simulations to sta-tistically estimate meaningful time lags and their uncertain-ties. The influence of uneven sampling and the measurementerrors were taken into account by using random subsets ofthe two spectral features time series and by adding randomGaussian fluctuations constrained by the measurement er- C o rr e l a t i o n C o e ff i c e n t
180 155 130Lag (d) 0.050.150.25 P r o b a b ili t y Figure 9.
Left.
Example of discrete correlation function betweenthe light curves of features − − − . Right.
Cross-correlation peak distribution for the peak of the DCF.
Table 3.
Time lags corresponding to the average centroids of theDCF peaks obtained from cross-correlation peak distributions.Uncertainty values correspond to 1 σ obtained from CCPDs.Features ( km s − ) τ (d) − − − ± − − ± rors. For two pairs of the spectral features ( − − vs − − and − − vs − − ) thecentroid of the DCF maximum was computed 1000 times.The DCF and cross-correlation peak distribution (CCPD)for one pair of features are shown in Figure 9. Table 3 sum-marizes the results for two pairs of features. The time de-lay of the peaks estimated between the two maser features( − − − ) in the timespan from ∼ MJD56200 to 56900 is 5 months. The delay of the peak of theredshifted emission ( − − ) has a too large uncer-tainty to be considered as reliable. A visual inspection ofthe light curves (Fig. 3) implies that the burst of feature − − around MJD 56200 was advanced by 147 drelative to the burst of feature at − − . This crudeestimation is in good agreement (within 20 per cent) withthat presented in Table 3. Our monitoring revealed that the methanol emission fromG111 has experienced complex and erratic variations for thelast decade. In the following we discuss the observed char-acteristics which may shed light on the processes causingvariability.One of the important factors producing the maser vari-ability can be changes in the pump rate. Since the 6.7 GHztransition is thought to be pumped by infrared photons(Sobolev et al. 1997; Cragg et al. 2005) we examine if themaser emission is related to the near infrared radiation.
MNRAS000
MNRAS000 , 1–7 (2019)
M. Durjasz et al.
51 0 5 5 0 0 0 5 5 5 0 0 5 6 0 0 0 5 6 5 0 0 5 7 0 0 0 5 7 5 0 0 5 8 0 0 0 024 ( a )
W 1 m e a s u r e dW 2 m e a s u r e d
IR relative flux
M J D m e a nm e a n
I R r e l a t i v e f l u x
W 2 N E O W I S EW 2 W I S Er = 0 . 6 5 2 ( p < 0 . 0 2 9 8 ) ( b )
W 1 N E O W I S EW 1 W I S Er = 0 . 8 5 1 ( p < 0 . 0 0 0 8 8 )
Figure 10. ( a ) Part of the maser light curve from Fig. 2 withthe IR relative fluxes superimposed. The maser relative velocity-integrated flux density and its 5 point averaged values are markedby grey circles and grey line, respectively. IR data are from WISEand NEOWISE observations at bands W1 (3 . µ m, blue crosses)and W2 (4 . µ m, red crosses). The corresponding relative averageIR fluxes are denoted by blue circles and red squares, respectively.( b ) Relationships between the maser and IR relative fluxes. Thesolid lines denote the linear least-square fits to the data. Data from the WISE and NEOWISE archives (Wright et al.2010; Mainzer et al. 2011) were retrieved for the target toconstruct the light curves at 3.4 and 4.6 µ m. In Figure 10awe show an overlay of the relative IR light curves with themethanol maser light curve. The WISE and NEOWISE pho-tometric observations were available for 11 sessions of lengthof 1 to 8 d (median 1.8 d) over time-span 2010 − > ∼ MJD 55200 − − − before ∼ MJD 57000 can be easily recognizedas due to blending effect of two features showing temporalintensity variations because velocity jumps were clearly seenover certain timespans (Fig. 5). A similar effect may occurfor the − − feature with FWHM ranging from 0.30to 0.45 km s − over the monitoring period (Fig. 6) and theslope of the drift changed to steeper after ∼ MJD 57000. Inthis case we cannot resolve blending components probablydue to the low signal-to-noise ratio. The opposite sense of thevelocity drifts reinforces the interpretation that the velocityof emission fluctuates rather than drifting steadily. One canconclude that the velocity drifts are caused by slow ( <
10 yr)variations in the flux density of a few features at very closevelocities.In G111 we sometimes observed correlated variationsacross two different features combined with short-lived(2 − − . It very unlikely that the methanol maser wouldsurvive such a shock. Random bursts indicate rather localchanges in the maser conditions modified for instance byturbulence (Sobolev et al. 1998); a dispersion of turbulencevelocity of about 1 km s − could account for the observedshort-lived bursts caused by changes in the velocity coher-ence. We report new 6.7 GHz methanol maser observations ofG111 obtained since 2013 which extend the previously pub-lished light curve to ∼
11 yr. The maser emission is charac-terised by short duration (2-3 months) bursts superimposedon long lasting ( > ACKNOWLEDGEMENTS
We thank the Torun CfA staff for assistance with the obser-vations. We appreciate Eric G´erard for carefully reading themanuscript and providing critical comments. This researchhas made use of the SIMBAD data base, operated at CDS(Strasbourg, France), as well as NASA’s Astrophysics DataSystem Bibliographic Services. This work has also made useof data products from the Wide-field Infrared Survey Ex-plorer, which is a joint project of the University of California,Los Angeles, and the Jet Propulsion Laboratory/CaliforniaInstitute of Technology, funded by the National Aeronauticsand Space Administration and from NEOWISE, which is aproject of the Jet Propulsion Laboratory/California Insti-tute of Technology, funded by the Planetary Science Divi-sion of the National Aeronautics and Space Administration.
MNRAS , 1–7 (2019) ethanol variability in G111.256 − The work was supported by the National Science Centre,Poland through grant 2016/21/B/ST9/01455.
REFERENCES
Aller M. F., Aller H. D., Hughes P. A., 2003, ApJ, 586, 33Bartkiewicz A., Szymczak M., van Langevelde H. J., 2016, A&A,587, A104Caratti o Garatti A., et al., 2017, Nature Physics, 13, 276Choi Y. K., Hachisuka K., Reid M. J., Xu Y., Brunthaler A.,Menten K. M., Dame T. M., 2014, ApJ, 790, 99Cragg D. M., Sobolev A. M., Godfrey P. D., 2005, MNRAS, 360,533Dodson R., Ojha R., Ellingsen S. P., 2004, MNRAS, 351, 779Edelson R. A., Krolik J. H., 1988, ApJ, 333, 646Ellingsen S. P., 2006, ApJ, 638, 241Fuhrmann L., et al., 2008, A&A, 490, 1019Goddi C., Moscadelli L., Alef W., Tarchi A., Brand J., Pani M.,2005, A&A, 432, 161Goedhart S., Gaylard M. J., van der Walt D. J., 2004, MNRAS,355, 553Goldreich P., Kwan J., 1974, ApJ, 190, 27Gray M. D., Mason L., Etoka S., 2018, MNRAS, 477, 2628Heeschen D. S., Krichbaum T., Schalinski C. J., Witzel A., 1987,AJ, 94, 1493Hu B., Menten K. M., Wu Y., Bartkiewicz A., Rygl K., ReidM. J., Urquhart J. S., Zheng X., 2016, ApJ, 833, 18Hunter T. R., et al., 2018, ApJ, 854, 170Lew B., 2018, Experimental Astronomy, 45, 81Mainzer A., et al., 2011, ApJ, 731, 53Massi F., Moscadelli L., Arcidiacono C., Bacciotti F., 2018, inTarchi A., Reid M. J., Castangia P., eds, IAU SymposiumVol. 336, Astrophysical Masers: Unlocking the Mysteries ofthe Universe. pp 289–290, doi:10.1017/S1743921317011644Moscadelli L., et al., 2016, A&A, 585, A71Moscadelli L., et al., 2017, A&A, 600, L8Pandian J. D., Menten K. M., Goldsmith P. F., 2009, ApJ, 706,1609Peterson B. M., Wanders I., Horne K., Collier S., Alexander T.,Kaspi S., Maoz D., 1998, PASP, 110, 660S´anchez-Monge ´A., L´opez-Sepulcre A., Cesaroni R., WalmsleyC. M., Codella C., Beltr´an M. T., Pestalozzi M., MolinariS., 2013, A&A, 557, A94Sanna A., Moscadelli L., Surcis G., van Langevelde H. J.,Torstensson K. J. E., Sobolev A. M., 2017, A&A, 603, A94Sobolev A. M., Cragg D. M., Godfrey P. D., 1997, A&A, 324, 211Sobolev A. M., Wallin B. K., Watson W. D., 1998, ApJ, 498, 763Szymczak M., Hrynek G., Kus A. J., 2000, A&AS, 143, 269Szymczak M., Olech M., Sarniak R., Wolak P., Bartkiewicz A.,2018, MNRAS, 474, 219Trinidad M. A., Curiel S., Torrelles J. M., Rodr´ıguez L. F., Mi-genes V., Patel N., 2006, AJ, 132, 1918Urquhart J. S., et al., 2013, MNRAS, 431, 1752Walsh A. J., Burton M. G., Hyland A. R., Robinson G., 1998,MNRAS, 301, 640Wright E. L., et al., 2010, AJ, 140, 1868Wu Y. W., Xu Y., Pandian J. D., Yang J., Henkel C., MentenK. M., Zhang S. B., 2010, ApJ, 720, 392de Villiers H. M., et al., 2015, MNRAS, 449, 119van der Walt J., 2005, MNRAS, 360, 153
APPENDIX A: ANALYSIS METHODSA1 Statistical variability measures
We used the variability index as proposed by Aller et al.(2003)
V I = ( S max − σ Smax ) − ( S min − σ Smin )( S max − σ Smax ) + ( S min − σ Smin ) , (A1)where S max and S min refer to the maximum and minimumflux densities, respectively, while σ Smax and σ Smin are thecorresponding uncertainties of S max and S min . This index isuseful to estimate the amplitude of the variability accountingfor the measurement uncertainties and is well determinedwhen variability is much greater than measurement errors.The second variability measure used is the fluctuationindex (Aller et al. 2003) F I = (cid:20)(cid:18) (cid:80) S w i − S (cid:80) S i w i N − − (cid:19) N (cid:80) w i (cid:21) . S , (A2)where N is the number of observations of the flux density S i measured with error σ i , weight is w i = σ − i and S is theaverage flux density. This index well estimates variability offeatures with low signal-to-noise ratio.We also examined variability by computing the reducedvalue of χ . χ = 1 N − N (cid:88) i =1 (cid:18) S i − Sσ i (cid:19) (A3) A2 Time scales of variability
We calculated the discrete structure function SF ( τ j ) = n (cid:88) i =1 [ S ( t i ) − S ( t i + τ j )] (A4)for each data set S ( t i ) following the procedure described inHeeschen et al. (1987). Here, n is the number of flux densitymeasurements obtained at epochs t i and t i + τ j . We binnedthe data in 31 d intervals to get evenly sampled dataset.Two linear fits were performed to estimate the slope ( aτ β )of the SF before the first maximum and the mean value of SF after the first maximum which estimates the saturationlevel ρ . The timescale of variability was then determinedfrom formula t sf = (cid:16) ρ a (cid:17) β (A5) This paper has been typeset from a TEX/L A TEX file prepared bythe author.MNRAS000