A multi-epoch, multi-wavelength study of the classical FUor V1515 Cyg approaching quiescence

Historically, FU Orionis-type stars are low-mass, pre-main sequence stars. The members of this class experience powerful accretion outbursts and remain in an enhanced accretion state for decades or centuries. V1515 Cyg, a classical FUor, started brightening in the 1940s and reached its peak brightness in the late 1970s. Following a sudden decrease in brightness it stayed in a minimum state for a few months, then started a brightening for several years. We present results of our ground-based photometric monitoring complemented with optical/NIR spectroscopic monitoring. Our light curves show a long-term fading with strong variability on weekly and monthly time scales. The optical spectra show P Cygni profiles and broad blue-shifted absorption lines, common properties of FUors. However, V1515 Cyg lacks the P Cygni profile in the Ca II 8498 \r{A} line, a part of the Ca infrared triplet (IRT), formed by an outflowing wind, suggesting that the absorbing gas in the wind is optically thin. The newly obtained near-infrared spectrum shows the strengthening of the CO bandhead and the FeH molecular band, indicating that the disk has become cooler since the last spectroscopic observation in 2015. The current luminosity of the accretion disk dropped from the peak value of 138 $L_{\odot}$ to about 45 $L_{\odot}$, suggesting that the long-term fading is also partly caused by the dropping of the accretion rate.


