© The Authors 2026

1 Introduction

The existence of variable sources in high-mass star-forming regions has been established in the past through the detection of regular and irregular changes in emission at radio and infrared wavelengths. Monitoring of water maser emission revealed increases in intensity by orders of magnitude on timescales as short as a few days (see e.g. Felli et al. 2007 and references therein). Although more stable, class II methanol masers were also found to have explosive behaviour in a handful of objects (Fujisawa et al. 2015; Hunter et al. 2017; Sugiyama et al. 2019; Hirota et al. 2022), but more than 30 sources are known to undergo regular variations with periods ranging from ~24 days to ~4.4 yr (see Tanabe et al. 2024 and references therein). Sometimes variations in IR intensity have been observed in association with increases in CH₃OH maser flux density (Stecklum et al. 2016; Garatti et al. 2017; Harajiri et al. 2026), as expected if the inversion of the level populations is due to radiative pumping.

Recently, the Near-Earth Object Wide-field Infrared Survey Explorer (NEOWISE) mission has provided the astronomical community with regular monitoring of the sky at 3.4 µm and 4.6 µm over a decade. These data triggered searches for variable young stellar objects (YSOs), establishing that variability is quite common over timescales from days to decades (Park et al. 2021; Kulkarni et al. 2026) and present in a fraction of the targets ranging from ~26% (Neha et al. 2025) to ~32% (Lu et al. 2024), depending on the selection criteria.

Besides these statistical studies on large samples of targets, it is also very instructive to focus on specific objects and use NEO-WISE data to constrain their physical properties and, possibly, understand their nature. A good target for this type of study is the protostar IRAS 20126+4104, a well-studied high-mass YSO driving a precessing jet and outflow surrounded by a Keplerian disk (see Cesaroni et al. 2025 and references therein). Its mass is estimated to be ~12 M_⊙ (Chen et al. 2016), consistent with a bolometric luminosity of ~10⁴ L_⊙ (Johnston et al. 2011; Cesaroni et al. 2023). The precession of the bipolar jet and outflow hints at a companion (Shepherd et al. 2000; Cesaroni et al. 2005; Caratti o Garatti et al. 2008), although much less massive than the primary (Cesaroni et al. 2023). Estimates of the distance to IRAS 20126+4104, measured with H₂O maser parallax, range from ${1.33}_{-0.15}^{+0.19}$ (Hirota et al. 2020) to 1.64±0.05 kpc (Moscadelli et al. 2011). For the present study, we adopted the value with the smallest error, i.e. 1.64 kpc.

The main features of the source are shown in Fig. 1. The IR emission traces a bipolar nebula oriented SE–NW, consistent with the bipolar outflow, emerging from a dusty core revealed by the millimetre continuum emission. The protostar is embedded in the core and lies at the centre of a Keplerian disk whose plane is almost perpendicular to the outflow axis.

A recent study by (Szymczak et al. 2024) reveals periodic variation in the 6 GHz Class II CH₃OH maser emission, showing significant correlation with the IR emission from the NEOWISE database. If variability is a clue to the physical mechanisms at work, periodic variability is even more so. Its association with a rotating circumstellar disk and precessing jet makes the picture more intriguing and prompted deeper analysis of the IR emission from IRAS 20126+4104 at different epochs. We describe the archival data used for our study in Sect. 2, show the results obtained in Sect. 3, and present our interpretation in Sect. 4. Finally, we draw our conclusions in Sect. 5.

Fig. 1 Overlay of the 3.6 µm Spitzer/IRAC image (contours) on the 2.2 µm LBT/SOUL image from Massi et al. (2023). The contour levels range from 30 to 430 in steps of 50 MJy/sterad. The dashed ellipse denotes the full width at half power of the 1.4 mm continuum emission mapped with ALMA by Cesaroni et al. (2025). The cross marks the expected position of the protostar.

2 Archival data analysis

The NEOWISE¹ (Mainzer et al. 2011; Mainzer et al. 2014) database was used to obtain multi-epoch images of the IRAS 20126+4104 region. In particular, we used single-exposure high-quality images (qual_frame = 5 or 10) in the W1 band (3.4 µm) containing the source position, spanning 2014–2024. We ignored the W2 band at 4.6 µm, because the emission is partially saturated. The images that were acquired within a few days have been averaged together with the ICORE tool² to improve the signal-to-noise ratio. In this way, we obtained 21 images at regular intervals of ~6 months. The corresponding dates are listed in columns 2 and 3 of Table 1, while the images are shown in Fig. 2. The same table also provides the emission peak positions from a 2D Gaussian fit, as well as the flux densities measured within the two polygons drawn in each panel of Fig. 2. These polygons encompass the emission from the NW blueshifted ($S_{\nu }^B$), and SE redshifted ($S_{\nu }^R$), lobes of the IRAS 20126+4104outflow.

To extend monitoring of the region over a longer time interval, we applied the same procedure to the ALLWISE data (Wright et al. 2010) at 3.4 µm acquired in 2010, obtaining two additional images. Moreover, we retrieved the Spitzer IRAC image at 3.6 µm (Werner et al. 2004; Fazio et al. 2004) obtained by (Qiu et al. 2008). The corresponding observation dates, positions, and flux densities are given in the first three rows of Table 1. We stress that the Spitzer data were smoothed to the same resolution as the ALLWISE and NEOWISE data before deriving the relevant parameters.

3 Results

