A Study of the Photometric and Spectroscopic Variations of the Prototypical FU Orionis–type Star V1057 Cyg

We performed the majority of our photometric observations in B, V, R_C , I_C , g', r', and i' filters at the Piszkéstető Mountain Station of Konkoly Observatory (Hungary) between 2005 and 2021. Three telescopes with three slightly different optical systems were used. In 2005–2007 we observed the star with the 1 m Ritchey–Chrétien-coudé (RCC) telescope, equipped with a 1300 × 1340 pixel Roper Scientific VersArray: 1300B CCD camera (pixel scale: 0306). The 60/90/180 cm Schmidt telescope, equipped with a 4096 × 4096 pixel Apogee Alta U16 CCD camera (pixel scale: 1027), was used in 2011–2019. In each of the BVR_C I_C filters, typically three images per night were taken. Since 2020 we started to use the Astro Systeme Austria AZ800 alt-azimuth direct drive 80 cm Ritchey–Chretien (RC80) telescope operating in fully autonomous mode. The optical setup with the effective focal length of F = 5600 mm yielded a pixel scale of 055 and a field of view of 188 × 188 for a 2048 × 2048 pixel FLI PL230 CCD camera. We obtained three images per night in BVg'r'i' filters.

The frames were calibrated for bias, dark, and flat field in the standard fashion. Photometry of V1057 Cyg and 12 comparison stars in its 8' vicinity was extracted using an aperture radius of 41 and sky annulus between 103 and 154 for RCC and Schmidt frames, and 55 and sky annulus between 11'' and 22'' for RC80 telescope frames. In order to eliminate system-related effects, photometric calibration was performed by fitting a color term using the magnitudes and colors of the comparison stars from the APASS DR9 catalog (Henden et al. 2016), after converting them from the Sloan to the Bessel system using transformations from Jordi et al. (2006). We note that many Schmidt observations actually targeted another, fainter young eruptive star, HBC 722 (Kóspál et al. 2011, 2016), and V1057 Cyg just happened to be in the field of view. As a consequence, V1057 Cyg saturated the detector in some of the R_C and I_C images, which were discarded from further analyses.

In addition to our national facilities, we occasionally used other telescopes. On 2006 July 20 and 2012 October 13 we obtained B, V, R_J , and I_J images of V1057 Cyg with the IAC80 telescope of the Instituto de Astrofísica de Canarias located at Teide Observatory (Canary Islands, Spain). It was equipped with the Tromsoe CCD Photometer (TCP) with a 92 × 90 field of view and a 0537 pixel scale. After the standard reduction steps for bias, dark, and flat-field correction, aperture photometry was done by using the same aperture and sky annulus size as for the Schmidt and RCC data. Photometric calibration was done using the same comparison stars, except for the two that fell outside the smaller field of view of the telescope. During 2019 August–September, in parallel with TESS, we additionally observed V1057 Cyg at the Northern Skies Observatory (NSO). We used the 0.4 m telescope equipped with BVI filters, operated remotely through Skynet. The calibration procedures and comparison stars were the same as above, but only the VI NSO filter data were of analysis quality.

We also observed V1057 Cyg with the 2.56 m Nordic Optical Telescope (NOT) at the Roque de los Muchachos Observatory, La Palma, in the Canary Islands (Plan ID 61–414, PI: Zs. M. Szabó). For optical imaging we used the Alhambra Faint Object Spectrograph and Camera (ALFOSC) on 2020 August 17. ALFOSC is a 2048 × 2064 pixel CCD231-42-g-F61 CCD camera with a field of view of 64 × 64 and pixel scale of 021. The Bessel BVR filter set was supplemented by an i interference filter, which is similar to the SLOAN i', but with a slightly longer effective wavelength of λ_eff = 0.789 μm. We obtained three images in each filter, with exposure times between 1.5 and 30 s. After the standard CCD reduction steps, we obtained aperture photometry using an aperture radius of 32 and a sky annulus between 64 and 86. Because of the small field of view, the magnitudes of V1057 Cyg were obtained based on only one comparison star.

Our photometric results are shown in Figures 1 and 2 and listed in Table 4 in Appendix A. The typical uncertainty of our measurements is 0.03 mag in B and 0.01 mag in all other filters.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

During 2019 August 15–October 7, V1057 Cyg was observed with 30-minute cadence with Camera 1 of the Transiting Exoplanet Survey Satellite (TESS; Ricker et al. 2015). The total coverage time of Sectors 15 and 16 of the satellite is 50.5625 days, but the run was interrupted three times, each for about 3.1–3.4 days to download the data to the MAST archive. ¹⁰ The calibrated full-frame images were processed in two main steps using the FITSH package (Pál 2012). First, the plate solution was derived based on the Gaia DR2 catalog—details of this complex procedure are described by Pál et al. (2020). As part of this step, we derived the flux zero-point with respect to the G_RP magnitudes of the matched Gaia sources, utilizing the similarities between the TESS and Gaia G_RP filter throughputs. By examining various TESS fields observed in the first two sectors, we found that the rms of our zero-level calibration is ∼0.015 mag. The photometry of the source was performed via differential image analysis using FITSH/ficonv and fiphot (Pál 2012). It requires a reference frame, which we constructed as a median of 11 individual 64 × 64 subframes obtained close to the middle of the observing sequence. As reference fluxes, required to correct for various instrumental and intrinsic differences between the target and the reference frames, we used the Gaia DR2 magnitudes. Data points affected by momentum wheel desaturation or significant stray light were flagged and removed, which caused three additional breaks of 1.2–1.3 days in the time coverage. The resulting typical formal uncertainties of the data are about 0.65 mmag. The TESS light curve of V1057 Cyg is presented in Figure 3.

Zoom In Zoom Out Reset image size

We complemented our work with data from the Wide-field Infrared Survey Explorer (WISE; Wright et al. 2010). We used data obtained in the 3.4 μm (W1) and 4.6 μm (W2) bands from 2010 up until the most recent data release in 2021 (Cutri 2012; Cutri et al.2014). Since V1057 Cyg was saturated, we corrected the data points using the saturation bias correction curves for the appropriate survey phase available in the WISE Explanatory Supplement. ¹¹ The corrected WISE data are shown in Figure 1 and listed in Table 5 in Appendix A.

We obtained NIR images in the J, H, and K_s bands at six epochs between 2006 July 15 and 2012 October 13 using the 1.52 m Telescopio Carlos Sanchez (TCS) at the Teide Observatory. This telescope is equipped with CAIN III, a 256 × 256 Nicmos 3 detector, which provided a pixel scale of 1'' in the wide optics configuration. Observations were performed in a five-point dither pattern in order to enable proper sky subtraction. The total integration time was typically 1 minute per dither position in each filter, split into 1.5–5 s exposures. The images were reduced using caindr, an IRAF-based data reduction package written by J. A. Acosta-Pulido and R. Barrena, ¹² as well as our own IDL routines. Data reduction steps included sky subtraction, flat-fielding, registration, and co-adding exposures by dither position and filter. To calibrate our photometry, we used the Two Micron All Sky Survey (2MASS) catalog (Cutri et al. 2003). The instrumental magnitudes of the target and all good-quality 2MASS stars in the field were extracted using an aperture radius of 2'' in all filters. We determined a constant offset between the instrumental and the 2MASS magnitudes by averaging typically 20–30 stars by means of biweight_mean—an outlier-resistant averaging method.