INTRODUCTION
A rather small, but spectacular class of the low-mass pre-main sequence stars are the FU Orionis type objects (briefly FUors). The group was defined by Herbig (1977) based on the properties of the first few outbursting sources observed in the optical, today known as classical FUors. However, these and the other optically visible sources represent only part of a more complex phenomenon that is conspicuous also at infrared wavelengths. FUors are characterized by their enormous inner disks brightness increase, which is caused by the enhanced accretion from the disk onto the protostar. Enhanced accretion is most likely caused by disk instabilities, and this stage can last for several decades, or likely even centuries (Paczynski 1976;Lin & Papaloizou 1985;Kenyon et al. 1988;Bell et al. 1995;Turner et al. 1997;Audard et al. 2014;Kadam et al. 2020). FUors show very similar optical spectra with the properties of F or G supergiants with broad absorption lines, P Cygni profile of Hα, wind components, and strong Li I 6707Å absorption in the optical regime.
This paper is the second in a series started by Szabó et al. (2021), with the aim of characterising past and current evolutionary stages of the so called classical FUors by means of archival and new photometric as well as spectroscopic observations. The target of the present study, V1515 Cyg, was the third discovered FUor. Unlike the first two classical FUors (FU Orionis and V1057 Cyg), V1515 Cyg did not reach the photometric maximum brightness within a year, but this process lasted significantly longer: it started in the 1940's and it reached the peak (13.7 mag in the B-band) at the end of the 1970's (Landolt 1977). In 1980 a significant, yet only temporary fading occurred, which was the sharpest decline yet observed during a FUor outburst (Kolotilov & Petrov 1981;Kolotilov 1983;. V1515 Cyg faded 1.5 mag in the B-band within a few months, but after reaching the minimum at the end of 1980, the source started to brighten again (Clarke et al. 2005). The fading observed in V1515 Cyg after reaching the light curve maximum is commonly observed among the classical FUors: both FU Ori and V1057 Cyg showed similar dips after reaching their maximum (Kolotilov & Kenyon 1997a;. Unlike V1057 Cyg, where the pre-outburst spectrum showed properties of a classical T Tauri-type star (CTTS; Herbig 1977;Herbig et al. 2003), no spectroscopic observations were carried out prior to the brightening of V1515 Cyg. The spectra obtained during the first fading event in the 1980's by Kolotilov (1983) showed absorption lines typical for high-luminosity G2-G5 star and despite the fading, the spectrum of V1515 Cyg has not fundamentally changed since the observations obtained around 1974-1975by Herbig (1977. All features common for FUors like broad Hα line having P Cyg profile, strong Li I 6707Å as well as strong Na D lines, were present. The presence of blueshifted absorption in the Hα and Hβ lines are common signatures for the outflowing wind in these objects. In the case of V1515 Cyg the blueshifted features were interpreted as an expanding envelope by Kolotilov (1983). Using photometric data,  proposed to explain the 1980 fading by a sudden dust grain condensation event in the outflowing wind. They proposed that the slow recovery from the minimum was most likely caused by the gradual dispersal (expansion) of a newly formed dust cloud.
Several years ago, we started a photometric monitoring program of a few classical FUors at the Piszkéstető Mountain Station of Konkoly Observatory in Hungary, including V1515 Cyg. After combination with other public-domain photometric data of this source, we realized that about 15 years ago the source started a longterm fading trend. Given that the last detailed photometric analysis of V1515 Cyg was presented by Clarke et al. (2005), while the last optical and infrared spectral analysis were performed by , Agra-Amboage & Garcia (2014) and Connelley & Reipurth (2018), we decided to increase the cadence of our photometric monitoring and obtain new spectroscopic observations to properly characterise this fading process. Although the FUor phenomenon can occur throughout the YSO evolution (see e.g., Hartmann & Kenyon 1996;Fischer et al. 2022), it is still an open question how an eruption ends, and how the properties of a system undergoing an eruption are evolving. Considering the longterm fading trend, V1515 Cyg could easily be the first FUor going back to the quiescent stage, therefore, it is crucial to gather as much data as possible to better un-derstand an important phenomenon in the evolution of young stars.
In Section 2 we briefly summarize our observations and the data reduction. We analyse the data in Section 3 and we discuss obtained results in Section 4. Lastly, we summarize our findings in Section 5.
2. OBSERVATIONS AND DATA REDUCTION 2.1. Ground-based optical and near-infrared photometry The bulk of photometric observations was obtained at Piszkéstető Observatory between 2005 and 2021 in B, V , R C , I C , g , r , and i filters. Three telescopes, each having somewhat different optical system, were used. Between 2005 October and 2011 May 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). In August 2011 -December 2019 we observed the star with the 60/90 cm Schmidt telescope (pixel scale: 1. 027) and from May 2020 we use the Astro Systeme Austria AZ800 alt-azimuth direct drive 80-cm Ritchey-Chrétien (RC80) telescope (pixel scale 0. 55).
In 2006 we also observed V1515 Cyg at Mt. Maidanak observatory by means of the 100 and the 60 cm Zeiss telescopes, both equipped with the same IMG1001E 1k×1k CCD camera and BV R C I C filters set: the larger reflector was used during 9 nights between 17 -27 June, while the smaller during 35 nights between 4 July -28 August. Due to different focal lengths, the larger telescope offered 6.35, while the smaller 11.74 arcmin field of view, with the pixel scale of 0.372 and 0.688 arcsec pix −1 , respectively. Typically 4 images in B, 2 -3 images in V , and single images in R C I C filters per night were obtained. After standard reduction steps on bias, dark and flat field, we extracted aperture photometry of V1515 Cyg using three comparison stars for the 100 cm, and five comparison stars for the 60 cm telescope. The results from separate images obtained during the same night were averaged. This resulted in typical uncertainty of 0.04 -0.05 mag in B, 0.01 -0.02 mag in V , and 0.01 mag in R C I C filters.
In July -September, 2021, we observed the star at the Mt. Suhora Observatory of the Cracow Pedagogical University (Poland). We used the 60 cm Carl-Zeiss telescope equipped with an Apogee AltaU47 camera, giving 1. 116 pixel scale and 19. 0 × 19. 0 field of view. Johnson BV and Sloan g r i filters were used.
Occasionally we also used other telescopes: on 2006 July 19 and 2012 October 14 we obtained B, V , and R J data with the IAC80 telescope of the Instituto de Astrofísica de Canarias located at Teide Observatory (Canary Islands, Spain) (0. 537 pixel scale). The data reduction for the RCC, Schmidt, RC80, Mt. Suhora and IAC80 observations was carried out in the same way as described in detail by Szabó et al. (2021). We show our photometric results in Table 4 in Appendix A and in Figs. 1 and 2. The typical uncertainty of our measurements is 0.02-0.03 mag in B and 0.01 mag in all other filters.
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. Lastly, we used the NOTCam instrument on the NOT on 2021 August 9. 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 detailed data reduction was carried out the same way as for V1057 Cyg using TCS as described in Szabó et al. (2021). In the same paper we also provided further details about the above instruments. The typical photometric uncertainties are of 0.01-0.03 mag. We present the JHK s photometry in Table 4 in Appendix A. We also observed V1515 Cyg with the 2.56 m Nordic Optical Telescope (NOT) located at the Roque de los Muchachos Observatory (Canary Islands; Plan ID 63-401, PI: Zs. M. Szabó) We used the NOTCam instrument in wide-field imaging mode on 2021 April 25. The instrument consists of a 1024×1024 pixel HgCdTe Rockwell Science Center 'HAWAII'. The field of view was 4 ×4 with a pixel scale of 0. 234.
We supplemented our ground-based optical photometry with publicly available data from the ASAS-SN (Shappee et al. 2014;Kochanek et al. 2017) and ZTF surveys (Masci et al. 2018), which were reduced by dedicated pipelines. We also used the DASCH database (Grindlay et al. 2012) for our work. The first ASAS-SN data were obtained on 2015 February 24, while the latest on 2022 January 4, whereas the ZTF observations range from 2018 March 28, to 2021 November 3. The DASCH data ranges between 1944 and 1989, and seasonal averages were used for plotting purposes in Fig. 1. We also use the 1982 -2003 Mt. Maidanak U BV R data set, already analysed in detail by Clarke et al. (2005).

Space-based optical and infrared photometry
We used the Transiting Exoplanet Survey Satellite (TESS, Ricker et al. 2015), observed V1515 Cyg between 2019 July 18 and September 11 (Sectors 14 and 15), and then between 2021 July 20 and August 20 (Sector 41). The photometry was extracted from the fullfield images using 1.5 pix aperture radius in the same way as described in Szabó et al. (2021). In the final step, the instrumental TESS magnitudes were aligned to our Schmidt-(direct) and RC80-telescope (transformed) measurements in the I C band, and the results are shown in Fig. 3. The seven breaks visible in the light curve in 2019, and the single break in 2021, are caused either by the regular data transfers to the ground or by scattered Earth and Moon light, preventing accurate sky subtraction -these affected data were discarded from the presented plots and from further analyses.
On 2018 August 29, V1515 Cyg was observed with the Stratospheric Observatory for Infrared Astronomy (SOFIA; Young et al. 2012) using the Faint Object in-fraRed CAmera for the SOFIA Telescope (FORCAST; Herter et al. 2013). We obtained mid-infrared imaging using the F056, F077, F088, F111, F242, F315, and F348 filters with effective wavelengths of 5. 6, 7.7, 8.8, 11.1, 24.2, 31.5, and 34.8 µm, respectively (Plan ID 06 0062, PI: J. D. Green). Total on-source integration times were between and 45 s and 87 s depending on the filter. The images were processed using the SOFIA pipeline and retrieved as Level 2 data products from the SOFIA Science Archive as ingested into the IRSA database. We obtained aperture photometry on the images with a series of apertures with increasing radii and selected an aperture for the final photometry where this curve of growth flattened. The resulting fluxes are: F 5.6 = 0.486 ± 0.071 Jy, F 7.7 = 0.502 ± 0.065 Jy, F 8.8 = 0.531 ± 0.069 Jy, F 11.1 = 0.981 ± 0.069 Jy, F 24.2 = 2.42 ± 0.13 Jy, F 31.5 = 2.06 ± 0.16 Jy, F 34.8 = 2.21 ± 0.18 Jy. We also used data from the Wide-field Infrared Survey Explorer (WISE, Wright et al. 2010). In Tab. 5 in Appendix A we summarize the saturation corrected WISE data for V1515 Cyg shown in Fig. 1.

Spectroscopy
We obtained a new optical spectrum of V1515 Cyg with the high-resolution FIbre-fed Echelle Spectrograph (FIES) instrument on the NOT on 2021 May 29. We used a fibre with an entrance aperture of 2. 5 which provided a spectral resolution R=25 000 between 370 -900 nm. The final spectrum is composed of three individual spectra, each with 2400 s exposure time. We reduced the spectra using the FIEStool software. We also observed V1515 Cyg 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). BOES provides R=30 000 (using a 300 µm fiber) in the wavelength range ∼4000-9000Å. We binned the spectra with 2×2 pixels to increase the S/N. We took four spectra of V1515 Cyg between 2015 and 2018, for the exact dates, see the spec- troscopic log of observations in Tab. 1. We reduced these spectra in a standard way within IRAF: after bias and flatfield corrections, the ThAr lamp spectrum was used for wavelength calibration, and continuum fitting was performed using the continuum task. Finally, heliocentric velocity correction was applied by the rvcorrect task and the published radial velocity of V1515 Cyg (−15 km s −1 ; Hartmann et al. 2004). Unlike the usual observing method (and our 2021 NOT spectrum) with several exposures, the BOES spectra were taken with one, long exposure. This made the spectrum sensitive for cosmic rays with no means to remove them using the usual procedures. As no telluric standard stars were observed either for FIES or 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.
On 2021 August 9, we used the NOTCam on the NOT to obtain a new near-infrared spectrum of V1515 Cyg and Iot Cyg (telluric standard, A5 V) in the JHK s bands (1.25-2.2 µm). We used the low-resolution camera mode (R=2500) with ABBA dither positions, and exposure times ranged from 60 to 100 seconds (Tab. 1). For each image, flat-fielding, bad pixel removal, sky subtraction, aperture tracing, and wavelength calibration steps were performed within IRAF. A Xenon lamp spectrum was used for the wavelength calibration. The Hydrogen absorption lines in Iot Cyg were removed by Gaussian fitting and then the spectrum of V1515 Cyg was divided by the normalized spectrum of Iot Cyg for telluric correction. Flux calibration was performed by applying the accretion disk model (Sec. 4.2). It is based on the NOT JHK s photometry taken on the same night, but the photometry was taken in broadband filters, hence gives fluxes only for the effective wavelengths of these filters. We therefore use the SED model (desribed in Sec. 3.7) fitted to the photometric points to interpo-late the fluxes for any given wavelength needed for the flux calibration of the spectrum.