As mentioned, Szymczak et al. (2024) were the first to detect a periodic variation in IRAS 20126+4104. In addition to discovering variability in the CH₃OH masers, they revealed a correlation between maser emission and mid-IR flux from the NEOWISE catalogue. The same authors found that the peaks of the different NEOWISE images are distributed along a straight line parallel to the outflow axis and that the peak position changes with time. Our results fully confirm this scenario, as illustrated in Fig. 3. The linear fit to the peaks has a position angle of ~−60°, consistent with that of the outflow (see e.g. Cesaroni et al. 2025); it is almost perpendicular to the projected major axis of the disk, whose position angle is ~49°, as shown in Fig. 3a. This same figure shows that the 1.4 mm continuum peak is offset from the outflow axis by ~0.″9 or ~1500 au. The (sub)millimetre continuum peak is a reasonable indication of the protostar position and such an offset suggests that it may be drifting from NE to SW, as suggested by Massi et al. (2023).

In Fig. 3b, we plot the peak position of the mid-IR emission as a function of time. Clearly, the peak oscillates between two extremes, and we fitted the pattern with a sinusoidal curve using the MFIT command from the GILDAS software³. The fitted function is given by s = ∆s sin[(2π/T) t + Φ] + s₀, where s is the offset and t is the time, while ∆s, T, Φ, and s₀ denote the free parameters of the fit. The best fit is obtained for a period T=6.7±0.2 yr, consistent within the errors with the 6.9±0.2 yr period derived by Szymczak et al. (2024) for the CH₃OH maser. In the following we assume T=6.8 yr, a reasonable compromise between the two estimates.

The explanation for the periodic pattern in Fig. 3b is shown in Fig. 2. The SE (redshifted) and NW (blueshifted) lobe intensities appear to vary but not simultaneously: when one brightens, the other dims. Consequently, the 2D Gaussian fit to the image shows the peak alternating between lobes along the straight line in Fig. 3a.

This anti-correlated variation is clearly shown in Fig. 4a – where we plot the two lobe fluxes as a function of time – and likewise emphasized in Fig. 4b, which shows their flux ratio. A striking result evident from this plot is that the Spitzer data (first epoch) intensity far exceeds that of the ALLWISE and NEOWISE data. We believe this is not due to the relative calibration between the two instruments, because the flux difference for the field stars between the ALLWISE and the Spitzer images is ~2%, compared to ~100% for the IRAS 20126+4104 outflow lobes. Such a ratio hints at an episodic outburst from the massive protostar, possibly due to an accretion event analogous to that observed, for example, in S255 NIRS3 (see Garatti et al. 2017).

After excluding the Spitzer data, we fitted two sinusoidal functions to the fluxes of the two lobes, as for the offsets, but this time the period was fixed to the assumed value of 6.8 yr. The fits yield a phase difference of 2.5±0.2 yr between the two light curves. Although a perfect anti-correlation would correspond to a phase difference of T/2=3.4 yr, we can conclude that the emissions from the two lobes are significantly out of phase.

It is worth noting that, despite the remarkable difference in intensity between the Spitzer and WISE data, the ratio of the Spitzer fluxes between the two lobes is consistent with the extrapolation of the fit to the WISE data only (see Fig. 4b). This suggests that the observed variability in the mid-IR emission from IRAS 20126+4104 is due to the overlap of two phenomena: a regular, periodic variation plus a random increase in intensity that may be related to accretion episodes.

An anti-correlated variation was also observed by Szymczak et al. (2024) among the different CH₃OH maser features. In particular, they divided the maser spots into two groups, one toward the SE and the other toward the NW (see their Fig. 8). They find the former to be correlated with the near-IR emission of the SE lobe observed by Massi et al. (2023); we confirm this result in Fig. 4a, where the CH₃OH maser intensity (dotted curve; scaled by an arbitrary amount to enable comparison with IR data) closely follows the best fit (red curve) to the 3.4 µm SE lobe flux.

In light of the NEOWISE results, we also re-analysed the near-IR data of Massi et al. (2023), who observed an anti-correlation between the 2.2 µm continuum emission from the two lobes and hypothesised regular variability with a period of 12 yr. In Fig. 5 we show the flux densities obtained by integrating the near-IR emission inside the polygon S1+S2 for the SE lobe and the polygons N1+N2+N3 for the NW lobe (see Fig. 2 of Massi et al. 2023 for a definition of these polygons). To extend the time interval as much as possible, we also added the flux densities obtained from the K-band image of Cesaroni et al. (1997) acquired on August 26, 1996 (MJD=50 321 days). As for NEOWISE, we fitted the data with two sinusoids, assuming a fixed period of 6.8 yr. In this case, the fits to the two lobe intensities and their ratio are satisfactory, confirming that the anti-correlated variation is also present in the near-IR.

To lend further support to this conclusion, we compare the flux ratio between the two lobes in the mid-IR with that in the near-IR in Fig. 6. The phase difference between the two sinusoidal fits is marginal (0.6±0.3 yr) and may be due to poor sampling of the near-IR light curve.

Overall, the results based on 28 years of monitoring confirm the existence of a mechanism that drives the IR emission from IRAS 20126+4104 to undergo regular variations with a period of ~6.8 yr. We note that Szymczak et al. (2024) report significant variability in the mid-IR emission from IRAS 20126+4104but find no clear periodicity in the 3.4 µm and 4.6 µm light curves shown in their Fig. 6. This is explained by their use of the data from the NEOWISE catalogue, which include the total emission from both lobes. Since the emissions from the lobes are anti-correlated, their sum makes it difficult to identify a periodic pattern. It is also worth noting that our analysis covers a longer time interval as it also includes the last three epochs (in 2023–2024) of NEOWISE monitoring as well as three previous measurements from ALLWISE and Spitzer.