We also used the NOTCam instrument on the NOT on 2020 August 29. The instrument includes a 1024 × 1024 pixel HgCdTe Rockwell Science Center "HAWAII" array, and for wide-field (WF) imaging it has a 4' × 4' field of view (pixel scale: 0234). We obtained nine images in each of the JHK_s bands with 3.6 s exposures. Because of the brightness of our target in the infrared, we used a 5 mm diameter pupil mask intended for very bright objects to diminish the telescope aperture, which gave about 10% transmission. The images were reduced using the same method as described above at the TCS data reduction. The instrumental magnitudes of the target and the comparison star in the field were extracted using an aperture radius of 33 and a sky annulus between 66 and 94. The photometric calibration was performed in the same fashion as the TCS images. Typical photometric uncertainties are of 0.01–0.03 mag, and we present the results of the optical and infrared photometry in Appendix A, Table 4.

We obtained a new optical spectrum of V1057 Cyg with the high-resolution FIbre-fed Echelle Spectrograph (FIES) instrument on the NOT on 2020 August 17. We used a fiber with a larger entrance aperture of 25, which provided a spectral resolution R = 25,000, covering the 370–900 nm wavelength range. We obtained two spectra, each with 1800 s exposure time. During our analysis, we used the spectra reduced by the FIEStool software.

V1057 Cyg was also observed with the Bohyunsan Optical Echelle Spectrograph (BOES; Kim et al. 2002) installed on the 1.8 m telescope at the Bohyunsan Optical Astronomy Observatory (BOAO). It provides R = 30,000 in the wavelength range ∼400–900 nm. The first spectrum was obtained on 2012 September 11 and the last on 2018 December 18. We reduced these spectra in a standard way within IRAF: after standard calibrations on bias and flat field, the ThAr lamp spectrum was used for wavelength calibration, and continuum fitting was performed by the continuum task. Finally, heliocentric velocity correction was applied by the rvcorrect task and the published radial velocity (RV) of V1057 Cyg (−16 km s⁻¹; Herbig et al. 2003).

As no telluric standard stars were observed for either FIES or BOES, we performed the telluric correction using the molecfit software (Kausch et al. 2015; Smette et al. 2015) by fitting the telluric absorption bands of O₂ and H₂O. This generally provided good correction except for the deepest lines, where the detected signal was close to zero.

We present the spectroscopic observing log in Table 1.

On 2020 August 29, we used the NOTCam on the NOT to obtain new NIR spectra of V1057 Cyg and Iot Cyg (A5 V) as our telluric standard star in the JHK_s bands. We used the low-resolution camera mode (R = 2500) with ABBA dither positions, and exposure times ranged from 25 to 60 s (Table 1). For each image, flat-fielding, bad pixel removal, sky subtraction, aperture tracing, and wavelength calibration steps were performed within IRAF. For the wavelength calibration, the Xenon lamp spectrum was used. The hydrogen absorption lines in Iot Cyg were removed by Gaussian fitting. Then, the spectrum of V1057 Cyg was divided by the normalized spectrum of Iot Cyg for telluric correction. Finally, flux calibration was performed by applying the accretion disk model obtained using the NOT JHK_s photometry (Section 4.1).

On 2018 September 6, we observed V1057 Cyg with the Stratospheric Observatory for Infrared Astronomy (SOFIA; Young et al. 2012) using the Faint Object infraRed CAmera for the SOFIA Telescope (FORCAST; Herter et al. 2013). We obtained mid-infrared imaging in a series of short exposures in band F111 (10.6–11.6 μm) totaling ∼30 s, a single exposure in F056 (5.6 μm) for 37 s and F077 (7.5–8 μm) for 42 s, and R ∼ 100–200 spectra with G063 (5–8 μm) and G227 (17–27 μm) (Plan ID 06_062, PI: J. D. Green). The spectra were processed using the SOFIA pipeline and retrieved as Level 3 data products from the SOFIA Science Archive as ingested into the IRSA database. ¹³ The program was only partially observed in SOFIA Cycle 6, and thus the data do not cover the full 5–25 μm spectral range. The observations were performed in "C2N" (two-position chop with nod) mode, using the 47 slit, and the NMC (nod-match-chop) pattern in which the chops and nod are in the same direction and have the same amplitude. In each case, an off-source calibrator was selected, using the observation closest in zenith angle and altitude to the science target, as previously done with FU Orionis in SOFIA Cycle 4 (Green et al. 2016). We did not use dithering.

To study the long-term variability of V1057 Cyg, we complemented our work with data published in the literature (Mendoza 1971; Rieke et al. 1972; Landolt 1975, 1977; Simon 1975; Welin 1975, 1976; Simon et al. 1982; Simon & Joyce 1988; Kenyon & Hartmann 1991; Kopatskaya et al. 2013). Our V1057 Cyg monitoring began in 2005 and overlapped with that of Kopatskaya et al. (2013). This enabled us to determine systematic shifts between filters utilized in these two data sets. We found systematic differences between the two sets of photometry, which may be due to different apertures, filters, detector throughputs, and different comparison stars used. For plotting purposes, we shifted our B-band light curves by +0.12 mag, V-band lights curves by +0.08 mag, R_C -band light curves by +0.05 mag, and I_C -band light curves by −0.14 mag to be consistent with the earlier papers. In Table 4 in Appendix A we present our original photometry, i.e., without these offsets. The resulting long-term light curves covering the 1965–2021 time period are shown in Figure 1, while in Figure 2 we show in detail our Piszkéstető optical monitoring (starting from 2011), complemented with V- and g-band observations from the All-Sky Automated Survey for Supernovae (ASAS-SN; Shappee et al. 2014; Kochanek et al. 2017). In order to align the ASAS-SN V-band observations with our data, we applied a −0.026 mag shift to the former ones. For consistency with the 1971–2019 data set, we also transformed our Sloan r'i' data obtained in 2020 and 2021 into the Johnson–Cousins R_C I_C system using the transformation equations given by Jordi et al. (2006). A brief summary of the data used for the construction of the long-term photometric light curve is presented in Table 2. The table includes the dates of the observations, filters used, status of the source, and the relevant papers.

Both the archival and our new light curves firmly indicate that the post-outburst brightness evolution of V1057 Cyg is exceptional as compared to other FUOrs. Kolotilov (1990) noticed that after the phase of exponential decay, in 1984–1988, the brightness of V1057 Cyg has stabilized at a nearly constant level in all used filters. This was the so-called "first plateau" phase, which lasted until 1995. As mentioned in Section 1, UBV measurements taken in 1995–1996 revealed a sudden fading by about 1 mag in these bands, and this process (indicated by the vertical line in Figure 1) stopped in 1997 (Ibragimov 1997; Kolotilov & Kenyon 1997; Kopatskaya et al. 2002). Since 1997, the average brightness of V1057 Cyg has remained practically constant in all bands, and this phase is known as the "second plateau" (Kopatskaya et al. 2013). This plateau is also still present in the infrared region, as inferred from comparison of our JHK_s observations with the latest data points found in the literature (Kopatskaya et al. 2013).

