A study of the photometric and spectroscopic variations of the prototypical FU Orionis-type star V1057 Cyg

Among the low-mass pre-main sequence stars, a small group called FU Orionis-type objects (FUors) are notable for undergoing powerful accretion outbursts. V1057 Cyg, a classical example of an FUor, went into outburst around 1969-1970, after which it faded rapidly, making it the fastest fading FUor known. Around 1995, a more rapid increase in fading occurred. Since that time, strong photometric modulations have been present. We present nearly 10 years of source monitoring at Piszk\'estet\H{o} Observatory, complemented with optical/near-infrared photometry and spectroscopy from the Nordic Optical Telescope, Bohyunsan Optical Astronomy Observatory, Transiting Exoplanet Survey Satellite, and the Stratospheric Observatory for Infrared Astronomy. Our light curves show continuation of significant quasi-periodic variability in brightness over the past decade. Our spectroscopic observations show strong wind features, shell features, and forbidden emission lines. All of these spectral lines vary with time. We also report the first detection of [S II], [N II], and [O III] lines in the star.


INTRODUCTION
Photometric and spectroscopic monitoring of premain sequence (PMS) stars over a broad spectral range is crucial to understand the mechanisms leading to the formation of stars and ultimately planets. A small, but spectacular class of low-mass young stars are known as FU Orionis-type stars (FUors), referring to the nova-turies (Paczynski 1976;Lin & Papaloizou 1985;Kenyon et al. 1988;Kenyon & Hartmann 1991a;Bell et al. 1995;Turner et al. 1997;Audard et al. 2014;Kadam et al. 2020).
The members of this group, currently about 30 objects (Audard et al. 2014) show very similar optical spectra: F or G supergiants with wide absorption lines, Hα P Cygni profiles, shell components, and strong Li I 670.7 nm absorption. During a FUor eruption, the disk outshines the luminosity of the central star. Assuming that the bolometric luminosity calculated from the observed spectral energy distribution (SED) is dominated by the accretion luminosity, the accretion rate during the FUor stage can directly be obtained. Observations showed that the accretion rate rises from the average rate of a typical T Tauri star (10 −9 -10 −7 M yr −1 ) up to 10 −5 -10 −4 M yr −1 in only a few months (Hartmann & Kenyon 1996).
V1057 Cyg became the second identified FUor in 1969, when it brightened by 6 mag in the V -band (Welin 1971a,b). The source is located in the North America Nebula (NGC 7000), which, together with the Pelican Nebula (IC 5070) form a large HII region (Wendker 1983;Rebull et al. 2011). Previous distance estimates for these regions vary between 520 and 700 pc (Laugalys et al. 2006;Skinner et al. 2009;Fischer et al. 2012). In a recent work, Kuhn et al. (2020) determined new distances for the members of the North America Nebula using Gaia DR2 astrometry (Gaia Collaboration et al. 2018). They found that the main parts of the North America and Pelican Nebula are located at ∼795 pc, however V1057 Cyg, as a part of a smaller group of stars is located somewhat farther away. In this paper we adopt the Gaia DR2 distance value of 897 pc from Bailer-Jones et al. (2018a,b) which was specifically determined for V1057 Cyg. Herbig (1977) studied V1057 Cyg in detail both photometrically and spectroscopically. They concluded that before its eruption, the object had shown the properties of a classical T Tauri-type star (CTTS). They also characterised a 1. 0 × 1. 5 ring-like nebula which appeared around the object after the outburst. Further observations showed that the ring faded with the central star in the following years, but its structure remained unchanged. This indicated that the ring was a reflection nebula: a structure already present before the eruption of V1057 Cyg, illuminated by the central source, and not material that had been blown out during the eruption.
Three decades later Herbig et al. (2003) presented another detailed spectroscopic study focusing on this star. These high-resolution spectra, taken in 1996-2002, confirmed some of the previously observed features, such as the 'doubling' of low-excitation absorption lines which became more apparent between the 1980's and 1994. In this subsequent study, Herbig (2009) pointed out that V1057 Cyg has a long-lasting, high-velocity wind, which manifests itself through strong blueshifted absorption components at various optical lines.
The last photometric analysis of V1057 Cyg was performed by Kopatskaya et al. (2013), who demonstrated that immediately after reaching the light maximum in 1970, the light curve started an exponential brightness decline until ∼1985, when the so called 'first plateau' phase started and lasted for about 10 years. After that, the source faded by ∼0.5-1 mag in the optical within a year, and started to show quasi-periodic variations. The authors found that the variations could be characterized with two different periods: a longer period 1631±60 d, dominating the BV R data, and a shorter one 523 ± 40 d, dominating the IJHK data. They initially concluded that these fluctuations reflected the binary nature of V1057 Cyg, which has also been proposed as a possible mechanism leading to enhanced accretion and the FUor phenomenon (e.g. Bonnell & Bastien 1992;Bell et al. 1995). Interestingly, using non-redundant aperture-masking interferometry, Green et al. (2016) detected a faint companion star of V1057 Cyg, located at a projected separation of 58 mas with a brightness difference of ∆K = 3.3 mag. Its distance from V1057 Cyg suggests that it could have triggered the original outburst with a close fly-by encounter (Vorobyov et al. 2021). Connelley & Reipurth (2018) published a nearinfrared (NIR) spectroscopic survey including V1057 Cyg with observations from 2015, the latest NIR spectroscopic data of the source. They concluded that the CO absorption band was much weaker than in 1986. In contrast, the first high-resolution NIR spectroscopic observations of Hartmann & Kenyon (1987a) showed that the CO features have not changed much compared to Mould et al. (1978). Biscaya et al. (1997) showed that the CO features became weaker in 1996 than in 1986 (Hartmann & Kenyon 1987a) and interpreted that this weakening might be related to the brightness decline in 1995.
The infrared excess emission apparent in the SED of V1057 Cyg is due to a flared disk and envelope geometry (Kenyon & Hartmann 1991b). The presence of an envelope was also confirmed by Green et al. (2006) based on 5 − 35 µm Spitzer/IRS observations with an estimated radius of 7000 au. Zhu et al. (2008) modeled the dust from Spitzer/IRS observations and found that an envelope typical of protostars is required for V1057 Cyg to match the observations. Green et al. (2013) found the observed Herschel spectra were generally brighter than model predictions, which indicated an underestimate of the large scale reservoir of cold dust surrounding FUors. These works also suggested the idea of a large bipolar cavity in the envelope. Fehér et al. (2017) surveyed northern hemisphere FUors with the Plateau de Bure Interferometer (PdBI) and the IRAM 30 m telescope. Based on 13 CO observations, they found a rotating envelope around V1057 Cyg which is roughly spherical with a radius of 5 (3000 au) and a total circumstellar mass of 0.21 M .
Despite significant fading, the last visual spectrum of V1057 Cyg obtained in 2012 by Lee et al. (2015) did not resemble that of a CTTS, thus, further monitoring is key in tracing the gradual return of V1057 Cyg to quiescence. We have occasionally observed our target in optical and infrared bands since 2005, but intensified our monitoring after 2011, due to increased telescope time.
We describe the new observations and our reduction methods in Section 2. Results obtained from the data analysis are presented in Sec. 3 and discussed in Sec. 4. We summarize our findings in Sec 5.