Fig. 2 NEOWISE images of the 3.4 µm emission towards IRAS 20126+4104 at 21 epochs. The number in the bottom right corner of each panel indicates the observation date given in the first column of Table 1. The offsets are relative to position α(J2000)=20h14m26.s0364, δ(J2000)=41°13′32.″516. The dashed and dotted patterns outline the polygons used to derive the mid-IR flux densities of the SE (redshifted) lobe and NW (blueshifted) lobe of the bipolar outflow.

Fig. 3 Result of the 2D Gaussian fit to the mid-IR emission from IRAS 20126+4104. (a) Distribution of the emission peaks at different epochs (from Table 1). The solid red line with PA=–60° represents a linear fit to the data. The solid green circle marks the 1.4 mm continuum peak, a good proxy for the protostar position. The dashed green line indicates the direction of the plane of the circumstellar disk. (b) Offsets of the mid-IR peaks, measured along the red line in the top panel, as a function of time. The horizontal green line indicates the offset corresponding to the projection of the green circle onto the red line in the top panel. The red sinusoid represents the best fit to the data.

4 Discussion

Our findings imply that models aiming to explain observed variability from IRAS 20126+4104 must satisfy four constraints:

1. The emission from the outflow lobes must be periodical.

2. The emissions from the two lobes must be anti-correlated.

3. The difference between the peak and the dip of the light curves must be ~40% of the peak (see Fig. 4a). We define this as the ‘dimming factor’ η=0.4.

4. The light curves must have a sinusoidal shape, meaning that the widths of the peaks must be similar to the widths of the dips.

With this in mind, we discuss below four possible mechanisms that could explain the observed variability: pulsations, precession, rotation, and revolution.

Fig. 4 Variation in mid-IR emission from the outflow lobes. (a) Integrated flux density over the SE (red points) and NW (blue points) lobes as a function of time (from Table 1). The sinusoids represent the best fits to the solid points of corresponding colour. The empty points correspond to the Spitzer data, which have not been accounted for in the sinusoidal fits. The dotted sinusoid represents the fit (with arbitrary intensity) to the CH3OH maser intensity obtained from Fig. 1 of Szymczak et al. (2024). (b) NW-to-SE lobe flux ratio as a function of time. The green curve represents the ratio of the blue-to-red sinusoids in the top panel.

4.1 Pulsations

By ‘stellar pulsations’ we refer generically to any periodic variability in stellar luminosity. This can be due to instabilities in a star internal structure or regular accretion events from the circumstellar disk. Whatever the reason for the periodicity in the stellar luminosity, this phenomenon cannot explain the light curves in Fig. 4, because it is incompatible with the time shift between the emission peaks from the two lobes. In fact, the photons emitted by the star would reach the lobes of the outflow at the same time. In this case, the observed ~2.5-year delay (see Sect. 3) between the emission peaks of the blue lobe and the red lobe must be equal to the light-travel time between the lobes along the line of sight. The latter is given by

$$t_{light}=\frac{2L\tan \alpha }{c}=2.5\text{ yr.}$$(1)

where L is the size of a lobe projected onto the plane of the sky, α is the angle between the outflow axis and the plane of the sky, and c is the speed of light.

From Fig. 2 we derive L < 10″ = 0.08 pc and hence obtain α > 78° from Eq. (1). Such a value is by far greater than the estimates for the inclination angle of the outflow obtained from the SiO(5–4) line (Cesaroni et al. 2025), the H₂ line (Massi et al. 2023), and the 22 GHz H₂O masers (Moscadelli et al. 2005), which range from 3° to 8°. We thus believe that stellar pulsations are not a viable explanation for the observed variability.

Finally, it is also worth noting that the model-predicted pulsational instability phase in bloated protostars (Inayoshi et al. 2013) has been tested and found consistent with the periodic variability observed in the massive protostars IRAS 19520+2759 (Pandey et al. 2025) and G353.273+0.641 (Harajiri et al. 2026). However, this cannot be the case for IRAS 20126+4104. In fact, as noted by Szymczak et al. (2024), for a period of 6.8 yr the same model would imply a stellar luminosity of 10⁶ L_⊙, which is two orders of magnitude greater than that of IRAS 20126+4104 (see Eq. (1) of Inayoshi et al. 2013).

Fig. 5 Same as Fig. 4 for the near-IR data taken from Cesaroni et al. (1997) and Massi et al. (2023).

Fig. 6 Comparison between the flux density ratios between the NW and SE lobes at mid-IR (black circles) and near-IR (magenta squares) as a function of time. The solid black and dashed magenta curves represent the best fits to the points of the same colour. The data and the fits are the same as in Figs. 4b and 5b.

4.2 Precession

Massi et al. (2023) and Szymczak et al. (2024) propose a scenario where the observed periodic variability is due to the precession of the innermost part of the circumstellar disk. However, this mechanism cannot explain any variability if the disk is axially symmetric: in this case only the orientation of the precessing disk would change with time, whereas its inclination and solid angle shadowed by the disk would remain constant. The model works only if the disk is, for example, optically thick on one side and optically thin on the diametrically opposite side. In this case stellar photon shadowing alternates in the two lobes, thus explaining the anti-correlation in their variability.

To describe the precessing disk effect on the observed luminosity in a quantitative way, we assumed a simple model where only half of the disk is opaque in the mid-IR, the outflow lobes are conical with opening angle θ₀, and the disk precesses about the outflow axis at an angle ψ₀ between the disk rotation axis and the outflow axis. As detailed in Appendix A, one can express the dimming factor as a function of θ₀ and ψ₀ using Eq. (A.12). Therefore, from η(θ₀, ψ₀)=0.4 (see constraint no. 3), we derived ψ₀ as a function of θ₀, which implies that the disk inclination must be ψ₀ ≥ 72° for any lobe opening angle (see Fig. A.2).