Light Curves
In order to study the brightness variations of V1515 Cyg in the long-term, we constructed light curves with data gathered from the literature (Herbig 1977;Landolt 1977;Kolotilov 1983;Herbst & Shevchenko 1999;Clarke et al. 2005;Grindlay et al. 2012) and we complemented the archival data with results from our monitoring. We present the long-term photometric light curves of V1515 Cyg in Fig. 1. In order to stay consistent between numerous data sets and to simplify further analyses, we transformed our RC80, Mt. Suhora, and ZTF Sloan r i magnitudes into Johnson-Cousins R C I C system by constructing the SED for each day observations and interpolate for the effective wavelengths of the R C and I C filters using a spline function in the log-log space. The resulting interpolated fluxes were transformed to magnitudes with the zero points of the R C I C filters. We stress that this approach gave more consistent results than the transformation formulas between Johnson and Sloan systems given by Jordi et al. (2006). Any remaining offsets between overlapping data sets in optical bands were then removed by constant shifts, as listed in Tab. 6 in the Appendix. We give a brief summary of the different data sets used for Fig. 1 in Tab. 2.
As mentioned in Section 1, the outburst of V1515 Cyg was unlike those any other classical FUors: it brightened for ∼30 years, before reaching the peak brightness in the late 1970's (Landolt 1975(Landolt , 1977, then in 1980 it started a quick fading event, reaching the minimum brightness, and re-brighten again for several years (Kolotilov & Petrov 1981;Kolotilov 1983;Clarke et al. 2005). Between 1984-2004 the light curve showed a plateau shape (Clarke et al. 2005). However, a significant long-term fading of V1515 Cyg started around 2005-2006, with a sudden drop in the brightness (∼1 mag in the BV R filters), a trend that is currently still ongoing with the mean rate of 0.085 mag yr −1 in the V -band. We note, that both the plateau and the longterm fading are visible in the optical BV R, the nearinfrared JHK, and the WISE light curves also confirm the fading trend at 3.4 and 4.6 µm.
The small amplitude light changes were thoroughly analysed by Clarke et al. (2005) between the 1980's and early 2000's. However, our careful look into the same Mt. Maidanak data set revealed a feature unnoticed earlier: between 1984 -2003 the small-scale variations visible on top of the major light curve plateau predomi- 1974 -1991 JHK Molinari et al. (1993) 1975, 1976V Herbig (1977) 1976 U BV Landolt (1977Landolt ( ) 1978Landolt ( -1983 U BV Kolotilov (1983-1995V Herbst & Shevchenko (1999-2005 U BV R Clarke et al. (2005) 1982V Kenyon et al. (1991) 1985JHK Kenyon et al. (1991) 1985-1998K Kolotilov (1990  nantly manifest themselves in the form of peaks. At first glance, they resemble accretion bursts typical for magnetospheric accretion in CTTS, however, the U BV R light curves show very similar variability amplitudes. In combination with the characteristic time scales of these brightness fluctuations (∼ 10 − 30 d) this strongly suggests that all these variations most likely originated in the inner protoplanetary disk. The INTEGRAL V -band data are sparse, thus the peak-like events could not be seen clearly. Instead, these data revealed at least one well-defined dip at JD ≈ 2454700 (2008 August) -perhaps a signature of gradual disk envelope cavity closure by increased dust condensation as the outburst is weakening. The dips are in fact more numerous in the data obtained in 2010s and 2020s.

Frequency Analysis
In order to search for periodic and quasi-periodic oscillations (QPOs) in the ground-based data, we utilised the V -and g-band light curves, providing both the best temporal coverage and showing variability present in other optical bands. In the first step, all outlier points were removed. As the ASAS-SN survey gradually switched to the g-band, we re-scaled the amplitude of these variations to match those observed in the V -band by the multiplicative factor of 0.95. We determined this value during the time period, when simultaneous V -and gband data were obtained, then, we aligned the scaled g-band light curve to the V -band by a constant shift.
In the first step we searched for long-term quasiperiodic variations in the 1982-2021 dataset superim-  1942 1946 1950 1954 1958 1962 1966 1970 1974 1978 1982 1986 1990 1994 1998 2002 2006  posed on the major light curve plateau. For this purpose we subtracted the major long-term trend from the light curve by fitting a 2nd order polynomial. Although the Lomb-Scargle (Zechmeister & Kürster 2009) diagram computed from the entire light curve reveals a forest of peaks, none of them represent stable period nor quasi-period. For this reason, in the second approach we decided to split the light curve based on the dominant variability morphology and analyse each part separately, i.e. the first part obtained before 2005 (JD = 2453500), dominated by peaks, and the second part obtained after, which is dominated by dips. No significant quasiperiodicity was found, except for the three consecutive dip-like 350-370 d variations visible directly in the light curve between JD = 2458000 − 2459000 (2017-2020).
In the third step we decided to analyse each wellsampled observing season separately. Prior to this analysis, any trends visible during certain seasons were removed by using low order (1-2) polynomials to enhance the shorter changes. We obtained that in the 1987 data set there is a significant 13.91 d QPO with false alarm probability of 4×10 −5 . This is similar to the 13.89 d period found by Clarke et al. (2005) in the same data. In addition, we found a significant 13.06 d peak with false alarm probability of 10 −3 in the 2003 data set. During the other seasons, the light curve usually exhibited only two consecutive oscillations of a similar pattern, quickly disappearing and giving way to light changes with different patterns.
In order to study the small-scale variability unavailable from the ground, we utilised the TESS light curves obtained in 2019 and 2021. We performed frequency f analysis using the procedure of Rucinski et al. (2008), where the amplitude a f errors are estimated by means of the bootstrap sampling technique. Prior to the analysis, the general downward trends in the brightness visible during both in 2019 and 2021 ( Fig. 3a,b) were removed by linear fits. Subsequently, the residual magnitudes were transformed to normalized flux units. We found that in both seasons the amplitude-frequency spectrum exhibits a Brownian random-walk nature (Fig. 4a,b), in which the amplitudes scale with the frequency as a f ∼ f −1 (Press 1978). Only in the 2021 data set there is a single peak at f = 0.439 c d −1 , which appears to represent a significant 2.28 d QPO.

Wavelet Analysis
It is known that Fourier analysis cannot trace temporal changes nor localize in time finite oscillatory packages such as those present in our data. As it was for   the first time proven by Rucinski et al. (2008Rucinski et al. ( , 2010, the wavelet analysis can be successfully applied to extract such information from the evenly-sampled space-based light curves of vigorously accreting CTTS and Herbig Ae stars. In particular, the Morlet wavelet allows us to separate individual oscillatory packages and directly characterise their temporal changes. Although the variability observed in the above mentioned classes of young stars is caused by the changing visibility, number, and the lifetime of the numerous hot spots resulting from the (most often) unstable or moderately stable magnetospheric accretion (Blinova et al. 2016), or in some particular cases the pulsed accretion (Tofflemire et al. 2017;Kospál et al. 2018;Tofflemire et al. 2019;Fiorellino et al. 2022), it was already theoretically proven that the hot spot behavior actually reflects the physical conditions in the inner disks of CTTS (see Blinova et al. 2016 and references within). In the case of the FUors, where the disk radiation is dominating over the stellar one, and the magnetosphere is likely totally crushed by the enhanced disk-plasma pressure so that hot spots typical for CTTS cannot form, effects brought by these inner disk instabilities can be directly observed in the visual light only in the objects with residual traces of the envelope (Siwak et al. 2013;Green et al. 2013a;Baek et al. 2015;Siwak et al. 2018Siwak et al. , 2020. We obtained that the Morlet-wavelet spectrum of the 2019 TESS data does reveal only a single, weak signature of QPO drifting from ∼6 to ∼5.5 d between BT JD − 2458000 = 690 − 720 (Fig. 4c). In accord with the visual examination of the light curve, we note random changes, usually showing only one cycle or two self-similar oscillatory cycles. The shortest variations of this kind have a period of 1.5 d, but most of them are longer than 2 d. We have to stress that this conclusion is unfortunately disturbed by the seven interruptions in the TESS data acquisition, strongly affecting the short periodic part of the spectrum, as indicated by the solid white lines in Fig. 4c,d. Moreover, there is a non-zero chance that any of these affected regions may propagate towards the longer periods inside the open cone(s), which can be imaginatively extended as indicated by the white continuous lines.
The shortest significant variations in the 2021 TESS wavelet spectrum are of 0.8-0.9 d and are mostly localised in the middle of the data set. The spectrum does also reveal that the 2.28 d peak obtained in Sec. 3.2 is unstable, and is seen only between BT JD − 2459000 ≈ 425 − 430, i.e. it lasted for barely two oscillatory cycles (Fig. 4d), it re-appeared around BT JD − 2459000 ≈ 441, but it is not well defined in the light curve itself. However, together with the former case it can produce an excess of power in the frequency spectrum, which leads to the spurious detection.
Similarly to Szabó et al. (2021), we also applied 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). In general, the analysis did not show any oscillatory packages other than those directly inferred from the best-sampled Mt. Maidanak light curves themselves, i.e. typically consisting of two self-similar oscillations (see in Sec. 3.2). We show only the results for the 1987 and 2003 seasons, although instead of the WWZ, we decided to present the ordinary Morlet wavelet spectra. We stress that the spectra were computed from the original data interpolated into evenly-spaced time grid with 1 day step, which is equal to the median spacing of the Mt. Maidanak original data (Fig. 4e,f). The signatures of the 13.1 d QPO are evidently best defined for the 2003 data set, but we stress that the detection efficiency is a function of the data sampling, as shown in the associated bottom panels. For that reason, the 13.9 d QPO from the 1987 season is less securely defined, as the maximum amplitude appears before the first dashed line, denoting edge effects, and then the amplitude and the period appear to change.

The Curious Case of the 2021 Phenomenon
The data coverage from ground-based telescopes was exceptional in 2021. The early 2021 season data showed brightening trend by 0.2 -0.4 mag, which stabilised for about 1 -1.5 month, then (in June 2021) it faded quickly and stayed nearly-constant for about one month, exactly during the TESS observations in July -August. In early September, the brightness started to rise and stabilized on the level expected from the long-term fading trend.
Taking advantage of the excellent coverage, we analysed this phenomenon in detail as follows. First we detrended the long-term fading trend by linear fit to the high brightness level, as indicated in the upper panel of Fig. 5. One can notice that brightness evolution in different filters occurred in slightly different way. This is better emphasized on color-magnitude diagrams, as shown in the bottom panels of Fig. 5. It is evident that the brightness drop followed the extinction path. Then, in the light curve minimum the colors became bluer, which suggests scattering by small grains, like in UXor stars (see e.g., Ábrahám et al. 2018). Finally during the rising branch, all color indices remained almost unchanged, i.e. the brightness went to the upper level almost on a grey path.
We performed the same analysis for the 2019 and 2020 observing seasons, but conclusive results could not be     reached due to the coarser sampling. Therefore, we only focused on the 2021 data set, since this is the most useful one to draw conclusions thanks to the good data sampling. We discuss in detail the possible physical interpretation of the dimming events of V1515 Cyg in comparison with the other classical FUors and the classical T Tauri star (CTTS) RW Aur in Sec. 4.1.

Spectroscopy
We monitored V1515 Cyg from 2012 to 2021 with optical and near-infrared spectrographs and obtained six spectra. We detected P Cygni profiles, broad blueshifted absorption lines, and metallic lines in absorption. For the identification of the lines the NIST Atomic Spectra Database 1 was used.

Optical Spectroscopy
Classical FUors share several common spectroscopic features in their optical spectra, such as the P Cygni profile of Hα, strongly blue-shifted absorption lines like the Li I 6707Å absorption, and the fact that their spectral type is wavelength-dependent (Hartmann & Kenyon 1996;Audard et al. 2014). Some of these spectroscopic characteristics are also seen in our observations, presented in the upper panel of Fig. 6. Apart from 2018 December spectrum, all the other spectra show remarkably similar profiles within the measurement uncertainties. The 2018 December spectrum in Fig. 6 is different in the way that the absorption profile was not as deep as in all the other epochs. We note that for plotting purposes in this section, we used our NOT/FIES spectrum obtained in 2021 May.
We observed several P Cygni profiles, specifically in Hα, Hβ, and two lines of the Ca infrared triplet (IRT), Ca II 8542, and Ca II 8662Å lines, but interestingly not in the Ca II 8498Å line. Fig. 7 shows observed P Cygni profiles in the 2021 spectrum. In the case of the Hα line, the previously published optical spectrum of V1515 Cyg presented by Herbig (1977) and  showed a pure absorption profile, but we observed Hα with a clear P Cygni profile in all our observations taken between 2015 and 2021 (middle panel in Fig. 6), including a clear emission component. Hα and Hβ show higher velocity of both blue-shifted absorption (∼ −300 km s −1 ) and red-shifted emission (∼100 km s −1 ) components compared to the Ca II lines (∼ −200 km s −1 and ∼50 km s −1 ).
Wind signatures in the form of broad blue-shifted absorption lines were also observed in the Na I D 5889 1 https://physics.nist.gov/PhysRefData/ASD/lines form.html and 5895Å, Mg I 5167 and 5184Å, and K I doublet lines (7664 and 7698Å) (Fig. 8). The highest velocity component extended up to about 250 km s −1 , similar to the blue-shifted absorption component of Hα and Hβ. Moreover, some lines show multiple velocity components, namely the Na I D 5889, Hα 6562, K I 7664, and Ca II 8542Å lines (Fig. 9). These lines share at least three velocity components among around −141, −65, −44, and −16 km s −1 , implying that the lines might form in the same region. Here we note that jet tracer lines, such as [S II], [N II], and [O III] were not detected in our multi-epoch spectra as opposed to V1057 Cyg (Szabó et al. 2021).
We also observed numerous metallic absorption lines during our observations. These lines include various Fe I lines between 5446 and 7511Å, Ca I 6122-6471Å, Ti I 5965, 7344Å, Ni I 6108, 6767Å, and the Li I 6707Å line (see some of these lines in Fig. 16 in Section 4.3). In order to estimate the rotational velocity (vsini), we convolved the stellar spectrum with a Gaussian profile to account for the rotational broadening of the lines and compared their spectra to V1515 Cyg. For this analysis, we adopted the standard stellar spectra from F8 to G8 supergiants used in Park et al. (2020), which are observed with the same observing setup as V1515 Cyg. First, we averaged four spectral lines (Fe I 6141.7, Fe I 6411.6, Ca I 6439.1, and Fe I 7511.0Å) by normalizing each peak intensity for V1515 Cyg and standard stars. Then, the stellar spectrum was convolved with a Gaussian to account for various rotational broadening from 0 to 30 km s −1 with a 1 km s −1 step. The best fit was found by χ 2 minimization, and the best fit rotational velocity is 15 km s −1 found using G5 Ib-II type stellar template (Fig. 10). The obtained rotational velocity is similar to the previously found ∼20 km s −1 (Hartmann & Kenyon 1985) and consistent with low disk inclination (Gramajo et al. 2014;Hillenbrand et al. 2019;Milliner et al. 2019). Furthermore, the best fit spectral type (G5 Ib-II) is also consistent with G supergiants ).

Near-infrared Spectroscopy
Our recent NIR spectrum taken in 2021 shows several absorption lines: Paβ, Al I, Mg I, FeH molecular bands, and strong CO absorption bandhead features are present, similarly like in the previous spectra of V1515 Cyg taken by Connelley & Reipurth (2018). In Fig. 11 we compare their 2015 June spectrum and our spectrum from 2021 August: in the first panel we show the full JHK spectra, the next three panels show the J, H, and K-band spectrum respectively. Although the NIR photometric observations are insufficient since 2012, V1515 Cyg has almost the same brightness in  2012 and 2020 observations. There is no significant flux change of the JHK s spectrum in 2015 and 2020. However, our observation shows interesting changes in some spectral features: the FeH (> 1.62 µm) and CO (> 2.29 µm) molecular bands became stronger. Both the CO rovibrational and the FeH molecular bands are stronger at lower temperatures (Wallace & Hinkle 2001;Bieging et al. 2002;Hargreaves et al. 2010), implying that V1515 Cyg became cooler than observed in the previous observation in 2015 by Connelley & Reipurth (2018).
The Paβ equivalent widths (EWs) became stronger in our recent NOTCam (1.33 ± 0.53Å) observation compared to the IRTF (1.03 ± 0.09Å) observations of Connelley & Reipurth (2018). However the uncertainties are higher for the new observations, therefore, this change may not be significant. The usual scaling relations between the fluxes of the hydrogen emission lines were established for T Tauri stars (see e.g., Alcalá et al. 2017). These relations work because the accretion in all of these  stars happen according to the magnetospheric accretion scenario (see e.g., Hartmann et al. 2016). This scenario cannot be applied to FUors, as indicated by their lack of emission lines. The absorption lines are associated with an active accretion disk, and their formation is explained by temperature inversion (hot midplane, colder disk atmosphere) (see e.g., Hartmann & Kenyon 1996). In V1515 Cyg we see the Paβ line in absorption, as most of the lines observed in the spectrum of FUors in general. Therefore, the hydrogen lines cannot be used for the mass accretion calculation unlike in T Tauri stars.

Accretion Disk Modeling
In this section we model the inner part of the system with a steady, optically thick and geometrically thin viscous accretion disk. We already applied this method for a number of young eruptive stars (Kóspál et al. 2016Ábrahám et al. 2018;Kun et al. 2019; Szegedi- Averaged metallic absorption profiles of V1515 Cyg in 2021 May (black) and standard star (red) convolved with a rotational velocity of 15 km s −1 .
Elek et al. 2020; Szabó et al. 2021) to separate the effects of variable accretion, extinction and study their long term evolution quantitatively. The mass-accretion rate in our model 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 . The inner radius of the disk mainly influences the emission in the optical bands. In practice, we usually cannot discriminate between a small stellar radius with higher line-of-sight extinction and a larger radius with lower extinction, using only broad-band optical photometry. In the subsequent modeling we prescribe that the aver- Figure 11. The NIR J, H, and K spectrum of V1515 Cyg. Gray and red colors present the spectra observed with IRTF (Connelley & Reipurth 2018) and NOTCam, respectively. First panel shows the full JHK spectrum, second the J, third the H, then the K-band spectrum.
age of our A V values obtained between 2012 and 2021 must comply with the A V =3.5±0.4 mag proposed by Connelley & Reipurth (2018) based on an infrared spectrum taken in 2015. This constraint sets R in = 1.4 R . The outer radius of the accretion disk can be estimated by fitting the mid-IR emission of the disk. We fixed R out = 0.45 au, which matches the L-and M -band observations of V1515 Cyg from 1989 . The inclination of the disk is also an important input parameter of our accretion disk model. Most studies agree that V1515 Cyg is seen close to poleon (see e.g., Gramajo et al. 2014;Milliner et al. 2019). In the following we adopt i=10 • .
Our models have only two free parameters: the product of the accretion rate and the stellar mass MṀ , as well as the line-of-sight extinction A V . Disk SEDs were calculated for a large range of MṀ values, and at each step the fluxes were reddened using a grid of A V values 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 quasi-simultaneous flux values between 0.4 and 2.5 µm. 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 reproduced the measurements well, with low reduced χ 2 values. The resulting temporal evolution of the accretion rate and extinction values, together with the optical V -band and the near-infrared K-band light curve, are plotted in Fig. 12, where the K-band light curve is shifted by +4.2 mag for plotting purposes. We discuss these results further in the next subsection and in the Sec. 4.  (Green et al. 2013b) and the WISE (Wright et al. 2010), and we also used data from the airborne SOFIA (Young et al. 2012) missions. The data points from the Herschel and the AllWISE catalogs were obtained within a year, thus we combined the two data sets into one SED. In the 2 Data were taken from "A Spitzer Legacy Survey of the Cygnus-X Complex" catalogue, available at the IRSA database https: //irsa.ipac.caltech.edu/data/SPITZER/Cygnus-X/ far-infrared domain (λ > 25 µm) the IRAS and ISO measurements were apparently contaminated by some extended emission in their large beams, thus were discarded from the plot. Contamination was confirmed by archival 70 µm and 160 µm images by the Herschel Space Telescope, which show that V1515 Cyg is located at the edge of a dust cloud, on a structured far-infrared background.

Spectral Energy Distribution
There are several properties of the SED of V1515 Cyg, which are discussed below in Sec. 4.2 4. DISCUSSION

Dimming events in V1515 Cyg
In Sec. 3.4, we described the interesting dimming event of V1515 Cyg observed in 2021. In the following we investigate the possible physical cause/interpretation of such dimming events.  and Kolotilov & Kenyon (1997b) assumed that the occultation events were due to dust formation in an expanding shell in the wind. However, Clarke et al. (2005) considered that in the case of V1515 Cyg the occultation events are interactions between the wind and the surrounding infalling envelope. They proposed that, for instance, V1057 Cyg has a strong wind and was able to clear the circumstellar material in response to the outburst. The wind in V1515 Cyg has always been weaker, therefore it has never been able to completely clear the line of sight of dusty material, which can explain the variations since the beginning of the outburst. The dust lifted by the wind could cause this variability. In contrast, FU Ori was able to clear a large cavity, hence its light curve is free from this kind of variability (Clarke et al. 2005). This could also explain why there is no correlation between the photometric and the spectroscopic variation suggested by Clarke et al. (2005). We tried to draw an analogy with the T Tauri star, RW Aur, where multiple, similar dimming events were observed with similarly good data sampling. Koutoulaki et al. (2019) suggested that the deep minima are most likely caused by a layer of dust obscuring the inner 0.05 -0.1 au of the system, while the accretion rate did not change significantly (see also other studies, e.g., Petrov et al. 2015;Facchini et al. 2016, and references therein). However, there are many obvious differences between the two systems. Firstly, that although V1515 Cyg is approaching the quiescence phase, it is still in outburst, hence exhibiting a more enhanced accretion rate compared to RW Aur. Therefore, the optical photometry of V1515 Cyg is still dominated by the accretion disk (and the star is negligible), as is usual in FUors in outburst. Moreover, in V1515 Cyg we see the continuous decrease of the accretion rate despite the dimming events. Another difference is the that RW Aur is a binary, while V1515 Cyg seems to be a single star.
In order to understand the observed dimming event in 2021, additional, similarly well sampled multi-filter light curves, ideally complemented by spectroscopic monitor-ing, are needed to look for similarities or differences between these two kinds of dimming events. Fig. 12 shows the temporal evolution of the accretion rate and the extinction together with the optical V -band and the near-infrared K-band light curves of V1515 Cyg. The long-term evolution of the accretion rate until 2010 is approximately consistent with an exponential decay, with an e-folding time of 20600 days (∼56 years). The dashed line in Fig. 12 was fitted to the V -band data after 2003 in order to determine the e-folding time of the long-term optical fading, corresponding to approximately 12 years.

Accretion Disk Modeling and SED
The peak and current MṀ values are 1.29×10 −5 and 4.17×10 −6 M 2 yr −1 , respectively. Assuming a stellar mass of 0.3 M , the real accretion is higher by a factor of about 3. The A V curve is rather constant: a weak gradual decrease is hinted before 2000, which might be responsible for the brightening of the system during this period, in spite of the decreasing accretion rates. The current luminosity of the accretion disk is about 45 L , dropped from a peak value of 138 L .
The SED of the source has several properties which we discuss further below. The SED is dominated by a peak at optical-near-NIR wavelengths. The clear excess above a blackbody-like extrapolation of this peak to infrared wavelengths implies the presence of a significant amount of circumstellar matter (the total circumstellar mass, based on IRAM observations, is 0.04 M (Fehér et al. 2017)). The spectral range around 10 µm and 18 µm exhibits broad emission features, which is attributed to submicrometer size silicate grains (see also e.g., Green et al. 2006). The strength of the short wavelength peak monotonically decreased with time, reflecting the evolution of the hot inner accretion disk modeled in Sec. 3.6. A temporal descrease of the mid-infrared flux is also suggested, although the data coverage is less complete than in the optical.
At far-infrared and submillimeter wavelengths the Herschel observations may represent the emission of the circumstellar disk, with no clear sign of a massive cold circumstellar envelope. There is no information on the temporal evolution of the system during the outburst at these wavelengths. The apparent consistency between the 450 µm JCMT/SCUBA (Sandell & Weintraub 2001a) single dish measurement from 1998 (Sandell & Weintraub 2001b) and the 500 µm Herschel/SPIRE observations from 2011 suggests that the flux change at these long wavelengths is not significant.
Due to the dispersion of the data points, the color evolution can be consistent with evolution along either con-stant accretion or constant extinction path, in general the color evolution is consistent with the evolution along constant accretion or constant extinction path, obtained from our accretion disk model in Sec. 3.6 and Sec. 4.2. Considering the small variability amplitude and the photometric uncertainty, it is not possible to decide which is the dominant mechanism. As the global long-term color evolution is accounted by the modeling, we do not present here the global color-magnitude diagrams but refer the reader to Sec. 3.6 and Sec. 4.2 for the general results.

Comparison between V1515 Cyg and V1057 Cyg
In this section we give a comparison between V1515 Cyg and V1057 Cyg. First we compare their light curves, then their spectroscopic properties using our previous study by Szabó et al. (2021) with the caveat that the spectra for the two sources were taken in different stages of their outburst evolution. Both sources are considered classical FUors (Herbig 1977;Herbig et al. 2003;Connelley & Reipurth 2018), and we were able to observe similar properties in the spectra of both targets, comparing them in the following.
Similarities. A similar aspect of both sources can be found in their light curves: both objects underwent a sharp dip event after reaching the peak brightness. In the case of V1515 Cyg this happened right after reaching the light maximum in 1980, recovering and then showing a plateau stage until the early 2000's ( Fig. 1). However, the light curve of V1057 Cyg is different, this source also experienced a similar, sharp dip, around 1995 after an exponential fading trend. Similarly, just like in V1515 Cyg, V1057 Cyg also recovered from this event and showed a plateau stage afterwards (called 'second plateau ' Szabó et al. 2021). Interestingly, FU Ori also showed a similar event, followed by a long-term plateau stage (see e.g., Hartmann et al. 2016).
In order to compare two FUors homogeneously, we used the latest observation of our NOT/FIES optical spectra for both. Our targets show four Balmer lines from Hα to Hδ in our observations with broad blueshifted absorption profiles. Among these four lines, relatively clear (higher S/N) Hα and Hβ are presented in Fig. 14. Hα and Hβ show P Cygni profiles with broad blue-shifted absorption component and red-shifted emission profile in both targets. The two FUors show P Cygni profiles of the Hα, Ca II 8542, and 8662Å lines as shown in Fig. 14. In addition, both FUors show blueshifted absorption line profiles formed by an outflowing wind in several neutral and singly ionized metallic lines (Fig. 15). However, V1057 Cyg shows higher velocity components of shell features than V1515 Cyg (e.g., Ba II 4934Å and Fig. 11 in Szabó et al. 2021). Differences.
Even though we see similarities in their light curves, perhaps the most striking difference between the two classical FUors is the brightening phase presented in their light curves. The outburst of V1057 Cyg was an abrupt one, happening within less than a year (see Fig. 1 in Szabó et al. 2021). The outburst of V1515 Cyg was different from the beginning compared to the other classical FUors, taking about four decades to reach the maximum brightness. This difference suggests that overall the outburst mechanisms are different.
Interestingly, in the optical spectrum of V1515 Cyg there is a lack of P Cygni profile in the Ca II 8498Å line compared to V1057 Cyg. In our previous study on V1057 Cyg, we found that all three lines of the Ca II IRT showed clear P Cygni profiles, whereas in the case of V1515 Cyg, we found P Cygni profiles only for the Ca II 8542 and Ca II 8662Å lines. The strength/depth of the absorption is also weaker than those observed in V1057 Cyg (Herbig et al. 2003;Szabó et al. 2021). Several works investigated the relation between the Ca II IRT lines in T Taur stars and found that if the lines are optically thin the ratio of the three lines is 1:9:5, and in the optically thick case it is 1:1:1 (see e.g., Herbig & Soderblom 1980;Hamann & Persson 1989Azevedo et al. 2006). This means that the absorption in the 8498Å line would be 9 times weaker than in the 8542Å line. The fact, that the absorption component is not detected (i.e., remains below our detection threshold) in the 8498Å line suggests that the gas in the outflowing wind causing the absorption component is probably optically thin.
Since the strength of the blue-shifted absorption component in a P Cygni profile is formed by an outflowing wind in the system, in the case of V1515 Cyg, the weaker strength could be an indication of the weaker outflowing wind in the structure. Compared to V1057 Cyg where the Ca IRT lines clearly varied with time, in V1515 Cyg, we did not observe strong temporal variation of the lines either. The strength of the P Cygni profiles varies over time, however, they stay within the same velocity range. In addition, no shell lines nor forbidden emission lines were observed in V1515 Cyg as opposed to V1057 Cyg (Szabó et al. 2021). Fig. 16 (Szabó et al. 2021). The absence of these lines further supports the idea of a weaker, uncollimated optically thin wind in this system. Fig. 17 shows the comparison of the NIR spectra between V1515 Cyg (dark gray and blue) and V1057 Cyg (light gray and red). Overall, the flux of V1515 Cyg is lower than the one of V1057 Cyg, and the flux change between 2015 and 2020/21 is also smaller in V1515 Cyg. Most of the observed lines are the same as those of V1057 Cyg in the NIR spectra. To compare spectral changes in V1515 Cyg itself and the two FUors, we measured the EW of the CO overtone bandhead features (2.293 -2.317 µm) with the same method described in detail in Szabó et al. (2021). Fig. 18 shows the CO bandhead of V1515 Cyg (upper spectra) and V1057 Cyg (lower spectra), and the measured EW is listed in Table 3. In the case of V1057 Cyg, the CO 2-0 bandhead feature became weaker with the decreasing brightness, and it is interpreted as the decreasing mass accretion rate (Szabó et al. 2021). On the other hand, in the case of V1515 Cyg, the bandhead features almost remain the same strength (Fig. 18), similarly to the monitoring results of the optical spectra (Fig. 6). In our observation, the EW of the CO increases, which is mostly affected by the lower temperature rovibrational transitions (∼2.31 µm) in Fig. 18. This result is consistent with the presence of the stronger FeH in our observation (Fig. 11). We estimated the gas temperature by comparing the FeH feature at 1.62 µm in V1515 Cyg and those of standard stars (Rayner et al. 2009) between M4 and L5 spectral types (Fig. 19). Both the target and the stellar spectrum were normalized. We found the best fit by eye due to the blended molecular nature of this feature. The previous spectrum taken in 2015 by Connelley & Reipurth (2018) was well-matched with an M6 -M7 dwarf spectrum (T eff ∼2700 K (GJ 406); Kuznetsov et al. 2019), while the 2021 spectrum is well-matched with a later type dwarf of M9 -M9.5 (T eff ∼2500 K (LHS 2924); Gagné et al. 2015). Together with the decreasing temperature from the FeH, and the increased EW of the CO bandhead in the 2021 observation (see Sec 3.5.2), we conclude that the disk of V1515 Cyg became cooler since 2015, further suggesting that the source is approaching a quiescence phase.
An accurate comparison is still difficult to make between the two sources. The two system have very different light curves, indicating that overall different energetics were present in the early stages of the outburst (see Fig. 1 in Szabó et al. 2021). While a long term fading trend is present in both cases, the decline is currently much faster and drastic in V1515 Cyg. As shown, the current mass accretion rate of V1515 Cyg a Spectrum from Connelley & Reipurth (2018) is much lower (4.17×10 −6 M 2 yr −1 ) than in V1057 Cyg (10 −4 M 2 yr −1 ), which can justify the lower mass ejection rate and the optically thin outflow in the former system. Overall, the weaker wind in V1515 Cyg further suggests that the outburst mechanisms are quite different in the two sources. It is important to mention that different geometries can also play a role in differences between systems, however the inclination of the two sources are similar as both systems are seen close to face-on (see e.g., Gramajo et al. 2014;Milliner et al. 2019).

CONCLUSIONS
In this paper we performed a multi-epoch, multiwavelength study of the classical FU Orionis-type star V1515 Cyg and we arrived to the following conclusions: • Our long-term photometric monitoring of V1515 Cyg supplemented with data obtained at Mt. Maidanak and data from publicly available archives such as INTEGRAL, ASAS-SN, and ZTF show a long-term fading trend which started around 2006. The fading is not monotonous: after the maximum brightness , the inner disk apparently experienced several obscurationlike events. The long-term evolution of the source's optical brightness after 2003 is approximately consistent with an exponential decay, with an e-folding time of approximately 12 years.
• The Mt. Maidanak photometry enabled us to investigate the light curve changes occurring in the timescale of 5-100 d during the maximum. Comparisons of these multi colour light curves suggest that the majority of these light curve changes did originate from the inner disk. We confirmed the 13.9 d QPO found in the 1987 data set by Clarke et al. (2005)       • Although TESS enabled investigation of shortperiodic, down to ∼20 min oscillations, the 2019 and 2021 light curves do not show rapidly occurring light variations on such time scales. The frequency spectra exhibit the random-walk character, being of the same type as the one found in our previous study of V1057 Cyg (Szabó et al. 2021). It is plausible that the dimming event studied in great detail during the 2021 season could be similar to what was observed in the CTTS RW Aur, however, there are many differences between the two system. Additional, similarly well sampled seasons are required to investigate and compare this phenomena further in the future.
• Our accretion disk modelling reveals that the current luminosity of the accretion disk dropped from the peak value of 138 L to about 45 L . The peak and current MṀ values are 1.29×10 −5 and 4.17×10 −6 M 2 yr −1 , respectively. This indicates that the long-term fading is also partly caused by the dropping of the accretion rate. The optical and near-infrared peak present in the SED indicates the presence of significant circumstellar matter surrounding the source. The long-term evolution of the accretion rate until 2010 is approximately consistent with an exponential decay, with an e-folding time of ∼56 years.
• We performed optical spectroscopic monitoring of the source between 2015-2021 and detected several absorption lines, which show small variations over time. Several high-velocity wind features were also observed in the form of absorption components of P Cygni profiles which similarly show small variations. The source lacks the P Cygni profile in the Ca II 8498Å line, part of the Ca infrared triplet (IRT), compared to V1057 Cyg and other classical FUors (Szabó et al. 2021). The strength of the blue-shifted absorption component in this P Cygni profile is formed by an outflowing wind. The weaker strength in this case could be an indication of the weaker, optically thin outflowing wind in the structure.
• During our spectroscopic monitoring despite the pure emission of the Ca II 8498Å line, we did not observe any other emission, nor forbidden emission lines in the spectra of V1515 Cyg. The absence of certain forbidden emission lines indicate that there are no jets/outflows present while these were detected in V1057 Cyg (Szabó et al. 2021).
• For our study, we obtained a new NIR spectrum of V1515 Cyg in 2021, and compared it with the one obtained in 2015 by Connelley & Reipurth (2018). In the new one, we found that both the CO overtone bandhead features and the FeH molecular bands became stronger since the previous observation. Together with the decreasing temperature from the FeH, and the increased EW of the CO bandhead in the 2021 observation, we conclude that the disk of V1515 Cyg became cooler by about 200 K since 2015, further suggesting that the source is approaching the quiescent phase.
• Despite the long-term fading of V1515 Cyg, our results from the spectroscopic monitoring and analysis still resemble the properties those of a classical/bona fide FUor.
The continuation of the photometric monitoring supplemented with additional spectroscopic snapshot observations can be crucial for our understanding of the late evolution of FUor outbursts and how they eventually end. This study is the second part of a series with the aim of revisiting the first few classical FUors after the famous 1936 eruption of FU Orionis, the prototype source. Our study highlights the importance of scrutinizing these objects which ultimately can help us to better understand the accretion processes in the early stages of Sun-like star formation.  Kausch et al. 2015), IRAF (Tody 1986(Tody , 1993 APPENDIX A. PHOTOMETRY OF V1515 CYG Table 4 contains our original optical and near-infrared photometry of V1515 Cyg before the shifts discussed in Sec. 2, while Table 5 contains the saturation corrected WISE data that we use in Fig. 1. Table 6 contains shifts applied for the different data sets in the optical bands described in Sec. 3.1