Ground-based optical photometry
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: 0. 306). The 60/90/180 cm Schmidt telescope, equipped with a 4096 × 4096 pixel Apogee Alta U16 CCD camera (pixel scale: 1. 027), was used in 2011-2019. In each of the BV R 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 0. 55 and a field-of-view of 18. 8×18. 8 for a 2048×2048 pixel FLI PL230 CCD camera. We obtained three images per night in BV g r i filters.
The frames were calibrated for bias, dark, and flatfield in the standard fashion. Photometry of V1057 Cyg and 12 comparison stars in its 8 vicinity was extracted using an aperture radius of 4. 1 and sky annulus between 10. 3 -15. 4 for RCC and Schmidt frames, and 5. 5 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(Kóspál et al. , 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. Except of 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 9. 2×9. 0 field of view and a 0. 537 pixel scale. After the standard reduction steps for bias, dark, and flatfield 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 BV I filters, operated remotely through Skynet. The calibration procedures and comparison stars were the same as above, but only the V I NSO filter data was 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 6. 4×6. 4 and pixel scale of 0. 21. The Bessel BV R 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 − 30 s. After the standard CCD reduction steps, we obtained aperture photometry using an aperture radius of 3. 2 and a sky annulus between 6. 4 and 8. 6. 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 Fig. 1 and 2, and listed in Tab. 4 in the Appendix The typical uncertainty of our measurements is 0.03 mag in B and 0.01 mag in all other filters.

Space-based optical and infrared photometry
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 1 . The calibrated full-frame images were processed in two main steps using the FITSH package (Pál 2012). Firstly the plate solution was derived based on the Gaia DR2 catalogue -details of this complex procedure are described by Pál et al. (2020). As the 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 filters throughputs. By examining various TESS fields observed in the first two sectors we found that the RMS of our zerolevel 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, what caused three additional 1.2-1.3 d breaks 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 Fig. 3.
We complemented our work with data from the Widefield 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 & et al. 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 2 . The corrected WISE data are shown in Fig. 1 and listed in Tab. 5 in the Appendix.

Near-infrared photometry
We obtained near-infrared 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 5-point dither pattern in order to enable proper sky subtraction. The total integration time was typically 1 min 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 3 , as well as our own IDL routines. Data reduction steps included sky subtraction, flat-fielding, registration, and coadding 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-ofview (pixel scale: 0. 234). We obtained 9 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 aperture radius of 3. 3 and a sky annulus between 6. 6 and 9. 4. 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, Tab. 4.

Optical spectroscopy
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 fibre with a larger entrance aperture of 2. 5 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 flatfield, the ThAr lamp spectrum was used for wavelength calibration, and continuum fitting was performed by continuum task. Finally, heliocentric velocity correction was applied by the rvcorrect task and the published radial velocity of V1057 Cyg (−16 km s −1 ; Herbig et al. 2003).
As no telluric standard stars were observed neither for FIES nor BOES, we performed the telluric correction using the molecfit software Kausch et al. 2015) by fitting the telluric absorption bands of O 2 and H 2 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 Tab. 1.