This result has important implications on the disk size. During the precession, the disk plane wobbles between −ψ₀ and +ψ₀ in half a period; therefore, the velocity at the maximum radius, R₀, of the precessing disk is equal to

$$v=4\frac{R_0{\psi }_0}{T}.$$(2)

For the disk to remain coherent, such a velocity cannot exceed the speed of sound, i.e. v ≤ c_s ≃ 1 km s⁻¹. With the additional condition ψ₀ ≥ 72°, we obtain

$$R_0=\frac{vT}{4{\psi }_0}<0.28\text{ au.}$$(3)

Such a small radius is not acceptable for two reasons. First, it is smaller than the expected dust-destruction radius for a 12 M_⊙ star (see Table 3 of Vaidya et al. 2009), which would make the disk optically thin to mid-IR radiation. And, most importantly, the rotation period around a 12 M_⊙ star at this radius is ~16 days, which is by far shorter than the 6.8 yr precession period. This result contradicts the hypothesis that a sector of the disk remains opaque, since the material would be azimuthally rehashed on a much shorter timescale than the precession period. In conclusion, we believe that precession of the circumstellar disk cannot explain the observed periodicity.

Incidentally, we stress that this conclusion does not rule out the possibility that part of the disk precesses on a much longer timescale – as suggested by the precession period of the associated outflow (~ 2 × 10⁴ yr; Cesaroni et al. 2005). Furthermore, the hypothesis that the disk might not be axially symmetric (as suggested by Cesaroni et al. 2014) remains valid, although insufficient on its own to account for the observed variability.

4.3 Rotation

Similar to Sun spots, the protostellar surface might present dimmer regions that, corotating with the protostar, could produce the observed periodical variations. This situation is similar to that of weak-line T Tauri stars, where cold spots covering up to ~45% of the stellar surface account for the observed long-term variability, which persists for many years (Grankin at al. 2008). In our case, the problem with this interpretation is that young early-type stars are known to have rotation periods of several days – far shorter than the 6.8 yr period of the IRAS 20126+4104 light curves. For example, Abt et al. (2007) found that B-type stars typically rotate at 25% of the breakup velocity, which for a zero-age main-sequence star of 12 M_⊙ with a radius of ~5 R_⊙ (Panagia et al. 1973) corresponds to a rotation period of ~1.5 days.

However, accreting protostars are expected to be bloated (Hosokawa et al. 2009, 2010) with radii of up to a few 100 R_⊙. In this case rotation should be slowed, as angular momentum is conserved. Since the momentum of inertia of a sphere is proportional to R², we conclude that the radius of the bloated protostar should be equal to $5R_{\odot }\sqrt{6.8\text{ yr}/1.5\text{ days}}\simeq 200R_{\odot }$, a value comparable to those theoretically predicted for massive protostars undergoing disk accretion (see Hosokawa et al. 2010). These models also predict a convective zone close to the surface of bloated protostars, consistent with the occurrence of spots in the stellar photosphere.

Based on these considerations, we propose a model in which the protostar features a spot covering half of the equator and extending by ±∆θ across it. The star rotation axis must be inclined by at least ∆θ with respect to the disk axis, so that during rotation the spot gradually moves from one half-space relative to the disk plane to the opposite one. This makes the emission dimmer for half a rotation period – alternatively in the two lobes – thereby explaining the anti-correlated variability of the light curves.

In this model the dimming factor η is equal to the ratio between the solid angle subtended by the spot with respect to the centre of the protostar, Ω_S = 2π sin ∆θ, and the solid angle of a hemisphere, 2π. Therefore, η = sin ∆θ = 0.4, which yields ∆θ ≃ 24°.

As demonstrated in Appendix B, we can express the flux densities of the two lobes using Eq. (B.5):

$$S_{\nu }^B(\varphi )=S_{\nu }^{\text{B,max}}\left(1-\frac{\pi -|\varphi (t;\Phi )|}{\pi }\eta \right),$$(4)

$$S_{\nu }^R(\varphi )=S_{\nu }^{R,max}\left(1-\frac{\pi -|\varphi (t;\Phi +\Delta \Phi )|}{\pi }\eta \right),$$(5)

where $S_{\nu }^{\text{B},\max }$ and $S_{\nu }^{\text{R},\max }$ are assumed to remain constant in time, Φ is the phase of the spot (see Eq. (B.6)), and a phase difference, ∆Φ, is allowed between the two lobes.

Using Eqs. (4) and (5), we simultaneously fitted the fluxes of the two outflow lobes. The free parameters of the fit are $S_{\nu }^{\text{B}}$, $S_{\nu }^{\text{R}}$, Φ, and ∆Φ. Moreover, we fixed the period to T=6.8 yr. The best fit to the light curves, obtained for $S_{\nu }^{\text{B}}=\text{0}.225\pm 0.006\text{\hspace{0.17em}Jy}$, $S_{\nu }^{\text{R}}=\text{0}.229\pm 0.005\text{\hspace{0.17em}Jy}$, and ∆Φ=4.0±0.3 yr is shown in Fig. 7, where we also show the corresponding ratio. Although the fitted curve is clearly very schematic due to the approximations adopted, its shape is well reproduced.

It is also worth noting that, in this model, any random variation in stellar luminosity (e.g. due to accretion outbursts) appears in $S_{\nu }^{\text{B}}$ and $S_{\nu }^{\text{R}}$ as a multiplicative factor that cancels out in the ratio between the two lobes’ fluxes. Therefore, the curve in Fig. 7b represents only the periodic variation of the emission. This interpretation is supported by the Spitzer data (empty points), which are consistent with the model in Fig. 7b, despite the large discrepancy seen in Fig. 7a.