The TESS light curve is presented in Figure 3, and we shifted Sector 15 and Sector 16 to match our light curve in the I_C band. We performed interpolation to shift Sector 15 by +0.07 mag and Sector 16 by −0.05 mag. The cause of the six major breaks in the data acquisition was described in Section 2. This precise light curve clearly shows the brightness changes occurring on a daily timescale, whose detailed investigation remains beyond the capabilities of the ground-based telescopes.

3.2.1. Long-term Variability as Seen from the Ground

As mentioned above, Kopatskaya et al. (2013) discovered wavelength-dependent periodic components during the "second plateau" in all bands but U. The authors initially interpreted this finding as caused by the presence of a stellar companion or a forming planet, but they strongly emphasized that future photometric observations will be essential to verify the driving mechanisms that they proposed. For this reason, we combined archival and new light curves to check if these oscillatory features are stable in time. In contrast to Kopatskaya et al. (2013), who for period analysis utilized detrended UBV data collected since 1995, in this study we use their BV data obtained since 1997 (HJD = 2,450,509, i.e., when the brightness level rested on that typical for the "second plateau") and the R_C I_C -filter data obtained since 2002. Afterward, we included data gathered with the Schmidt (BVR_C I_C ), RC80 (BV), and NSO (VI) telescopes, as well as the public-domain ASAS-SN Johnson-V data. The new and archival light curves were aligned to 0.002–0.005 mag by means of constant shifts to form uniform 19–23 yr long time series. To ensure linearity during period analysis, the light curves were transformed from magnitudes to flux units and were then normalized to unity at the mean brightness level of the complete 19–23 yr light curve.

Three period analysis techniques were used: as the light curves do not generally exhibit sine-like brightness variations, we decided to rely on the phase dispersion minimization (PDM) method (Stellingwerf 1978). We confronted these results with those obtained by means of the Fourier analysis, in which the mean standard errors of the amplitudes are conservatively evaluated using the bootstrap sampling technique (Rucinski et al. 2008). Finally, in order to check for stability of these oscillatory features in time, we used the weighted wavelet Z-transform (WWZ; Foster 1996), designed for analysis of unevenly sampled time series and available within the Vartools package (Hartman & Bakos 2016).

Results obtained by means of the PDM technique are shown in Figure 4(a). Only the significant parts of the periodograms, showing periods covered at least three times and longer than 100 days, are presented. The most significant peaks for BV filters are centered at 1707 ± 70 days. In spite of the formally inconclusive value (0.6) of the θ statistic, both the archival and the new BV filter phased data clearly show periodic behavior (Figure 4(b)). HJD₀ ^BV = 2,454,410—the best-defined minimum in BV filters that occurred at the end of 2007—was assumed during phase calculation.

Zoom In Zoom Out Reset image size

At first sight the PDM diagram obtained from R_C -filter data may appear to be inconclusive. First, because the primary ∼2000-day peak is poorly defined, the full extent of the long period lies outside the plotted portion of this periodogram. Second, the derived primary period is a multiple of the identified 497-day periodicity, which is also seen clearly in the R_C -filter data. Note that within the measurement uncertainty 497 days is indistinguishable from 502 days obtained from I_C I-filter data (see below) and 523 days obtained from the archival R_C I_C JHK-filter data (Kopatskaya et al. 2013). After rejecting the 2000-day peak, following the authors in Figure 4(b), we plot the R_C -filter light curve phased with a 1707-day period to examine the wavelength-amplitude evolution of this quasi-periodicity. We note that this peak is fairly well defined (although shifted to 1750 days) in the R_C -band PDM diagram as well.

For the same reasons as above, we adopted 1707 days for I_C -filter data phasing (Figure 4(b)). This peak is visible between the major ones at 1500 and 2000 days (θ = 0.4), which are the multiples of the dominant 502-day quasi-period (θ = 0.6). According to Figure 1 in Schwarzenberg-Czerny (1997), the false-alarm probability of this 502-day quasi-period is ≤few percent.

In Figure 4(c) we show the light curves phased with the 502-day quasi-period. HJD₀ ^I = 2454698—the best-defined minimum in the I_C filter that occurred in 2008 (288 days after the best-defined minimum in BV filters)—was assumed as the reference moment during phase calculation. In order to prepare these light curves, we applied a custom procedure to clear the original "second plateau" observations from the 1707-day quasi-periodic oscillation (QPO) variability: the specific shape of each light curve shown in Figure 4(b) was approximated by an ordinary seventh-to-ninth-order polynomial fit and then periodically subtracted. Thanks to our Piszkéstető data being of a higher precision, the presence of the 502-day component was for the first time directly confirmed in the V-band light curve.

The wavelet analysis of the entire 19–23-day long BVR_C I_C light curves confirms the above results: the WWZ spectra indicate the broad, fairly stable in time, primary ∼1700–2000-day period for BVR_C filters (Figures 4(d)–(f)) and the strong ∼502-day period for the I_C filter (Figure 4(g)) ¹⁴ Interestingly, signatures of the 502-day signal are noticeable in the form of a few isolated features in the archival and new BVR_C I_C data, although surprisingly in certain bands this quasi-period appears to evolve or even to be suppressed. We stress that even though WWZ is designed for analysis of unevenly sampled data, the resulting spectrum does strongly depend on photometric quality and data density. A mixture of these effects makes existing quasi periods impossible to disentangle at all times. This limitation allows us to firmly detect the 502-day QPO only in periods with good temporal sampling. However, these problems can be partially overcome, as described in the discussion of the color index variations (Section 3.2.2).

Finally, we calculated Fourier spectra to check the PDM and WWZ results and to investigate the relationship of the amplitudes (a) in the frequency (f) space (Figures 5(a)–(b)), which is carrying information about the nature of these small-scale oscillations. Except for the rough confirmation of the PDM results, we found that the ground-based spectra that "feel" the longer family of QPOs only show the stochastic flicker-noise nature characterized by a_f ∼ f^−1/2 (Press 1978). We will return to this issue in Section 4.6.

Zoom In Zoom Out Reset image size

3.2.2. Periodic Color Index Variations

The light curves themselves are affected by secular light changes, which in turn worsen the above-obtained PDM, WWZ, and Fourier results. Therefore, we decided to reexamine the above-obtained quasi-periods by means of the color index (CI) variations. In other words, analysis of light curves formed from the CIs can be treated as a counterpart of the usual whitening that is sensitive to the nonperiodic and is intrinsic to the disk's environment gray variability factors. The majority of these undesirable effects are expected to be removed, while the pure quasi-periodic variations driven by the not-yet-well-understood mechanisms should still be preserved.