Near-infrared spectroscopy
On 2020 August 29, we used the NOTCam on the NOT to obtain new near-infrared spectrum of V1057 Cyg and Iot Cyg (A5 V) as our telluric standard star in the JHK s bands. We used the lowresolution camera mode (R=2500) with ABBA dither positions, and exposure times ranged from 25 to 60 seconds (Tab. 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 (Sec. 4.1).
3. RESULTS AND ANALYSIS

Light curves
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;Welin 1975Welin , 1976Landolt 1975Landolt , 1977Simon 1975;Simon et al. 1982;Simon & Joyce 1988;Kenyon & Hartmann 1991b;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 utilised in these two data sets.
We found systematic differences between the two sets of photometry, which may be due to different apertures, filters, detectors throughputs, and different comparison stars used. For plotting purposes, we shifted our B band light curves by +0.12 mag, V band by +0.08 mag, R C band by +0.05 mag and I C band light curves by −0.14 mag to be consistent with the earlier papers. In Appendix A Tab. 4 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 Fig. 1, while in Fig. 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 Tab. 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 nearly constant level in all used filters. This was the so called 'first plateau' phase, which lasted until 1995. As mentioned in Section 1, U BV 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 Fig. 1) stopped in 1997 (Kolotilov & Kenyon 1997;Ibragimov 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).