With this in mind, we can elaborate more on Eqs. (4) and (5), by introducing an additional dependence on time:

$$S_{\nu }^{\text{B}}(t)=S_{\nu }^{*}(t)f_{\text{B}}\left(1-\frac{\pi -|\varphi (t;\Phi )|}{\pi }\eta \right),$$(6)

$$S_{\nu }^{\text{R}}(t)=S_{\nu }^{*}(t)f_{\text{R}}\left(1-\frac{\pi -|\varphi (t;\Phi +\Delta \Phi )|}{\pi }\eta \right),$$(7)

where $S_{\nu }^{*}(t)$ denotes the flux density emitted by the star and f_B and f_R denote the fractions of that flux directed into the two lobes. In our previous fit, we assumed $S_{\nu }^{*}$ to be constant, in order to derive a mean shape of the periodic light curve. However, the stellar luminosity is likely to vary more or less randomly in time – unlike f_B and f_R, which depend on lobe structure and vary on much longer timescales than T. Thus, by taking the ratio of Eqs. (6)–(7), we obtain

$$S_{\nu }^{{\text{B}}^{\prime }}=\frac{f_{\text{B}}}{f_{\text{R}}}{S}_{\nu }^{{\text{R}}^{\prime }},$$(8)

where we define

$$S_{\nu }^{{\text{B}}^{\prime }}=\frac{S_{\nu }^{\text{B}}}{1-\frac{\pi -|\varphi (t;\Phi )|}{\pi }\eta },$$(9)

$$S_{\nu }^{{\text{R}}^{\prime }}=\frac{S_{\nu }^{\text{R}}}{1-\frac{\pi -|\varphi (t;\Phi +\Delta \Phi )|}{\pi }\eta }.$$(10)

In practice, these correspond to the flux densities of the lobes after correcting for the dimming caused by the star spot. Equation (8) indicates that the corrected fluxes must be correlated. This is exactly what one would expect if there were no spot on the stellar surface, because both lobes would be illuminated at the same time by the star.

To test this prediction, we compare in Fig. 8 the flux densities of the red and blue lobes before (black points) and after (red points) the corrections. Although they show a large spread, it is clear that the observed fluxes are anti-correlated, consistent with the analysis of the light curves in Sect. 3. Conversely, after applying the correction $S_{\nu }^{\text{B}\prime }$ and $S_{\nu }^{\text{B}\prime }$ are roughly proportional to each other, as in Eq. (8). We note that the Spitzer fluxes are not included in this plot because of their large deviation from the mean; however, they would obviously reinforce the correlation between $S_{\nu }^{\text{B}\prime }$ and $S_{\nu }^{\text{B}\prime }$.

Fig. 7 Same as Fig. 4, with the solid blue, red, and green lines representing the fits obtained from Eqs. (4) and (5).

Fig. 8 Flux densities of the red vs blue lobe measured at 23 epochs with ALLWISE and NEOWISE. The Spitzer measurement is not shown. The black symbols indicate the observed fluxes, $S_{\nu }^{\text{B}}$ and $S_{\nu }^{\text{R}}$. The red symbols denote the fluxes $S_{\nu }^{\text{B}\prime }$ and $S_{\nu }^{\text{R}\prime }$, corrected for periodic oscillations using Eqs. (9) and (10). The solid lines are linear fits to the points of the same colour, and the dashed lines indicate the corresponding uncertainties.

4.4 Revolution

The very same model fit presented in Sect. 4.3 can be obtained with a body orbiting the star and intercepting the stellar photons emitted into the solid angle given by Eq. (B.3). Such a body could be a dense, dusty filament distributed along an orbit of radius ~8 au, which corresponds to a revolution period of 6.8 yr around a 12 M_⊙ star. This orbit must be inclined by at least ∆θ = 24° with respect to the disk plane, as in the model in Sect. 4.3. We stress that any compact body with a smaller ∆ϕ would produce too steep a flux variation, which could not fit the smooth intensity changes observed between the peaks and the dips of the light curves (see constraint no. 4).

Although the existence of such a filament cannot be ruled out a priori, it is difficult to justify the presence of a significant amount of dust at only 8 au from a 12 M_⊙ star. In fact, according to Vaidya et al. (2009), the dust sublimation radius due to circumstellar disk heating around a 10 M_⊙ star is at least 10 au. We conclude that periodic shadowing by a body orbiting around the protostar is not a plausible explanation for the observed flux variations.

5 Summary and conclusions

We presented robust evidence for periodic variability in mid-IR emission from the massive protostar IRAS 20126+4104. The flux densities of the two outflow lobes obtained from the Spitzer, ALLWISE, and NEOWISE databases appear to vary over a couple of decades following a sinusoidal light curve with period of ~6.8 yr. In particular, the emissions from the two lobes are anti-correlated, with a phase difference of 2.5±0.2 yr. These results have also been confirmed by the the near-IR data from Massi et al. (2023) and Cesaroni et al. (1997).

We propose four possible interpretations of the results, based on periodic variability in stellar luminosity, circumstellar disk precession, protostellar rotation with a large spot, and stellar photon shadowing by a body orbiting at 8 au from the protostar. The only acceptable model that satisfies all the characteristics of the observed light curves is that involving a dark spot covering 20% of the entire stellar surface. The 6.8 year-long rotation period can be explained with angular momentum conservation if the proto-star has expanded up to a radius of ~200 R_⊙, as theoretically predicted. Therefore, our study supports the existence of bloated massive protostars and suggests additional methods for characterizing these sources. These include detecting slow rotation and periodic emission modulation due to star spots, presumably associated with a phase of active surface convection.