We have performed time variability analysis of B − V, V − R_C , V − I_C , and R_C − I_C CIs with the PDM and WWZ techniques. To ensure homogeneity, we decided to use only the archival and the Schmidt telescope data. The PDM analysis confirmed the previous value of 502 days. We also found that the long-periodic variations weakened to a large extent. This is most visible in the new wavelet spectra dominated by the 502-day oscillation, which now appear as persistent and stable for the entire "second plateau" (Figures 6(a)–(c)). In Figures 6(d)–(e) we show associated archival and new CI light curves phased with the 502-day quasi-period and assuming HJD₀ ^I . Note that this approach reveals 502-day variations in the B − V light curve, but only in the part composed of the precise Schmidt data (Figure 6(e); see also Figures 4(d)–(e)). This analysis clearly shows that in the zero phase all CIs are significantly bluer than during the light maximum. This is in line with the CIs that can be directly inferred from light curves alone (Figure 4(c)), but opposite compared with the 1707-day period (Figure 4(b)), in which the respective CIs are redder at times when the disk is fainter.

Zoom In Zoom Out Reset image size

3.2.3. Short-term Variability as Seen by TESS

As already shown by Kopatskaya et al. (2013) and confirmed in Figures 4(b)–(c), the two QPOs observed in the "second plateau" show very different amplitude–wavelength dependencies. We also noted that these amplitudes evolve in time in our observations. To characterize this effect more profoundly, for each Johnson filter we determined the amplitudes by sinusoidal fits to the phase-folded (Figures 4(b)–(c)) light curves constructed from the archival (1997–2011) and from the new (2011–2020) data only. This approach minimizes the nonperiodic overlapping effects.

In the case of the 1707-day period, the amplitudes decrease with increasing effective wavelength of a filter: for the archival data we obtained 0.143(9), 0.122(6), 0.088(6), and 0.068(7) mag for BVR_C I_C filters, respectively. The errors shown in parentheses represent the 1σ uncertainty obtained from the least-squares fits. New data show the same well-defined amplitude–wavelength trend, but the resulting 0.066(9), 0.033(6), and 0.021(6) mag for BVR_C filters, respectively, clearly indicate (within 3σ) that the amplitudes are systematically becoming smaller in all bands, to the point that no variability has recently been detected in I_C band.

The amplitudes associated with the 502-day period are gradually increasing with the wavelength: no variability has been detected in B band, in both the archival and the new light curves. There is no evidence of variability in the archival V-band data, and only archival red and NIR data show significant variation—0.068(5) (R_C ), 0.131(5) (I_C ), 0.161(31) (J), 0.146(33) (H), and 0.130(33) mag (K), respectively. We find variation of 0.026(7), 0.057(6), and 0.080(7) mag for VR_C I_C filters, respectively, in the accurate Schmidt telescope data. Recent JHK data are too sparse to estimate the current amplitudes. We conclude that, unlike the 1707-day QPO, there is no obvious sign of time evolution of the amplitudes associated with the 502-day QPO.

Evolution of the color indices during the first post-outburst stages has already been investigated by Kopatskaya et al. (2013). The authors found that after the gradual color evolution along the extinction path in the phase of the exponential decay (1971–1985), during the "first plateau" (1985–1995), when the source became fainter than V ∼ 11.5 mag, the color index showed a "blueing effect," which can be observed in the young UX Ori–type objects. According to the authors, this effect has no longer been obviously present in the "second plateau."

Here we continue investigation of the color index evolution during the "second plateau." We utilize the archival data combined with the new one obtained in BVR_C I_C and BV filters with the Schmidt and RC80 telescopes, respectively. We show obtained results in Figure 7. Data obtained during individual years are marked by different colors and symbols.

Zoom In Zoom Out Reset image size

In our figures, the majority of the CI variations most closely follow the extinction path (dark solid line), which is calculated by our accretion disk model assuming the mean extinction law (R_V = 3.1; see also Section 4.1 for more details). Both the uncertainty related to the true level of I_C -band photometry and simplified assumptions about the disk photosphere radiation function are potential sources of the differences between the observed and synthetic color–magnitude diagrams. However, a more detailed look into the 2019 data (panels (b), (d), (f)) does reveal a different relationship. The same is valid for the 2020 (B − V) − V diagram (Figure 7(d)) and for the associated (V − r) − V and (V − i) − V diagrams. These trends cannot be explained only by variable accretion (represented by the red solid line)—these brief "blueing events" (similar to those observed during the first plateau) are currently observed when the target is at the minimum brightness during the 502-day quasi-period, and the CI variations related with the 1707-day QPO are relatively constant (see Figure 8 for illustration of mutual relations between both QPOs during 1997–2011). Searching for similar events, we also examined archival 2002–2011 data. Only the (B − V) − V diagrams obtained in 2005, 2006, and 2010 exhibit signatures of the CI reversal, but they are absent in the associated (V − R_C ) − V and (V − I_C ) − V diagrams.

Zoom In Zoom Out Reset image size

We also investigated color–magnitude diagrams from the 2019 data gathered simultaneously with TESS. The spacecraft coincidentally observed V1057 Cyg during the major brightness increase (phases 0.98–0.1 according to ephemeris adopted for the 502-day QPO). Thus, the associated diagrams show the same well-defined CI reversal evidence characteristic of the entire 2019 data set. In addition, we performed analysis of two specific color–magnitude diagrams, constructed from data obtained during the fainter and the brighter stages. Several (but not all) diagrams indicated variations along the extinction path, suggesting that the small-scale light changes noticed by TESS are just scaled-down counterparts of the major ones observed from the ground. However, given the limited precision of ground-based data and these relatively small brightness changes, respective correlation rank numbers are not high enough to confirm this behavior with high certainty.

We detected several emission and absorption lines in the spectra of V1057 Cyg. We used the NIST Atomic Spectra Database for the line identification (Kramida, A., Ralchenko, Yu., Reader, J., and NIST ASD Team 2020, NIST Atomic Spectra Database, ver. 5.8 ¹⁵ ). The lines detected in the spectra of V1057 Cyg are also listed in Appendix B, Tables 6 and 7.

3.5.1. Optical Spectroscopy

Classical FUors show several common optical spectroscopic characteristics: P Cygni profile of Hα, strongly blueshifted absorption lines, Li i absorption, and spectra similar to F/G supergiants/giants (Hartmann & Kenyon 1996; Audard et al. 2014). These spectroscopic features are also seen in our observations, and most of the features vary with time.

P Cygni profiles of several lines of hydrogen and metallic lines are found in the spectra of V1057 Cyg. The blueshifted absorption component of these profiles is formed by an outflowing wind (Hartmann & Kenyon 1996; Hartmann 2009; Herbig 2009; Reipurth & Aspin 2010). The strength of the blueshifted absorption component in the P Cygni profile of the Hα line is related to mass loss in the wind (Herbig et al. 2003). In our observations, P Cygni profiles of Hβ 486.2, Hα 656.3, and the Ca ii infrared triplet (849.8, 854.2, and 866.2 nm) lines are identified. Figure 9 shows examples for the observed P Cygni profiles. The blueshifted absorption component of all P Cygni profiles strongly varies with time. The high-velocity component of the wind was observed in all P Cygni profiles, and the highest-velocity component was extended to about −300 to −350 km s⁻¹ in 2018 December.