Period analysis
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 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 utilised detrended U BV data collected since 1995, in this study we use their BV data obtained since 1997 (HJD = 2450509, 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. Afterwards, we included data gathered with the Schmidt Optical and infrared light curves of V1057 Cyg. We complemented our light curves with optical and infrared data prior to 2012 from Mendoza (1971); Rieke et al. (1972); Welin (1975Welin ( , 1976; Landolt (1975Landolt ( , 1977; Simon (1975); Simon et al. (1982); Simon & Joyce (1988); Kenyon & Hartmann (1991b); Kopatskaya et al. (2013); Green et al. (2016).
(BV R C I C ), RC80 (BV ) and NSO (V I) telescopes, as well as the public-domain ASAS-SN Johnson-V data.
The new and the archival light curves were aligned to 0.002-0.005 mag by means of constant shifts to form uniform 19-23 years long time series. To ensure linearity during period analysis, the light curves were transformed from magnitudes to flux units, and were then normalised to unity at the mean brightness level of the complete 19-23 years 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 Fig. 4a. Only the significant parts of the peri-odograms, 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 θ statistic, both the archival and the new BV -filter phased data clearly show periodic behaviour (Fig. 4b).
HJD BV 0 = 2454410 -the best defined minimum in BVfilters that occurred at the end of 2007 -was assumed during phase calculation.
At first sight the PDM diagram obtained from R Cfilter data may appear to be inconclusive. First, because the primary ∼2000 d 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 d periodicity, which is also seen clearly in the R C -filter data. Note that within the measurement uncertainty, 497 d is indistinguishable from 502 d obtained from I C I-filter data (see below) and 523 d obtained from the archival R C I C JHK-filter data (Kopatskaya et al. 2013). After rejecting the 2000 d peak, following the authors in Fig. 4b, we plot the R Cfilter light curve phased with 1707 d period to examine the wavelength-amplitude evolution of this quasiperiodicity. We note that this peak is fairly well defined (although shifted to 1750 d) in the R C -band PDM diagram as well. For the same reasons as above, we adopted 1707 d for I C -filter data phasing (Fig. 4b). This peak is visible between the major ones at 1500 and 2000 days (θ = 0.4), which are the multiples of the dominant 502 d quasiperiod (θ = 0.6). According to Fig. 1 in Schwarzenberg-Czerny (1997), the false alarm probability of this 502 d quasi-period is ≤ few %.
In Fig. 4c we show the light curves phased with the 502 d quasi-period. HJD I 0 =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 d QPO variability: the specific shape of each light curve shown in Fig. 4b was approximated by ordinary 7-9th order polynomial fit and then periodically subtracted. Thanks to our Piszkéstető data being of a higher precision, the presence of the 502 d component was for the first time directly confirmed in the V -band light curve.
The wavelet analysis of the entire 19-23 long BV R C I C light curves confirms the above results: the WWZ spectra indicate the broad, fairly stable in time, primary ∼1700-2000 d period for BV R C -filters ( Fig. 4d-f) and the strong ∼502 d period for the I C -filter ( terestingly, signatures of the 502 d signal are noticeable in the form of a few isolated features in the archival and new BV R 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, data density and a mixture of these effects makes existing quasi-periods impossible to disentangle at all times. This limitation allows us to firmly detect the 502 d 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 (Sec 3.2.2). Finally we calculated Fourier spectra to check the PDM and WWZ results and to investigate the relation- Figure 5. Amplitude-frequency spectra in log-log scale represented by solid (BRC-filters) and dotted-dashed lines (V IC-filters), calculated from the 'second plateau' data. The amplitude errors are marked by small dots; the significant peaks are located to the left from ≈0.01 c d −1 . The short marks indicate the frequencies corresponding to the periods determined by means of the PDM method. No dominant period can be indicated in TESS spectrum (panel c). The red flicker-noise spectrum slope is indicated by two parallel dashed lines: they show a f ∼ f −1/2 relation for the ground-based, and a f ∼ f −1 for TESS data. ship of the amplitudes (a) in the frequency (f ) space ( Fig. 5ab), 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 quasi-periodic oscillations (QPOs) only, show the stochastic flicker-noise nature characterised by a f ∼ f −1/2 (Press 1978). We will return to this issue in Sec. 4.6.

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 non-periodic and intrinsic to the disk's environment gray variability factors. The majority of these undesirable effects is expected to be removed, while the pure quasi-periodic variations driven by the not yet wellunderstood 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 PDM and WWZ technique. To ensure homogeneity, we decided to use only the archival and the Schmidt-telescope data. The PDM analysis confirmed the previous 502 d value. 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 d oscillation, which now appear as persistent and stable for the entire 'second plateau' (Fig. 6abc). In Fig. 6d-e we show associated archival and new CI light curves phased with the 502 d quasi-period and assuming HJD I 0 . Note that this approach reveals 502 d variations in the B − V light curve, but only in the part composed of the precise Schmidt data (Fig. 6e, see also in Fig. 4de). This analysis clearly shows that in the zero phase, all CI's are significantly bluer than during the light maximum. This is in line with the CIs that can be directly inferred from light curves alone (Fig. 4c), but opposite compared with the 1707 d period (Fig. 4b), in which the respective CIs are redder at times when the disk is fainter.

Short-term variability as seen by TESS
To gain insight into the variability occurring on the time scales of hours and days, we performed Fourier analysis of TESS data obtained with 30 min sampling. In accordance with the visual inspection of the light curve itself (Fig. 3), the amplitude-frequency spectrum does not show any dominant peaks (Fig. 5c). We also found that this spectrum exhibits a Brownian randomwalk, described by a f ∼ f −1 (Press 1978).
We also performed wavelet analysis of the TESS data. We do not report these results here, as the analysis is strongly affected by the previously mentioned (Sec. 2) six breaks in the data acquisition, which have a duration comparable to the characteristic time scale of observed light changes.

Amplitude-wavelength dependency of the QPOs
As already shown by Kopatskaya et al. (2013) and confirmed in Fig. 4bc, 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 ( Fig. 4bc) light curves constructed from the archival (1997-2011) and from the new (2011-2020) data only. This approach minimizes the non-periodic overlapping effects.
In the case of 1707 d 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 BV R C I C filters, respectively. The errors shown in parentheses represent the 1σ uncertainty obtained from the least-square 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 BV R 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 d period are gradually increasing with the wavelength: no variability has been detected in B-band, both in the archival and the new light curves. There is no evidence of variability in the archival V -band data, and only archival red and near-infrared data show significant variation -0.068(5) (R C ), 0.131(5) (I C ), 0.161(31) (J) 0.146(33) (H) 0.130(33) mag (K), respectively. We find variation of 0.026 (7), 0.057(6) and 0.080(7) mag for V R 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 d, there is no obvious sign of time-evolution of the amplitudes associated with the 502 d QPO.

Color-magnitude diagrams
Evolution of the color indices during the first postoutburst stages has already been investigated by Kopatskaya et al. (2013). The authors found that after the gradual colour evolution along the extinction path in the phase of the exponential decay (1971)(1972)(1973)(1974)(1975)(1976)(1977)(1978)(1979)(1980)(1981)(1982)(1983)(1984)(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 colour index evolution during the 'second plateau'. We utilize the archival data combined with the new one obtained in BV R C I C and BV -filters with the Schmidt and RC80 telescopes, respectively. We show obtained results in Fig. 7. Data obtained during individual years are marked by different colors and symbols.
In our figures, the majority of the CI variations most closely follow the extinction path (dark continuous line), which is calculated by our accretion disc model assuming the mean extinction law (R V = 3.1, see also Sec. 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 different relationship. The same is valid for the 2020 (B − V ) − V diagram (Fig. 7d) and for the associated (V − r) − V and (V − i) − V diagrams. These trends cannot be explained only by variable accretion (represented by continuous red 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 d quasi-period and the CI variations related with the 1707 d QPO are relatively constant (see Fig. 8 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 d QPO). Thus the associated diagrams show the same well-defined CI reversal evidence characteristic of the entire 2019 dataset. In addition, we performed analysis of two specific colormagnitude 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.

Spectroscopy
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 (version 5.8) 6 ). The lines detected in the spectra of V1057 Cyg are also listed in Appendix B, Tab. 6, and Tab. 7.

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 6 https://www.nist.gov/pml/atomic-spectra-database 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. Fig. 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 ∼ −350 km s −1 in 2018 December.
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 Fig. 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 (Fig. 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.
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 as the Li I 670.7 nm line, are selected. Four representative lines which show clear spectral profiles are presented in Fig. 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 be-  tween 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 September 2012. We also detected several forbidden emission lines in the spectra of V1057 Cyg, such as [N II] (Takagi et al. 2018;Park et al. 2020), and V346 Nor , and the [N II] emission line was only found in V346 Nor .
These forbidden emission lines are generally associated with spatially resolved jets or outflows in Class II YSOs (Cabrit et al. 1990;Hartmann 2009  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 km s −1 ) and low-velocity (a few tens of km s −1 ) (Hartigan et al. 1995;Hartmann 2009). In our observations, [O I] 630.0 nm line shows two velocity components which are both greater than 91 km s −1 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 −1 , and the relatively lower velocity peaks are at around −91 to −117 km s −1 (Fig. 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.
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.

Near-infrared Spectroscopy
We detected several absorption and emission lines in the near-infrared spectrum. Fig. 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, 1.315 µm, and strong water ab- sorption 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, 1.644 µm lines are detected in emission. The Mg I 1.588, 1.741 µm absorption lines are also detected. The Brγ 2.165, Ti I 2.228, 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 over- tone 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 near-infrared 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 Sec. 4.4).

Accretion disk modeling
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 towards 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Á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 (Eq. 1 in Kóspál et al. 2016). We neglect any contribution from the star itself, assuming that all optical and near-infrared 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 analysed 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 (1. 00×0. 59, PA=84 • ) by a beam of 0. 87×0. 50, 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 1970's. 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 lineof-sight extinction as opposed to larger radius with lower extinction using our broad-band 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Ṁ , and the line-of-sight extinction A V . We calculated disk SEDs for a large range of MṀ , and at each step the fluxes were reddened using a large grid of A V val-ues assuming the standard extinction law from Cardelli et al. (1989) with R V = 3.1. Finding the best MṀ -A V combination was performed with χ 2 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 χ 2 values below four.
The resulting temporal evolution of the accretion rate and extinction values, together with the V and J-band light curves, are plotted in Fig. 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 −3 M M yr −1 to ∼2.5×10 −4 M M yr −1 , with an e-folding time of 4300 days (∼12 yr). During this fading phase (1971)(1972)(1973)(1974)(1975)(1976)(1977)(1978)(1979)(1980)(1981)(1982)(1983)(1984)(1985)(1986)(1987) the extinction first slowly increased by ∼1 mag, while after 1983 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-95 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Ṁ , 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Ṁ 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 (Fig. 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 −4 M M yr −1 .

Spectral energy distribution
In Fig. 16 we plot the spectral energy distribution of V1057 Cyg at several epochs since the outburst. The optical and near-infrared points are from Fig. 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.
The gradual decrease in the short wavelength part reflects the evolution of the hot inner accretion disk as modeled in Sect. 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 near-infrared part stayed constant. The SEDs after 2003 ('second plateau' phase) were very similar; their slight differences reflect only the periodic behavior described in Sect. 3.2. Between 5 µm 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)  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 µma. This fading was also seen at farinfrared wavelengths by comparing the Akari and Herschel photometric points to earlier IRAS and ISO, 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.

Optical spectroscopy
As described above in Sec. 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 (Fig. 9) and wind features (Fig.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 variability 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 IRT P Cygni profiles are the strongest in 2018 December, while the other lines behave differently. All of the absorption and emission component of P Cygni profiles change with time, but there is no tight correlation between the two components, except for Hα (Sec. 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 (Fig. 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 non-axisymmetric components (Powell et al. 2012;Sicilia-Aguilar et al. 2020).
We also detected several forbidden emission lines:  (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 classical T Tauri stars 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 due to the contrast. On several epochs, [O III] 495.9, 500.7, [N II] 654.8, 658.3, and [S II] 671.6, 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.

Variation of the CO first overtone bandhead
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 1987b;Biscaya et al. 1997), and this trend continued in recent observations. Fig. 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 Tab. 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 (Fig. 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 (Fig. 15) confirms the proposed explanations of decreasing brightness and therefore likely decreasing mass accretion rate and midplane temperature. 4.5. About the nature of the two quasi-periodic components in the 'second plateau' With 23 years of coverage of the 'second plateau', in Sec. 3.2 we refined the values of the associated quasi-periods to 1707±70 d and 502±20 d. Precise Schmidt and RC80 telescope data enabled detection of the shorter period in the B and V -band 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 d period all CIs are becoming bluer, but   found the opposite (redder) for the 1707 d quasi-period. The amplitudes related to the longer period decrease, while the amplitudes related to the shorter period increase with increasing effective wavelength (BV R C I C ), respectively (Sec. 3.3). Moreover, our precise Schmidttelescope data revealed that as the time progresses, the amplitudes related with the 1707 d period are becoming smaller in all these filters. No significant evolution of amplitudes related to the 502 d 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 BV R 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 d) 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-2003 is associated with the fall back of dusty material to small radii and the subsequent passage of dust condensations across the line of sight to the inner accretion disk (Sec. 4.5.2 in their paper). However, at the time of this study the periodic behavior was unknown. Enriched with this knowl-edge, we conclude that this event can be uniquely associated with the 1707 d period. We base this result on respective amplitudes obtained in Sec. 3.3, which are decreasing with increasing effective wavelengths of consecutive filters roughly in line with the mean extinction law (Fig. 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 either by dispersion of this dusty region in time, and/or lower dust production rates as the disk wind is the subject of weakening due 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 (Sec. 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.
According to Güver &Özel (2009), N H [g cm −2 ]= 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 Fig. 4b), 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 Fig. 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 d QPO and the absorption components of the wind lines. In Fig. 19ab we present Hα and Hβ P Cygni profiles in function of phase calculated for this QPO (Fig. 4b). The radial velocities (RV) 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 d), with mean amplitude of ∆RV= 5.4±1.6 km s −1 and mean velocity of γ = −81.4 ± 1.2 km s −1 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.
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 doublepeaked broad emission peaks appear with change in accordance to the 1707 d 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 d period correlates with radial velocity changes of the emission component of the lithium line (Kopatskaya et al. 2013). The reversed color index behavior associated with the 502 d 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 localised at the distance of 1 au from the 0.5 M star (Gramajo et al. 2014) to be responsible for the 502 d QPO. Using the model presented in Siwak et al. (2018), but assuming blackbody instead of supergiant spectral intensities 7 , R in = 4.6 R and i=62 • (Sec. 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 Sec. 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 (Fig. 20). The obtained hot spot size and temperature is in principle consistent with that 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 −1 , which is very different from those reported by the authors (from −9 to −19 km s −1 ). We found no correlations between the 502 d 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 of astronomical units (Reggiani et al. 2014;Biller et al. 2014).
The 502 d quasi-periodic phenomenon can also be caused by obscuration of certain disk annuli, i.e. those mostly emitting in the red and near-infrared bands, by dust cloud. In order to reproduce the observed amplitudes, we first calculated magnitudes for an unshaded disk assuming steady accretion, 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 − 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 circular orbit can immediately be excluded. Occultation of the innermost disk annuli as in the 1707 d scenario proposed above does not apply because it would produce a strong signal at U BV wavelengths (in fact only barely detected in BV -bands). The remaining possibilities are a dusty cloud -a remnant of the envelope -rotating on 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 it is not clear the mechanism leading to periodic and continuous brightness variations. Thus, based on the data available so far, we conclude that obscuration scenarios considered above fail to explain the 502 d period and its spectral properties in a manner as self-consistent as a locally heated disk at 1 au.

About the intra-day and weekly variability observed by TESS
While the nature of the quasi-periodic light variations observed from the ground since 1997 seem 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 −1 ). These different relationships could be initially interpreted by combination of extinction, accretion and light-scattering processes. We arrive at this conclusion due 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(Siwak et al. , 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.

CONCLUSIONS
In this paper, we reported on a multi-epoch, multiwavelength study of V1057 Cyg, a classical FU Orionistype object. We arrived at the following conclusions: • The Gaia DR2 distance of V1057 Cyg of 897 pc is significantly larger than previous estimates in the literature, making this object more luminous than previously thought. We constructed multi-epoch SEDs, which we modeled with a simple geometrically thin, optically thick accretion disk model, with A V and MṀ as the two free parameters.
Our results show that the accretion rate reached 1×10 −3 M M yr −1 at the peak of the outburst in 1971 and is still about 1 × 10 −4 M M yr −1 . This makes it the most heavily accreting FUor ever discovered.
• Our long-term photometric monitoring shows the continuation of the second photometric plateau, showing a year to year variability in the optical bands limited to only a few tenths of a magnitude. A period analysis reveals a longer 1707±70 d and a shorter 502±20 d period. Our study detected the shorter period in the BV -filters for the first time. The amplitudes related to the longer period decrease with wavelength, while the ampli-tudes related to the shorter period increase with wavelength. Our data revealed that the amplitudes related to the 1707 d period have decreased during the last two decades. No clear evidence for evolution of the amplitudes related to the 502 d period was observed.
• Based on optical color-magnitude diagrams, we conclude that during the 'second plateau' the color index variations generally follow the extinction path. Due to mutual overlap of the two periodic components having different spectral characteristics, color index variations observed at certain times (especially in 2019) show different relationship, which cannot be explained either by changing extinction along the line of sight, or variable accretion, or combination of both.
• The origin of the 1707 d periodicity might be related to an orbiting dust structure periodically eclipsing the central part of the disk. This may likely arise from a fall back of dusty material from the envelope to small disk radii. For the 502 d period, the most plausible interpretation is a 0.3 × 0.2 au 'bubble' heated to 3500-4000 K, located at 1 au from the star.  In classical T Tauri stars, these lines are jets/outflows tracers, but they are not common in FUors. This is the first detection ever of these lines in a classical FUor. These lines are also variable in time, suggesting that the jet/outflow activity is not constant in V1057 Cyg.
• We obtained a new NIR spectrum of V1057 Cyg in 2020 and compared it with previous spectra from the literature. In this wavelength range we also detected various absorption and emission features. The strength of the CO first overtone bandhead absorption has been decreasing since its first observation, weakest in 2020. This can be interpreted as a sign of decreasing mass accretion rate.
• Our spectroscopic analysis shows that the properties of V1057 Cyg still mostly resemble those of a classical FUor. The photometric monitoring also indicates that it has not yet returned to quiescence, therefore, the FUor outburst of V1057 Cyg is still ongoing.