We conclude that in all likelihood IRAS 20126+4104 is a slowly rotating, bloated massive protostar. Its low temperature and faint Lyman continuum are consistent with the absence of a detectable HII region around IRAS 20126+4104. Although the existence of a large spot on the surface of a such a star may seem unusual, it remains the most likely explanation for the observed periodic anti-correlated variability. As Sherlock Holmes’ famously said, ’[When] you have eliminated the impossible, whatever remains, however improbable, must be the truth’ (Doyle 1890, The Sign of the Four).

Acknowledgements

I am indebted to Daniele Galli for insightful discussions and valuable suggestions which improved the quality of this article, as well as to Rino Bandiera for creating Fig. A.1 thanks to his unrivalled expertise in the Mathematica software. I also thank Fabrizio Massi for kindly providing me with the near-IR flux densities used in this study and Eugenio Schisano for useful suggestions on the Spitzer and WISE data. This publication makes use of data products from the Wide-field Infrared Survey Explorer, which is a joint project of the University of California, Los Angeles, and the Jet Propulsion Laboratory/California Institute of Technology, and NEOWISE, which is a project of the Jet Propulsion Laboratory/California Institute of Technology. WISE and NEO-WISE are funded by the National Aeronautics and Space Administration. This work is also based in part on archival data obtained with the Spitzer Space Telescope, which was operated by the Jet Propulsion Laboratory, California Institute of Technology under a contract with NASA.

References

Appendix A Dimming factor due to precessing disk

We want to estimate the effect of shadowing of a circumstellar disk on the stellar photons emitted inside an outflow lobe. In our assumptions the lobe is a cone with opening angle 0 ≤ θ₀ ≤ π/2 and vertex coincident with the star, and the disk axis forms an angle 0 ≤ ψ₀ ≤ π/2 with the outflow axis (see Fig. A.1). Since the disk can intercept part of the photons from the star only if it intersects the lobe of the outflow, we impose the condition ψ₀ > π/2 − θ₀.

The dimming factor, namely the ratio between the flux intercepted by the disk and the total flux emitted by the star inside the conical lobe, is equal to the ratio between the solid angle between the disk plane and the cone, Ω, and the total solid angle of the cone, Ω_cone. In spherical coordinates, the latter is given by the expression

$${\Omega }_{cone}={\displaystyle {\int }_0^{2\pi }d}\varphi {\displaystyle {\int }_0^{{\theta }_0}\sin }\theta d\theta =2\pi (1-\cos {\theta }_0)$$(A.1)

while to calculate the former we need to find the intersection between the cone and the plane. This is obtained from the equations

$$z=\tan {\psi }_{0}x$$(A.2)

$$z=\cot {\theta }_0\sqrt{x^2+y^2}$$(A.3)

where we have assumed that the plane of the disk contains the y axis. We remark that we consider only one of the lobes, for z > 0, and the intersections lie in the hemispace x > 0. The solution of these equations gives

$$y=\pm \frac{\sqrt{{\tan }^2{\psi }_0-{\cot }^2{\theta }_0}}{\cot {\theta }_0}x$$(A.4)

which are the projections of the intersections onto the x, y plane (red lines in Fig. A.1). The corresponding azimuthal angles are

$${\varphi }_{\pm }=\pm \arctan \left(\frac{\sqrt{{\tan }^2{\psi }_0-{\cot }^2{\theta }_0}}{\cot {\theta }_0}\right).$$(A.5)

We also need to express the angle ψ between the plane x, y and a generic line lying on the plane of the disk and passing through the origin, as a function of the azimuthal angle ϕ. For a generic point (x, y, z), the projection of its distance from the origin onto the x, y plane is given by the expression

$$\sqrt{x^2+y^2}=\cos \psi \sqrt{x^2+y^2+z^2}.$$(A.6)

Substituting y = x tan ϕ and Eq. (A.2) into Eq. (A.6) we obtain

$$\cos \psi (\varphi )=\sqrt{\frac{1+{\tan }^2\varphi }{1+{\tan }^2\varphi +{\tan }^2{\psi }_0}}$$(A.7)

and therefore

$$\sin \psi (\varphi )=\frac{\tan {\psi }_0}{\sqrt{1+{\tan }^2\varphi +{\tan }^2{\psi }_0}}$$(A.8)

Now, we can calculate the solid angle between the plane and the cone:

$$\text{Ω}={\displaystyle {\int }_{{\varphi }_{-}}^{{\varphi }_{+}}d\varphi }{\displaystyle {\int }_{\frac{\pi }{2}-\psi (\varphi )}^{{\theta }_0}\sin \theta \text{\hspace{0.17em}d}\theta }$$

$$\begin{array}{l}\Omega ={{\displaystyle \int }}_{{\varphi }_{-}}^{{\varphi }_{+}}d\varphi {\displaystyle {\int }_{\pi /2-\psi (\varphi )}^{{\theta }_0}\sin }\theta d\theta \\ =2\tan {\psi }_0{{\displaystyle \int }}_0^{{\phi }_{+}}\frac{d\phi }{\sqrt{1+{\tan }^2\phi +{\tan }^2{\psi }_0}}-2\cos {\theta }_0\arctan \left(\frac{\sqrt{{\tan }^2{\psi }_0-{\cot }^2{\theta }_0}}{\cot {\theta }_0}\right)\\ =2\tan {\psi }_0\frac{\sqrt{{\tan }^2{\psi }_0\cos \left(2{\phi }_{+}\right)+{\tan }^2{\psi }_0+2}}{\sqrt{2}\tan {\psi }_0\cos {\phi }_{+}\sqrt{{\tan }^2{\psi }_0+{\tan }^2{\phi }_{+}+1}}\times \arctan \left(\frac{\tan {\psi }_0\sin {\phi }_{+}}{\sqrt{1+{\tan }^2{\psi }_0{\cos }^2{\phi }_{+}}}\right)-2\cos {\theta }_0\arctan \left(\frac{\sqrt{{\tan }^2{\psi }_0-{\cot }^2{\theta }_0}}{\cot {\theta }_0}\right).\end{array}$$(A.9)