Zoom In Zoom Out Reset image size

The strength of the emission component of the lines with P Cygni profiles also varies with time. Although there is no tight correlation between the variation of the absorption and emission components in most lines, they show similar trends in the case of Hα: when the blueshifted absorption component of the Hα P Cygni profile was at the highest velocity (2018 December), the strength of the redshifted emission component was also the strongest, and vice versa (the weakest in 2020 August).

Strongly blueshifted absorption profiles caused by wind (Bastian & Mundt 1985; Herbig et al. 2003; Hartmann 2009; Miller et al. 2011) are also observed in V1057 Cyg. Some of the strongest examples are Fe ii 492.3 nm, Fe ii 501.8 nm, Mg i 518.3 nm, and the Na D doublet (588.9 and 589.5 nm), and these are plotted in Figure 10. All of the observed blueshifted absorption lines vary with time, and the variation trend is similar to that of the blueshifted absorption component of the P Cygni profiles (Figure 9). Among the observed blueshifted absorption lines, the Fe ii 501.8 nm and the Mg i 518.3 nm lines show the same velocity variation over time and thus likely originate from the same location in the structure.

Zoom In Zoom Out Reset image size

Several shell features are also found in the spectra of V1057 Cyg. A total of eight shell features in the range of 493–671 nm, showing similar velocity variations with time to the Li i 670.7 nm line, are selected. Four representative lines that show clear spectral profiles are presented in Figure 11: Ba ii 493.4 nm, Ti i 499.9 nm, Fe i 511.0 nm, and Li i 670.7 nm. Since various atomic lines show the same velocity distribution, the correlation between atomic properties (lower energy level E_i , upper energy level E_k , and transition probability A_ki ) and line profiles of shell features was investigated. However, no correlations between the line profiles and the different atomic parameters were found. All the detected shell features also vary with time during our observations. The highest-velocity and strongest absorption profile is detected in 2017 May (green line) when the width of the blueshifted absorption component of the P Cygni profile and the wind features are the narrowest (lower velocity). As noted in previous studies (Herbig et al. 2003; Herbig 2009; Kopatskaya et al. 2013), a weak emission component of the Li i 670.7 nm line was also observed in 2012 September.

Zoom In Zoom Out Reset image size

We also detected several forbidden emission lines in the spectra of V1057 Cyg, such as [N ii] 654.8 and 658.3 nm, [S ii] 671.6 and 673.1 nm, [O i] 630.0 and 636.3 nm, [O iii] 495.9 and 500.7 nm, and [Fe ii] 715.5 nm, which are rarely detected in classical FUors. Among the several forbidden emission lines, the relatively weak [N ii] 654.8 and 658.3 nm, [S ii] 671.6 and 673.1 nm, and [O iii] 495.9 and 500.7 nm lines are detected for the first time in the spectra of V1057 Cyg. The [S ii] emission line was previously found in only three known FUors: V2494 Cyg (Magakian et al. 2013), V960 Mon (Takagi et al. 2018; Park et al. 2020), and V346 Nor (Kóspál et al. 2020), and the [N ii] emission line was only found in V346 Nor (Kóspál et al. 2020).

These forbidden emission lines are generally associated with spatially resolved jets or outflows in Class II YSOs (Cabrit et al. 1990; Hartmann 2009). The [N ii], [S ii], and [O iii] emission lines are relatively narrow, and the peak velocity is located around systemic velocity (Figure 12). These emission lines are detected in most epochs, and their strengths also changed. The [O iii] emission lines are relatively stronger than the [N ii] and [S ii] lines. Most of the [O iii] lines are detected during our observations, except [O iii] 500.7 nm in 2015 December. The strength of the [S ii] emission lines is weaker than those of the [N ii] emission lines. The strength of [N ii] emission lines is very weak in 2020 August, and the [S ii] 673.1 nm line is not detected in 2018 October, and neither of the [S ii] emission lines was observed in 2020 August, indicating that the jet/outflow is also showing variability in time.

Zoom In Zoom Out Reset image size

In the case of T Tauri stars, the [O i] 630.0 nm emission line is often observed as two components: high velocity (a few hundred kilometers per second) and low velocity (a few tens of kilometers per second) (Hartigan et al. 1995; Hartmann 2009). In our observations, the [O i] 630.0 nm line shows two velocity components that are both greater than 91 km s⁻¹ wide. The high-velocity component can be formed by the outflowing wind (Hartmann 2009). The relatively higher velocity peaks are at around −140 to −213 km s⁻¹, and the relatively lower velocity peaks are at around −91 to −117 km s⁻¹ (Figure 13). The velocity variation of these components is similar to those of lower-velocity components of shell, wind, and P Cygni profiles. Therefore, the lower-velocity component of these lines can be formed at the same place of the structure. The strength of the forbidden emission lines varies slightly, but less than that of the wind features.

Zoom In Zoom Out Reset image size

In contrast with previous studies (Kenyon et al. 1988; Hartmann & Kenyon 1996; Hartmann 2009), we did not detect double-peaked line profiles in our observations.

3.5.2. Near-infrared Spectroscopy

We detected several absorption and emission lines in the NIR spectrum. Figure 14 shows the comparison between our NOTCam spectrum observed in 2020 August 29 (red) and that of IRTF (Connelley & Reipurth 2018) observed in 2015 June 26 (black). Similarly to Connelley & Reipurth (2018), we also detected Paβ 1.281 μm, Al i 1.312 and 1.315 μm, and strong water absorption bands, although our spectrum is not corrected well around 1.35 μm in the J band because of strong telluric absorption features. In the H band, the 19–4, 15–4, and 13–4 lines of the Br series are detected in broad absorption, and the [Fe ii] 1.533 and 1.644 μm lines are detected in emission. The Mg i 1.588 and 1.741 μm absorption lines are also detected. The Brγ 2.165 μm, Ti i 2.228 μm, and Ca i 2.265 μm lines are detected in absorption in the K band. The Brγ appeared as a weak P Cygni profile in the previous study (Connelley & Reipurth 2018), but it appeared as an absorption line in this study. The difference between the two spectra is the detection of the [Fe ii] emission lines and the shape of the CO overtone bandhead features. Emission lines are rarely detected in classical FUors. However, as in the optical spectra, we also detected a few [Fe ii] emission lines in the NIR spectrum. Compared to earlier observations of V1057 Cyg, the strength of the CO first-overtone bandhead feature appears to be the weakest in 2020 (see Section 4.4).

Zoom In Zoom Out Reset image size

