Seven-year periodic variations in the methanol maser line displayed by the massive protostar IRAS 20216+4104 ⋆

,


Introduction
Astrophysical masers are known to be very sensitive to changes in the physical conditions in their host clouds, and primarily those caused by radiation and collision processes.For class II methanol masers, diverse variability patterns occur from very complex alteration of the spectra with the appearance and disappearance of features, bursts of individual features and their groups on a timescale of a few days to years to a steady rise and/or decline in feature peaks over several years (e.g., Goedhart et al. 2004;Szymczak et al. 2018).While most methanol maser sources do not show predictable variability patterns, there is a growing group of 6.7 GHz methanol maser sources with strictly periodic variability.About 30 sources with period ranging from 24 to 1260 d with a median of 220 d have been discovered so far (Araya et al. 2010;Goedhart et al. 2003Goedhart et al. , 2004Goedhart et al. , 2009Goedhart et al. , 2014;;Fujisawa et al. 2014;Maswanganye et al. 2015Maswanganye et al. , 2016;;Olech et al. 2019Olech et al. , 2022;;Proven-Adzri et al. 2019;Sugiyama et al. 2015Sugiyama et al. , 2017;;Szymczak et al. 2011Szymczak et al. , 2015Szymczak et al. , 2016;;Tanabe et al. 2023).
The case of a wide span of periods and diversity of light curves from sinusoidal-like to intermittent bursts with a relative amplitude of 0.2 to more than 120 has brought on competing hypotheses on the nature of the maser periodicity phenomenon.van der Walt (2011, see also van der Walt et al. 2016) proposed changes in the seed photon flux due to the orbital modulation of the free-free emission of the hot post-shock gas in a ⋆ The reduced spectra are available at the CDS via anonymous ftp to cdsarc.cds.unistra.fr(130.79.128.5) or via https:// cdsarc.cds.unistra.fr/viz-bin/cat/J/A+A/682/A17colliding-wind binary (CWB) system as a primarily periodicity mechanism.This simple model satisfactorily explains the maser light curves with an asymmetric flare profile, namely, with a fast rise and slow decay in flux density, which is observed in about a half of the periodic 6.7 GHz masers (Szymczak et al. 2015).A more sophisticated model of a CWB aptly describes the flare profile of several sources and light curves over a decade (van den Heever et al. 2019).
Methanol maser periodicity can be driven by luminosity variations in a binary system with a massive primary (Araya et al. 2010;Parfenov & Sobolev 2014) or the pulsation of a single massive star (Inayoshi et al. 2013), which affect the dust temperature and, consequently, the pumping rate.Luminosity variations can be due to periodic modulations of accretion rate in binary systems (e.g., Artymowicz & Lubow 1996;Muñoz & Lai 2016;Parfenov & Sobolev 2014).In contrast to the CWB model, these mechanisms do not predict a specific flare profile or light curve.Gray et al. (2020) performed 3D model simulations for a maser cloud, ultimately concluding that a flare profile driven by a change in pump rate is qualitatively different from that caused by background variation.After applying the observed IR light curves of two sources as the driving function in their model, they rejected the variable background mechanism for the flaring 6.7 GHz methanol masers.
A scarcity of long (>2 yr) period masers may be due to the shortness of monitoring observations, which usually lasts 3-4 yr (Goedhart et al. 2004;Szymczak et al. 2018) because in those programmes, a periodicity that is on timescales of less than 220 d could be sampled effectively.The recent finding of 1260 d periodicity in G5.900−0.430has resulted from a 10 yr monitoring (Tanabe et al. 2023).Here, we report the 6.7 GHz maser periodicity in IRAS 20126+4104 that is longer by a factor of 2.
IRAS 20216+4104 is one of the most frequently studied highmass young stellar objects (HMYSOs) with the disc-outflow system (Cesaroni et al. 2007 for review).The central protostar has a mass of ∼12M ⊙ (Cesaroni et al. 2014(Cesaroni et al. , 2023;;Chen et al. 2016) and a luminosity of ∼10 4 L ⊙ (Johnston et al. 2011;Cesaroni et al. 2023) at distance of 1.64 kpc (Moscadelli et al. 2011).A resolved accretion disc traced by several molecular lines undergoes Keplerian rotation (Cesaroni et al. 2005), but its structure has no axial symmetry in the CH 3 CN(12-11) K = 8 line emission, probably due to the effect of nearby stellar companions (Cesaroni et al. 2014;Massi et al. 2023).Multi-epoch photometry of the outflow cavities reveals possibly periodic variability of the NIR continuum that may be due to the wobbling of the inner disc caused by a nearby low-mass star companion (Massi et al. 2023).The water maser emission at 22 GHz traces a bipolar outflow of velocity up to ∼100 km s −1 (Moscadelli et al. 2011).A sparsely sampled light curve of the 22 GHz water maser from 1989 to 2007 is reported in Felli et al. (2007); it does not show any apparent periodicity or trends.The 6.7 GHz methanol maser emission does not delineate any disc-like structure; it has appeared only in the northern parts of the disc at 390-530 au from a central star whose position was derived from the water maser modelling (Moscadelli et al. 2011;Cesaroni et al. 2013).The methanol maser at 6.7 GHz was observed at two non-overlapping 3-4 yr intervals (Goedhart et al. 2004;Szymczak et al. 2018).Goedhart et al. (2004) reported that the source may be variable, but the signal-to-noise ratio was low.We note, however, a weak variability of −6.1 km s −1 feature in both studies.Here, we use another data set to get a longer light curve.

Data
The HartRAO data used here were taken as part of a monitoring programme published in Goedhart et al. (2004), spanning from January 2000 to February 2003 (MJD 51567−52681).The target was observed on average twice a month with a spectral resolution of 0.112 km s −1 and 3σ sensitivity of 1.2 Jy.We obtained data from the Torun radio telescope between June 2009 and September 2023 (MJD 55006−60208) at a weekly cadence.The data from the initial period (June 2009 to February 2013) were published in Szymczak et al. (2018), where the instrument, observations, and data reduction methods were described in detail.In short, a spectral resolution was 0.09 km s −1 , a typical 3σ sensitivity of 0.9 Jy and calibration accuracy of ∼10% In addition, we searched for observations of the source in the literature and found eight data points spanning from October 1991 until June 2016.The dates, velocity-integrated flux density and instruments are listed in Table A.1.Thus, the time range of observations of the target was extended to 31.9 yr.Since these spectra were obtained with different sensitivity and resolution, we used the series of integrated flux densities across the entire velocity range and time interval.We used more homogeneous data obtained with the Torun telescope to analyze the light curve of individual features.observational data, scattered over time and from different independent instruments, are consistent with the pattern of variability displayed since 2009.The Lomb-Scargle periodogram (Fig. 2) indicates that the period of variability is 2520±70 d, the uncertainty is estimated as the half width at half maximum of the Gaussian function fitted to the top half of the fundamental peak.

Results
Since June 2009 (MJD 55006), the object was observed weekly with four gaps lasting more than two months.The overall picture of variability over 14.2 yr is shown in Fig. A.1.There is a clear indication (Fig. 1) that the average flux in the last four maxima increases by about 15% in each successive cycle, while at the minima, the flux is probably constant.Similar changes of the flare amplitude fitted with a low-degree polynomial were observed, over 3-8 yr, for some spectral features of periodic (0.4−0.6 yr) masers such as G30.400−0.296and G108.76−0.99 (Olech et al. 2019).
The average and typical spectra in the high and low emission states as well as the normalized χ 2 plot are displayed in Fig. 3. Significant differences exist in the spectra shapes in the low and high states.The light curves of seven features are presented in Fig. 4 and their variability properties are listed in Table 1.Except for the −6.1 km s −1 feature, all six others show a sinusoidal-like A17, page 2 of 9  variability with a plateau minimum equally lasted 1.5±0.2yr in both cycles.In the plateau state, the flux density dropped below our detection limit (∼0.9 Jy) for features −5.57and −4.96 km s −1 .The relative amplitude of these six features varies from 2 to >89.The index of their variability estimated with the χ 2 r is significant, while the −7.68 km s −1 feature is highly variable; its intensity changed by a factor of 6-9 on timescales of 1.5-2.0yr.The feature −6.10 km s −1 shows relatively weak variability, and its intensity is anti-correlated with all the rest of the features.Figure 5 compares the light curves of the most and most minor variable features at −7.68 and −6.10 km s −1 , respectively.There is moderate anti-correlation with the −0.55 coefficient for the entire period, but when the interval is limited to around MJD 55600-57100, the coefficient improves to −0.63 then declines to −0.45.A few episodes of even stronger anti-correlation lasting longer than one month occurred.All but −6.10 km s −1 features show quasi-periodic (∼3 weeks) variations in the flux density with the amplitude higher than 20-40%.There was no time lag between spectral components within an accuracy of 5 days.

Discussion
According to most hypotheses, periodic variations of the maser intensity occur in binary systems with a massive star.In such Notes.V p is the peak velocity, S p is the average flux density, R a is the relative amplitude, χ 2 r is the reduced χ 2 , and χ 2 99.9% is the corresponding value for 99.9% significance level of variability.systems, there is a well-defined periodicity with which the interaction of stellar winds near a periastron raises the level of radio background radiation, which may increase the number of seed photons (van der Walt 2011; van der Walt et al. 2016;van den Heever et al. 2019).In close binaries, the accretion rate can periodically vary (Artymowicz & Lubow 1996;Muñoz & Lai 2016), causing changes in stellar luminosity which modulate the dust temperature and affect the pumping rate (Araya et al. 2010;Parfenov & Sobolev 2014).
A sinusoidal light curve prompts us to consider explaining the cause of variability due to the pulsation of a massive protostar with a high accretion rate (Inayoshi et al. 2013).For a 2520 d period, their model implies a luminosity of the central star in IRAS 20216+4104 of ∼10 6 L ⊙ , that is two orders of magnitude higher than what was measured from the spectral energy distribution fit (e.g., Cesaroni et al. 2023).Therefore, we can rule out this hypothesis.
Spectroscopic and imaging optical observations have shown that 70% to 90% massive stars have at least one companion (Sana et al. 2014;Lanthermann et al. 2023).A survey of radial velocities of 16 embedded high-mass stars revealed that at least 20% of the targets are close binaries (Apai et al. 2007).IRAS 20216+4104 is probably no exception, as suggested by A17, page 3 of 9 some hints.X-ray observations revealed that the source is associated with an embedded, possibly very young stellar population (Montes et al. 2015).Shepherd et al. (2000) suggested that the jet precession can be due to tidal interactions between the disk and a companion in a non-coplanar orbit.Interferometric maps of molecular lines implied the Keplerian pattern of the disk.Still, the high energy (513 K) CH 3 CN emission presents two distinct peaks, which one located slightly offset from the centre may be due to a lower mass companion orbiting a high-mass central star (Cesaroni et al. 2005).The authors have speculated that the observed axial asymmetry of the disc (Cesaroni et al. 2014) may be caused by a nearby, unseen, low-mass (2 M ⊙ ) companion.
The recent finding of anti-correlated variability of the NIR continuum emission of the N-E and S-E outflow cavities gives another hint at a nearby low-mass companion which may cause the wobbling of the inner disc that alternatively eclipses the two hemispheres of the central massive star (Massi et al. 2023).Six-epoch NIR photometry suggests a possible periodicity of 12−18 yr (Massi et al. 2023, their Fig. 4); we note, however, that the NIR continuum emission of a large part (S1 and S2 polygons in their Fig. 2) of the S-E outflow cavity closely follows ∼7 yr periodic behaviour of the 6.7 GHz methanol maser line (Fig. 1).In contrast, the emission from the neighbouring region (polygon PD) nearest to the central star shows no similar variation.The brightness variations of the N-W lobe are anti-correlated with the maser intensity with a possible increase of the degree of anti-correlation for the farthest area (polygon N3).As the data are sparse, we can only hypothesize that the NIR continuum emission from the farthest (≥10 000 au) portions of lobes is more strongly (anti-) correlated with the maser emission.
To further investigate the source variability, we used the W1 (3.4 µm) and W2 (4.6 µm) observations from the NEOWISE (Mainzer et al. 2014) single exposure database 1 .The high-quality (ph_qual = AA) exposures are only used to obtain the average fluxes at eighteen epochs since May 2014.The W1 and W2 average magnitudes of each epoch superimposed with the maser integrated flux density are shown in Fig. 6.Significant variability exists with the relative amplitude of 0.28 and 0.37 at 3.4 µm and 4.6 µm wavelengths, respectively.There is a time lag of 0.12/0.20Pbetween the light curve maxima and minima in the 1 https://irsa.ipac.caltech.edumaser and NIR bands, where P denotes the period.High angularresolution observations have revealed that the maser regions are 390-530 au away from the central star (Moscadelli et al. 2011;Cesaroni et al. 2013;Aberfelds et al. 2023).This could imply a velocity of 1300-3000 km s −1 for the propagation infrared radiation, which drives the variations of the maser emission and is one order of magnitude lower than that postulated for mal induced by events (Burns al. 2020).Therefore, it is less probable that bursts of the central star the maser variability in IRAS 20216+4104.The maser curve of our appears to be similar that reported in G331.13−0.24,but with a 509 d period (Goedhart et al. 2004(Goedhart et al. , 2014) ) and is not consistent with a simple CWB model (van der Walt 2011).
Intriguingly, the NEOWISE NIR light centroids show a timedependent shift of about 5. ′′ 5 (Fig. 7) that corresponds to 9000 au, namely, it is equal to the separation of the brightest inner spots of outflow lobes seen in K s band (Cesaroni et al. 2013).Furthermore, the position angle of this displacement of −62 • is in perfect agreement with the outflow axis (Cesaroni et al. 2013).Epochs 1-3 of NEOWISE observations were likely at maximum NIR radiation and in the last part of the rising phase of the maser, the centroids are located in the middle of the plot.The NIR light centroids started to shift in the S-E direction after the maser maximum, reaching the extreme displacement at epochs 5-7, and then going back to the centre during the decay phase of the maser light curve.In the plateau minimum of the integrated maser emission (epochs 10-14), the centroids were shifted toward the northwest, and then when the maser emission increased in epoch 15, the centroids returned to the central position and started to move toward the southeast.In general, the observed shifts of NIR emission centroids can be an effect of light echo, triggered by luminosity bursts, in outflow cavities oriented outside the plane of the sky (Caratti o Garatti et al. 2017;Stecklum et al. 2018).In the case of IRAS 20216+4104, the jet or outflow average inclination angle to the plane of the sky is as small as 8±1 • (Massi et al. 2023), so the light echo effect is unlikely.We suggest that periodic luminosity variations of the powering protostar might not drive the methanol maser variability in the target.
The waveform light curves were observed in G338.92−0.06(Goedhart et al. 2004(Goedhart et al. , 2014) ) and G358.460−0.391(Maswanganye et al. 2015) and interpreted as an eclipsing A17, page 4 of 9 Szymczak, M., et al.: A&A, 682, A17 (2024) Fig. 8. Map of the maser spots in IRAS 20216+4104 obtained on 2012 March 4 with VLA-C (Hu et al. 2016).The colour and size of symbols correspond to the velocity scale, as shown in the wedge and the logarithm of intensity, respectively.The circles and squares denote the maser spots whose emission synchronously varies with the NIR continuum of the SE and NW lobes (Massi et al. 2023), respectively.The dashed arrows show the direction of the large-scale (∼10 ′′ ) outflow (Cesaroni et al. 2013).The asterisk marks the nominal position of the protostar derived from the model fit to the water masers (Moscadelli et al. 2011).The cross marks the position of the 230 GHz continuum emission peak and its size corresponds to the uncertainty (Cesaroni et al. 2014).The dotted black line marks the direction of the disc plane (Cesaroni et al. 2005).The dotted magenta line is orthogonal to the axis of large-scale jet or outflow and divides the maser spot distribution into two parts, which show anti-correlated variability.effect because they very much resemble that of some eclipsing binary systems (Maswanganye et al. 2015).They have a symmetric shape around the minima, which lasts 0.04-0.09P.Data from the Torun telescope indicate that the curve of IRAS 20216+4104 has a flat minimum lasting about 0.22P.During the minima, the red-shifted emission at a velocity higher than −5.6 km s −1 drops below 3σ detection level, while the emission at lower velocity is commonly above the sensitivity threshold (Fig. A.1).In general, there is a tendency for the blue-shifted (≤ −5.6 km s −1 ) emission to come from the NE and N sides of the distribution (Fig. 8).The flat and long minima of the emission from these regions are not affected by our sensitivity, in contrast to the red-shifted (≥ −5.6 km s −1 ) emission from the southern region, which was detected only with the VLA-C (Hu et al. 2016), but generally unseen in VLBI maps (Aberfelds et al. 2023).
A plausible explanation can be that an eclipse of the central star by an inner wobbling disk can modulate the scattered light in the outflow lobes.Moreover, this eclipsing can vary the maser transition's pumping rate by cyclic shading portions of the dust disk and indirectly affect the maser intensity.The inner disk cannot eclipse the maser itself since it is 390−530 au from the massive star in the N-E part of the gas disk (Fig. 8).The shadowed portions rotate over the dust disk surface and prograde with wobbling disk motion.A simple visualisation shows that a maximum of NIR radiation seen by the observer occurs for PA ≃ 57 • of the major axis of the wobbling disk (Fig. A.2b).Since the disk rotates counterclockwise (Cesaroni et al. 2014), this maximum of NIR radiation seen by NEOWISE precedes the maximum pumping radiation seen through the maser region (Fig. A.2c).Following the discussion in Massi et al. (2023), we refine the system's parameters.Assuming the masses of primary and secondary protostars are 12M ⊙ and 2M ⊙ , respectively (Chen et al. 2016;Cesaroni et al. 2014) and with the secondary's orbit inclination of 30 • relative to the disk plane and orbital radius of 14 au, we obtain a disk radius of 12 au.Such values fall well within the ranges observed for pre-main sequence stars (Terquem et al. 1999).Interpretations suggesting the periodic variability of the maser emission is attributed to an eclipsing effect are supported by a characteristic feature of the maser light curve in IRAS 20216+4104, namely: the rising and falling edge adjacent to the minimum is relatively symmetric standard for eclipsing binary systems.
The maser emission observed with the VLA-C in March 2012 (Hu et al. 2016) is elongated in the NE−SW direction with a velocity gradient (Fig. 8) and it partially overlaps with the NE part of the gas disk of ∼1.′′ 5 in size, as traced by the thermal lines of CH 3 CN at 230 GHz (Cesaroni et al. 2014).The extremely northwesterly components around −6.1 km s −1 are marked as the squares in Fig. 8, lying about 90 au away from the rest of the maser clouds; they do not follow a general trend in terms of their variations.This is probably because they could be located in the outer regions of the jet as indicated a proper motion study (Moscadelli et al. 2011).Consequently, these regions may not have been exposed (or only weakly) to eclipse effects by the wobbling of the inner disk.

Conclusions
In this paper, we characterize the long-term variability of the 6.7 GHz methanol maser emission of the high-mass protostar IRAS 20126+4104 using both new and archival data from a 31 yr period.The flux density of all but one spectral feature shows a variation with a period of 2520 days (∼7 yr).This is the longest observed period, exceeding twice the period of the previous known longest periodic source.The maser variations strongly correlate with the NIR continuum emission of the southeastern outflow cavity measured six times over 20 yr by Massi et al. (2023).Data from 18 NEOWISE observations over 8.5 yr indicate variations of the IR emission, whose peak precedes by one-fifth of the period the maser maximum.The eclipsing effect of the central star by a wobbling inner disk caused by a nearby orbiting low-mass companion can possibly explain the characteristics of both maser and infrared variability.
This serendipitous finding of long-term periodic variability of the methanol maser in the extensively studied massive star with the disk or jet system could trigger multi-band monitoring over a flare cycle to constrain maser models as well as contribute to our understanding of binary and multiple systems with a massive primary.

Figure 1 Fig. 1 .Fig. 2 .
Figure1shows the maser light curve of IRAS 20216+4104 over 31.9 yr and the best fit with a sinusoidal function.Noticeably, the

Fig. 3 .
Fig. 3. Comparison of 6.7 GHz maser spectra obtained with the Torun telescope.The average spectrum is for 669 observations from MJD 55006 to 60208.The spectra marked by the cyan and red lines are typical for high (H) and low (L) activity states, respectively.The black line represents the normalized χ 2 values as a measure of variability.The dashed vertical line indicates the systemic velocity(Cesaroni et al. 2005).

Fig. 4 .
Fig. 4. Light curves of the main maser features in IRAS 20216+4104 taken with the Torun telescope.The features are labeled with the peak velocity.The flux density scale of the subsequent curves is shifted by 10 Jy.

Fig. 5 .
Fig. 5. Light curves of the −7.68 and −6.10 km s −1 features.The thick lines show the average curves.The magenta boxes mark the intervals of anti-correlated short-term variability.The inset shows the anticorrelation of these features' flux density, where the symbol's colour denotes the time since MJD 55000.

Fig. 6 .
Fig. 6.NEOWISE light curves at 3.4 µm (W1) and 4.6 µm (W2) of the target superimposed with the 6.7 GHz maser velocity-integrated flux density.The colours of circles and squares correspond to the average date of NEOWISE exposures in 18 numbered epochs.The bars show the standard deviation.The thick black line shows the average maser light curves.

Fig. 7 .
Fig. 7. Offsets from the overall average position for 18 epochs of NEOWISE observations of IRAS 20216+4104.The numbered symbols correspond to the epochs in Fig. 6.The bars indicate the standard errors of the mean.

Fig. A. 2 .
Fig. A.2. Sketch of variability scenario in IRAS 20216+4104.The inner part of the disk of massive protostar wobbles due to the interaction of a nearby low-mass companion with the orbital period of ∼14 yr and alternatively eclipses the two hemispheres of the central protostar.The eclipse modulates the scattered light in the outflow lobes (blue and red arrows) and the dust temperature, causing the maser (green ellipse) intensity variations.The shades of grey encode the relative irradiance of the outer disk.The orientation of the major axis of the disk and outflow lobes and the maser location is from the observations (see the caption of Fig. 8), but the sizes of the model components are not shown to scale.

Table 1 .
Variability characteristics of the spectral features.