Fig. A.1 Sketch of the model consisting of a conical lobe intersected by the plane of the disk. The axis of the cone coincides with the z axis and the plane contains the y axis. The black lines are the intersections between the plane and the cone, while the red lines are their projections onto the x, y plane. Also indicated are the angle ψ0 between the plane x, y and the plane of the disk, the opening angle θ0 of the cone, and the angles ϕ± between the x axis and the red lines.

Here we have used the expression of sin ψ in Eq. (A.8) and the property sin ψ(ϕ) = sin ψ(−ϕ).

From Eq. (A.5) one obtains

$$\cos {\varphi }_{+}=\frac{1}{\sqrt{1+{\tan }^2{\varphi }_{+}}}=\frac{\cot {\theta }_0}{\tan {\psi }_0}$$(A.10)

$$\sin {\varphi }_{+}=\sqrt{1-{\cos }^2{\varphi }_{+}}=\frac{\sqrt{{\tan }^2{\psi }_0-{\cot }^2{\theta }_0}}{\tan {\psi }_0}$$(A.11)

By substituting these and Eq. (A.1) into Eq. (A.9), after some algebra one derives the expression for the dimming factor:

$$\eta =\frac{\Omega }{{\Omega }_{cone}}=\frac{1}{\pi (1-\cos {\theta }_0)}\left[\arctan \left(\sqrt{\frac{{\tan }^2{\psi }_0-{\cot }^2{\theta }_0}{1+{\cot }^2{\theta }_0}}\right)-\cos {\theta }_0\arctan \left(\frac{\sqrt{{\tan }^2{\psi }_0-{\cot }^2{\theta }_0}}{\cot {\theta }_0}\right)\right].$$(A.12)

We can now relate θ₀ to ψ₀ for a given value of η. In Fig. A.2 we plot ψ₀ as a function of θ₀ for η=0.4, the value obtained for IRAS 20126+4104 (see Sect. 4). We conclude that in our case ψ₀ > 72°.

Fig. A.2 Inclination angle of the precessing disk as a function of the opening angle of the outflow lobes, for a dimming factor η = 0.4.

Appendix B Light curve of rotating star with spot

We propose a simple model to describe the emission inside one of the outflow cavities due to a rotating star with a dark spot on its surface. We assume that the spot extends along half of the equator of the star and by ±∆θ in latitude. The rotation axis of the star must be inclined by at least ∆θ with respect to the disk axis, so that in one rotation period the spot moves from one side with respect to the disk plane, to the other side.

The flux emitted in one of the lobes can be expressed as

$$S_{\nu }(\varphi )=S_{\nu }^{max}\left(1-\frac{\Omega (\varphi )}{2\pi }\eta \right)$$(B.1)

where $S_{\nu }^{\max }$ is the maximum flux density emitted when the spot lies on the opposite side of the disk plane with respect to the lobe, Ω is the solid angle subtended by the portion of the spot that lies on the same side as the lobe, and ϕ is the azimuthal angle that we assume to lie in the range −π < ϕ < +π, with ϕ = 0 when the spot lies entirely on the same side as the lobe (i.e. when the minimum intensity is reached). The solid angle can be computed as

$$\Omega =\Delta \varphi {\int }_{\frac{\pi }{2}-\Delta }^{\frac{\pi }{2}+\Delta }sin\theta d\theta =2\Delta \theta sin\Delta \theta .$$(B.2)

Here, ∆ϕ is a function of ϕ as one can see from the sketch in Fig. B.1, where the ∆ϕ is the angle subtended by the thick black solid arc. In this fugure, we show four template cases with different orientations of the spot, for ϕ ≥ 0. It is easy to see that in all cases ∆ϕ = π−ϕ. Analogously, for ϕ < 0 one obtains ∆ϕ = π+ϕ, so that in general we can write

$$\Omega (\varphi )=2(\pi -|\varphi |)\sin \Delta \theta .$$(B.3)

For our purposes it is useful to express ∆θ as a function of the measurable quantity η, the dimming factor that by definition is equal to (see constraint n. 3 in Sect. 4)

$$\eta =\frac{S_{\nu }^{max}-S_{\nu }^{min}}{S_{\nu }^{max}}=\frac{S_{\nu }(\pi )-S_{\nu }(0)}{S_{\nu }(\pi )}=\frac{\Omega (0)}{2\pi }=\sin \Delta \theta$$(B.4)

where $S_{\nu }^{\min }$ is the minimum observed value of the flux. By substituting Eqs. (B.4) and (B.3) into Eq. (B.1), the expression of the flux takes the form

$$S_{\nu }(\varphi )=S_{\nu }^{max}\left(1-\frac{\pi -\left|\varphi \right|}{\pi }\eta \right).$$(B.5)