While the long-term light curve of V1057 Cyg suggests a general decay of the accretion rate after the outburst peak in 1971, changing extinction toward the source might also play a role. In this section we attempt to separate the effects of variable accretion and extinction and study their long-term evolution quantitatively. Following the method we have successfully applied on several young eruptive stars (Kóspál et al. 2016, 2017; Ábrahám et al. 2018; Kun et al. 2019; Szegedi-Elek et al. 2020), we model the inner part of the system with a steady, optically thick and geometrically thin viscous accretion disk, whose mass accretion rate is constant in the radial direction (Equation (1) in Kóspál et al. 2016). We neglect any contribution from the star itself, assuming that all optical and NIR emission in the outburst originates from the hot accretion disk. We calculated synthetic SEDs of the disk by integrating the blackbody emission of concentric annuli between the stellar radius and R_out.

A fundamental input parameter of the model is the inclination of the accretion disk. Estimates in the literature, mainly based on SED fitting, range between 0° (pole-on) and 30° (for a review, see Gramajo et al. 2014). In order to derive a value based on observations, we analyzed the 1.3 mm continuum observations of Liu et al. (2018) obtained with the Submillimeter Array (SMA). Deconvolving the measured size of their continuum source (100 × 059, PA = 84°) by a beam of 087 × 050, PA = 76°, the resulting ratio of the minor and major axes implies an inclination of i = 62°. This result indicates a more edge-on view of the disk than previously thought. While this inclination value was derived from measurements of the whole disk, including both the outer cold regions and the hot inner disk, we will adopt it for the subsequent modeling of the accretion disk. This assumption is independently confirmed based on comparison of our Na i doublet spectra with those obtained from disk wind models by Milliner et al. (2019).

The outer radius of the accretion disk, another input parameter, mainly affects the mid-IR emission. We fixed it to R_out = 1 au, which matches the early L-band observations of V1057 Cyg in the 1970s. The inner radius of the disk, equal to the stellar radius, mainly influences the optical bands. However, we cannot discriminate between the cases of smaller stellar radius with higher line-of-sight extinction and those of larger radius with lower extinction using our broadband optical photometry. In order to circumvent this problem, we prescribed that the A_V value computed for 2020 August must comply with the A_V = 3.9 ± 1.6 mag proposed by Connelley & Reipurth (2018) based on an infrared spectrum taken in 2015. This constraint set R_in = 4.6 R_⊙.

With this model setup, only two free parameters remain: the product of the accretion rate × stellar mass $M\dot{M}$, and the line-of-sight extinction A_V . We calculated disk SEDs for a large range of $M\dot{M}$, and at each step the fluxes were reddened using a large grid of A_V values assuming the standard extinction law from Cardelli et al. (1989) with R_V = 3.1. Finding the best $M\dot{M}$–A_V combination was performed with χ² minimization, by taking into account all measured flux values between 0.4 and 2.5 μm. Preferentially we performed our modeling when both optical and infrared data were available for the same night, but we also included epochs when only JHK photometry was taken but optical data were available within 10 days; thus, interpolation in the optical fluxes was acceptable. The formal uncertainties of the data points were set to a homogeneous 5% of the measured flux value, which also accounted for possible differences among photometric systems. The model fits usually reproduced the measurements reasonably well, with typical reduced χ² values below 4.

The resulting temporal evolution of the accretion rate and extinction values, together with the V- and J-band light curves, is plotted in Figure 15. The initial decay of the source, between the outburst maximum and 1987, can be explained by an exponential drop of the accretion rate from 10⁻³ M_⊙ M_⊙ yr⁻¹ to ∼2.5 × 10⁻⁴ M_⊙ M_⊙ yr⁻¹, with an e-folding time of 4300 days (∼12 yr). During this fading phase (1971–1987), the extinction first slowly increased by ∼1 mag, while after 1983 it slightly decreased again, suggesting a rearrangement of the circumstellar structure, and/or a change in the dust size distribution in the line of sight, leading to a different extinction law. Then, between 1987 and 1993 both the accretion rate and the extinction stayed constant. In 1994–1995 A_V suddenly rose by ∼0.6 mag (probably causing the sudden drop of optical brightness at the same time). In the "second plateau" phase no long-term trend can be seen in $M\dot{M}$, and only a weak initial decay in A_V . On top of this relatively constant behavior in the "second plateau," correlated oscillations can be seen in the $M\dot{M}$ and A_V curves. These are probably due to the fact that the unusual optical−infrared color variations, caused by the superposition of two periodic processes of very different wavelength dependencies (Figure 8), cannot be simply reproduced by our simple analytical model, and thus these variations should not be overinterpreted. The current luminosity of the accretion disk is about 330 L_⊙, but its value depends on the disk inclination value. Since we adopted a more edge-on orientation than before in the literature, our inferred luminosity also increased. The current accretion rate of V1057 Cyg in our model, also slightly dependent on the inclination and the stellar radius, is about 10⁻⁴ M_⊙ M_⊙ yr⁻¹.

Zoom In Zoom Out Reset image size

In Figure 16 we plot the SED of V1057 Cyg at several epochs since the outburst. The optical and NIR points are from Figure 1, while longer-wavelength photometry was collected from different space-borne (IRAS, ISO, Spitzer, WISE, Akari, Herschel) or airborne (SOFIA) missions. The data points from Herschel and the AllWISE catalog were taken within a year; thus, we combined the two data sets into a single SED. The SOFIA spectra were smoothed and scaled to simultaneous SOFIA photometry.

Zoom In Zoom Out Reset image size

The gradual decrease in the short-wavelength part reflects the evolution of the hot inner accretion disk as modeled in Section 4.1. The difference between the 1993 and 1995 SEDs displays how the fading in 1995 became apparent first in the optical regime, while the NIR part stayed constant. The SEDs after 2003 ("second plateau" phase) were very similar; their slight differences reflect only the periodic behavior described in Section 3.2. Between 5 and 100 μm, V1057 Cyg also faded, although significantly less than at optical wavelengths (part of this flux drop might be related to the improving spatial resolution and thus smaller aperture size of the subsequent telescopes).

Based on a comparison of IRAS and ISO measurements, Ábrahám et al. (2004) claimed that below 25 μm the flux was variable, while at longer wavelength it remained constant. Extending the temporal baseline of this study with subsequent Spitzer, Akari, WISE, Herschel, and SOFIA measurements, almost a factor of 3 systematic fading was observed between IRAS and SOFIA at ∼25 μm. This fading was also seen at far-infrared wavelengths by comparing the Akari and Herschel photometric points to earlier IRAS and ISO results, although the fading was less pronounced (less than a factor of 2). Ábrahám et al. (2004) concluded that the outer part of the system, responsible for the long-wavelength SED, has an energy source different from the central star. However, our new results imply that the circumstellar medium does react to the changing irradiation by the central source, and thus the origin of the energy emitted by the envelope is more likely the outbursting star than an external source.