Fig. B.1 Sketch of the rotating star with an equatorial spot. The cyan circle is the star, while the blue and red colours indicate, respectively, the equator and the rotation axis of the star. The brown line denotes the plane of the disk seen edge on. Dashed patterns indicate a position under the disk with respect to the lobe under consideration. The thick black pattern marks the extension of the spot along the equator, with the solid and dashed arcs representing, respectively, the part of the spot lying inside the lobe of interest and that on the opposite side. The pictures inside the boxes show four template positions of the spot for the azimuthal angles given in the top right corner of the boxes. The point of view of the four panels is along the rotation axis, as indicated by the arrow in the figure to the right, which shows the star as seen by an observer whose line of sight lies in the plane of the disk and is perpendicular to the star rotation axis.

To fit the data we must express S_ν as a function of time. We thus define the angle

$${\varphi }^{\prime }(t;\Phi )=2\pi \frac{t}{T}+\Phi$$(B.6)

where T = 6.8 yr is the rotation period and Φ is the phase. Finally, ϕ(t; Φ) is obtained by reducing ϕ′(t; Φ) into the range between −π and +π.

All Tables

All Figures

	Fig. 1 Overlay of the 3.6 µm Spitzer/IRAC image (contours) on the 2.2 µm LBT/SOUL image from Massi et al. (2023). The contour levels range from 30 to 430 in steps of 50 MJy/sterad. The dashed ellipse denotes the full width at half power of the 1.4 mm continuum emission mapped with ALMA by Cesaroni et al. (2025). The cross marks the expected position of the protostar.
In the text

Fig. 2 NEOWISE images of the 3.4 µm emission towards IRAS 20126+4104 at 21 epochs. The number in the bottom right corner of each panel indicates the observation date given in the first column of Table 1. The offsets are relative to position α(J2000)=20h14m26.s0364, δ(J2000)=41°13′32.″516. The dashed and dotted patterns outline the polygons used to derive the mid-IR flux densities of the SE (redshifted) lobe and NW (blueshifted) lobe of the bipolar outflow.

In the text

Fig. 3 Result of the 2D Gaussian fit to the mid-IR emission from IRAS 20126+4104. (a) Distribution of the emission peaks at different epochs (from Table 1). The solid red line with PA=–60° represents a linear fit to the data. The solid green circle marks the 1.4 mm continuum peak, a good proxy for the protostar position. The dashed green line indicates the direction of the plane of the circumstellar disk. (b) Offsets of the mid-IR peaks, measured along the red line in the top panel, as a function of time. The horizontal green line indicates the offset corresponding to the projection of the green circle onto the red line in the top panel. The red sinusoid represents the best fit to the data.

In the text

Fig. 4 Variation in mid-IR emission from the outflow lobes. (a) Integrated flux density over the SE (red points) and NW (blue points) lobes as a function of time (from Table 1). The sinusoids represent the best fits to the solid points of corresponding colour. The empty points correspond to the Spitzer data, which have not been accounted for in the sinusoidal fits. The dotted sinusoid represents the fit (with arbitrary intensity) to the CH3OH maser intensity obtained from Fig. 1 of Szymczak et al. (2024). (b) NW-to-SE lobe flux ratio as a function of time. The green curve represents the ratio of the blue-to-red sinusoids in the top panel.

In the text

	Fig. 5 Same as Fig. 4 for the near-IR data taken from Cesaroni et al. (1997) and Massi et al. (2023).
In the text

	Fig. 6 Comparison between the flux density ratios between the NW and SE lobes at mid-IR (black circles) and near-IR (magenta squares) as a function of time. The solid black and dashed magenta curves represent the best fits to the points of the same colour. The data and the fits are the same as in Figs. 4b and 5b.
In the text

	Fig. 7 Same as Fig. 4, with the solid blue, red, and green lines representing the fits obtained from Eqs. (4) and (5).
In the text

Fig. 8 Flux densities of the red vs blue lobe measured at 23 epochs with ALLWISE and NEOWISE. The Spitzer measurement is not shown. The black symbols indicate the observed fluxes, $S_{\nu }^{\text{B}}$ and $S_{\nu }^{\text{R}}$. The red symbols denote the fluxes $S_{\nu }^{\text{B}\prime }$ and $S_{\nu }^{\text{R}\prime }$, corrected for periodic oscillations using Eqs. (9) and (10). The solid lines are linear fits to the points of the same colour, and the dashed lines indicate the corresponding uncertainties.

In the text

Fig. A.1 Sketch of the model consisting of a conical lobe intersected by the plane of the disk. The axis of the cone coincides with the z axis and the plane contains the y axis. The black lines are the intersections between the plane and the cone, while the red lines are their projections onto the x, y plane. Also indicated are the angle ψ0 between the plane x, y and the plane of the disk, the opening angle θ0 of the cone, and the angles ϕ± between the x axis and the red lines.

In the text

	Fig. A.2 Inclination angle of the precessing disk as a function of the opening angle of the outflow lobes, for a dimming factor η = 0.4.
In the text

Fig. B.1 Sketch of the rotating star with an equatorial spot. The cyan circle is the star, while the blue and red colours indicate, respectively, the equator and the rotation axis of the star. The brown line denotes the plane of the disk seen edge on. Dashed patterns indicate a position under the disk with respect to the lobe under consideration. The thick black pattern marks the extension of the spot along the equator, with the solid and dashed arcs representing, respectively, the part of the spot lying inside the lobe of interest and that on the opposite side. The pictures inside the boxes show four template positions of the spot for the azimuthal angles given in the top right corner of the boxes. The point of view of the four panels is along the rotation axis, as indicated by the arrow in the figure to the right, which shows the star as seen by an observer whose line of sight lies in the plane of the disk and is perpendicular to the star rotation axis.

In the text