As described above in Section 3.5, we detected several wind features in the spectra of V1057 Cyg. The velocity of the blueshifted absorption component and the strength of the redshifted emission component of P Cygni profiles vary from year to year. The highest velocity of the blueshifted absorption component was observed in 2018 December in all P Cygni profiles (Figure 9) and wind features (Figure 10), and blueshifted velocity components of P Cygni profiles and wind features change similarly with time. From our observations, we can confirm that the year-to-year variabilities of strongly blueshifted absorption components of P Cygni profiles and wind features are similar to those observed by Herbig et al. (2003), suggesting variability over time in the strength of the wind.

The emission components of the Hα and Ca ii infrared triplet (IRT) P Cygni profiles are the strongest in 2018 December, while the other lines behave differently. All of the absorption and emission components of P Cygni profiles change with time, but there is no tight correlation between the two components, except for Hα (Section 3.5.1).

The shell features were variable during our observations, but the data do not suggest a well-defined trend. The shell features were the strongest in 2017 May, when the blueshifted components of wind and P Cygni profiles were the lowest velocity, and also when the system was close to the minimum light in that year (Figure 2). Since both depth and velocity change over time discontinuously, this variation over time can be interpreted as a rapidly changing wind or rotation of nonaxisymmetric components (Powell et al. 2012; Sicilia-Aguilar et al. 2020).

We also detected several forbidden emission lines: [N ii] 654.8 and 658.3 nm, [S ii] 671.6 and 673.1 nm, [O i] 630.0 and 636.3 nm, [O iii] 495.9 and 500.7 nm, and [Fe ii] 715.5 nm. These are rarely found in FUors: so far, only three known FUors show these properties (Magakian et al. 2013; Takagi et al. 2018; Kóspál et al. 2020; Park et al. 2020). In contrast, these lines are generally found in CTTSs as tracers of outflow or jets (Cabrit et al. 1990; Hartmann 2009). The lack of forbidden emission lines in FUors could be due to the lack of detailed spectroscopic studies and the small number of FUors known at this point. In addition, typically, the continuum of FUors is very bright, which makes it hard to detect the forbidden emission lines owing to the contrast. On several epochs, [O iii] 495.9 and 500.7 nm, [N ii] 654.8 and 658.3 nm, and [S ii] 671.6 and 673.1 nm lines are detected for the first time in the spectra of V1057 Cyg. All of the detected forbidden emission lines also vary with time, but less than the wind features. However, the variation of these emission lines suggests that any jets/outflows in the system also change with time.

The strength of the CO bandhead feature in V1057 Cyg decreased and the equivalent width (EW) increased in these epochs, according to the original studies (Mould et al. 1978; Hartmann & Kenyon 1987; Biscaya et al. 1997), and this trend continued in recent observations. Figure 17 shows the recent observations of the CO v = 2–0 2.293 μm first-overtone bandhead with the NOTCam (red) and the IRTF (black; connelley). We measured the EW of the CO feature from 2.293 to 2.317 μm (blue dashed line), which is the same region used by Biscaya et al. (1997) (see their Table 3). The EW of the NOTCam (27.03 ± 0.45 Å) and the IRTF (22.75 ± 0.39 Å) data were estimated with a Monte Carlo method. The EWs were measured 1000 times with random Gaussian errors multiplied by the observation errors. The standard deviation derived from all 1000 EW measurements was adopted as the uncertainty of the EW. Our results, together with values from the literature, are listed in Table 3. The measured EW is stronger in 2015 and 2020 than in 1986 (from 2.293 to 2.305 μm) and 1996, as the K-band magnitude decreases (Figure 1). We suggest that the weakened strength of the CO overtone bandhead features in our observation of V1057 Cyg is also caused by the decrease in brightness (Biscaya et al. 1997; Connelley & Reipurth 2018), which can then be related to decreasing mass accretion rate and disk midplane temperature. Our modeling of the disk (Figure 15) confirms the proposed explanations of decreasing brightness and therefore likely decreasing mass accretion rate and midplane temperature.

Zoom In Zoom Out Reset image size

With 23 yr of coverage of the "second plateau," in Section 3.2 we refined the values of the associated quasi-periods to 1707 ± 70 days and 502 ± 20 days. Precise Schmidt and RC80 telescope data enabled detection of the shorter period in the B and V bands for the first time. Similarly, we obtained the first marginal detection of the longer period in the I_C band. Furthermore, we obtained that during the light minimum associated with the 502-day period all CIs are becoming bluer, but we found the opposite (redder) for the 1707-day quasi-period. The amplitudes related to the longer period decrease, while the amplitudes related to the shorter period increase with increasing effective wavelength (BVR_C I_C ), respectively (Section 3.3). Moreover, our precise Schmidt telescope data revealed that as the time progresses, the amplitudes related to the 1707-day period are becoming smaller in all these filters. No significant evolution of amplitudes related to the 502-day period is observed.

In order to scrutinize possible mechanisms driving these quasi-periods, we searched for their signatures during the earlier post-peak epochs. We calculated the residual BVR_C light curves obtained by subtraction of the general trend during the exponential decay and the "first plateau." No significant peaks other than those closely related to the breaks in the data acquisition (340–410 days) were found (see also Clarke et al. 2005; Kopatskaya et al. 2013, for similar attempts). This suggests that both quasi-periods are not permanent features of V1057 Cyg but have a close relationship with the mechanism that led to the brightness drop in 1995–1996.

Two mechanisms driving these periodic variations have been proposed. Clarke et al. (2005) concluded that the erratic photometric variability observed in V1057 Cyg between 1996 and 2003 is associated with the fallback of dusty material to small radii and the subsequent passage of dust condensations across the line of sight to the inner accretion disk (Section 4.5.2 in their paper). However, at the time of this study the periodic behavior was unknown. Enriched with this knowledge, we conclude that this event can be uniquely associated with the 1707-day period. We base this result on respective amplitudes obtained in Section 3.3, which are decreasing with increasing effective wavelengths of consecutive filters roughly in line with the mean extinction law (Figure 18). Assuming Keplerian rotation of this inhomogeneous region, this dust condensation scene could be located 1.9–2.8 au from the 0.3–1 M_⊙ mass star, respectively. In these circumstances the recently observed amplitude decrease could also be explained by dispersion of this dusty region in time and/or lower dust production rates, as the disk wind is the subject of weakening owing to the decreasing accretion rate. On the other hand, this dimming and its subsequent evolution could also be caused by the disk warp localized at 1.9–2.8 au, as the disk is seen more edge-on (Section 4.1) than previously thought. According to the unified models of innermost disk warps (McGinnis et al. 2015), the maximum warp height is 20%–30% of the disk radius at which it originates (i.e., up to 0.5–0.8 au assuming that this model scales to these distances), and it may vary by 10%–20% during a single rotation.

Zoom In Zoom Out Reset image size

According to Güver & Özel (2009), N_H [g cm⁻²] = 0.00367A_V . Assuming that whatever feature localized at 2.2 au (for the 0.5 M_⊙ star) causing the extinction changes occupies a quarter orbit (based on Figure 4(b)), its linear size measured along the disk plane is about 3.5 au. Using an arbitrary height of 0.8 au, the mass of this structure would be 0.0004 M_⊕ for ΔA_V = 1 mag during the "second plateau," at most (see the bottom panel in Figure 15). Considering that V1057 Cyg is accreting about 33 Earth masses per year, this represents a negligible fraction of the total disk mass and an order of magnitude less than estimated for an analogous phenomenon in V582 Aur by Ábrahám et al. (2018).

In spite of the very limited number of spectra, we decided to search for correlations between the spectroscopic and the known photometric variability. We tentatively find correlation between the 1707-day QPO and the absorption components of the wind lines. In Figures 19(a)–(b) we present Hα and Hβ P Cygni profiles as a function of phase calculated for this QPO (Figure 4(b)). The RVs measured consistently at a depth of 2/3 of the virtual line bisector of the first (usually the deepest) absorption peak (Herbig 2009) are possibly periodic (i.e., 1707 days), with mean amplitude of ΔRV = 5.4 ± 1.6 km s⁻¹ and mean velocity of γ = −81.4 ± 1.2 km s⁻¹ The reason for measuring this peak is that it represents the last interaction of the inner disk light with the dusty environment that modifies it and is probably the most important factor in determining what kind of photometric variations will ultimately be seen by the observer. With a limited number of spectra, the uncertainties are too large for meaningful considerations of the distinct lines. If the above spectroscopic–photometric connection is real, the obtained γ value weighs against the hypothesis of occultations caused by a disk warp, as those should produce RV variations around the mean systemic velocity.

Zoom In Zoom Out Reset image size

Similar analysis performed for shell lines and the broad and narrow emission peaks in forbidden lines do not reveal variations correlated with phase. However, the lower- and higher-velocity components of the double-peaked broad emission peaks appear with change in accordance to the 1707-day quasi-period.

With regard to the shorter period, the most plausible interpretations invoked binarity or even planet(s) forming and/or tidally disrupting in the very inner disk (Lodato & Clarke 2004; Clarke et al. 2005). This possibility has been considered more seriously once new observations revealed that the 502-day period correlates with RV changes of the emission component of the lithium line (Kopatskaya et al. 2013). The reversed color index behavior associated with the 502-day period is similar both to the "blueing effect" observed in UXors and also to the effect caused by the rotation of a locally warmer plasma bubble, as proposed for FU Ori's inner disk (Siwak et al. 2018). If the latter scenario is valid and the disk rotation is Keplerian, the "hot spot" in V1057 Cyg's disk would have to be localized at the distance of 1 au from the 0.5 M_⊙ star (Gramajo et al. 2014) to be responsible for the 502-day QPO. Using the model presented in Siwak et al. (2018), but assuming blackbody instead of supergiant spectral intensities, ¹⁶ R_in = 4.6 R_⊙ and i = 62° (Section 4.1), we were able to reproduce the amplitudes observed in R_C I_C JHK filters at the same time by Kopatskaya et al. (2013) and listed in Section 3.3. A hot spot of T_eff ≈ 3500–4000 K (vs. ≈800 ± 100 K predicted for these disk annuli assuming steady accretion) with 30° azimuthal extent and 200–240 R_⊙ (0.93–1.12 au) is necessary to qualitatively reproduce these observations (Figure 20). The obtained hot spot size and temperature are in principle consistent with those numerically obtained for luminous optically thick shocks on the circumplanetary disks around giant forming planets (Szulágyi & Mordasini 2017). It is debatable whether such a shock could produce the associated RV variations observed in the Li line emission component by Herbig (2009) and Kopatskaya et al. (2013). The single piece of evidence for it can be found in our 2012 spectrum at +5 km s⁻¹, which is very different from those reported by the authors (from −9 to −19 km s⁻¹). We found no correlations between the 502-day quasi-period and shell or P Cyg or forbidden line profiles, but this may be due to the poor temporal coverage. The weakness of the hot spot scenario is that it fails to explain why the shorter QPO became observable simultaneously with the longer one. Disk instabilities ignited during the enhanced accretion are not well understood and could be the cause, or possibly the disk fragmentations that led to the FUor phenomenon (Vorobyov & Elbakyan 2018). Hot spot structures often attributed to forming planets have been directly imaged in disks of young stars at distances from a dozen to several dozen astronomical units (Biller et al. 2014; Reggiani et al. 2014).

Zoom In Zoom Out Reset image size

The 502-day quasi-periodic phenomenon can also be caused by obscuration of certain disk annuli, i.e., those mostly emitting in the red and NIR bands, by a dust cloud. In order to reproduce the observed amplitudes, we first calculated magnitudes for an unshaded disk assuming steady accretion and then compared the results with the synthetic magnitudes calculated assuming total eclipsing of certain disk annuli, contained within the a priori chosen azimuthal angle of 120°. This parameter has negligible effect on the final result. We approximated this phenomenon by setting T_eff = 0 in the shaded area. We obtained that a substantial part of the inner disk contained roughly between 20 and 70 R_⊙ must be periodically obscured to produce the observed effect. As the Keplerian period at 45 R_⊙ is about 50 days, a cold spot in a disk photosphere or a dusty cloud corotating with the disk on a circular orbit can immediately be excluded. Occultation of the innermost disk annuli as in the 1707-day scenario proposed above does not apply because it would produce a strong signal at UBV wavelengths (in fact, only barely detected in BV bands). The remaining possibilities are a dusty cloud—a remnant of the envelope—rotating on an inclined and eccentric orbit, and a "dust streamer" structure caused by interaction of forming planets and elevating a substantial amount of dusty material high above the disk midplane (Loomis et al. 2017). However, the first possibility seems to be unlikely, as this cloud would likely be absorbed after the first attempt to break through the protoplanetary disk midplane. A dust streamer could act as an occulting screen of certain disk annuli only, although the mechanism leading to periodic and continuous brightness variations is not clear. Thus, based on the data available so far, we conclude that obscuration scenarios considered above fail to explain the 502-day period and its spectral properties in a manner as self-consistent as a locally heated disk at 1 au.

While the nature of the quasi-periodic light variations observed from the ground since 1997 seems to be now better understood, the nature of these observed from space is still not clear to us. We obtained different Fourier spectrum slopes for the ground-based (a_f ∼ f^−1/2) and TESS data (a_f ∼ f⁻¹). These different relationships could be initially interpreted by combination of extinction, accretion, and light-scattering processes. We arrive at this conclusion owing to the lack of time-coherent QPOs in the TESS light curve of V1057 Cyg, similar to those observed in the disk-only FU Ori and V646 Pup (Siwak et al. 2018, 2020). This lack could be due to strong reprocessing of time-coherent inner disk light variations on their way to the observer. On the other hand, analysis of the accretion-dominated public-domain light curves of FU Ori itself and highly accreting CTTSs obtained during several seasons by several spacecraft revealed that the steeper, random walk Brownian nature observed by TESS during its 56-day run can also be an effect of a particular realization of stochastic accretion-related processes in the inner disk. Longer or more frequent observing runs with photometric precision provided by space telescopes are required to clarify these assumptions.