A Short-lived Rejuvenation during the Decades-long Changing-look Transition in the Nucleus of Mrk 1018

Long-term spectroscopic observations of Mrk 1018 were conducted utilizing the Yunnan Faint Object Spectrograph and Camera (YFOSC) attached to the Lijiang 2.4 m telescope (LJT). YFOSC is quipped with a back-illuminated 2048$\times$2048 pixel CCD, featuring a pixel size is 13.5 $\mu$m, a pixel scale is 0.283${}^{\prime\prime}$ per pixel, a field-of-view is $10^{\prime}\times 10^{\prime}$), and a series of filter and Grism (Fan2015; Wang2019). For spectroscopy of Mrk 1018, Grism 14 was employed, which spans wavelengths from approximately $\sim$3600 Å to 7460 Å and offers a dispersion of 1.8 Å  pixel${}^{-1}$. Following our earlier studies (Lu2021a; Lu2021b; Lu2022), a long slit with a projected width of $2.5^{\prime\prime}$ was used, taking into account the average seeing conditions at the observatory.

The spectroscopic monitoring of Mrk 1018 started on November 5, 2017, and to date, seven-season observations have been conducted from November 2017 to February 2024 (with Modified Julian Days ranging from 58063 to 60324). Similar to the long-term spectroscopic monitoring of NGC 5548 (Lu2022), we rotated the long slit within the field of view to simultaneously capture spectra of Mrk 1018 and a nearby stable comparison star. This observational approach is commonly employed in reverberation mapping campaigns, where the spectra of the comparison star can provide high-precision spectral calibration, including flux calibration and telluric absorption corrections (e.g., see Hu2015; Lu2019a; Lu2021b and below). Standard neon and helium lamps were used for wavelength calibration. In total, we acquired 96 spectroscopic observations over an approximate duration of seven years.

We processed the two-dimensional spectroscopic data using the standard procedures in PyRAF. This involved bias subtraction, flat-field correction, wavelength calibration, and the removal of cosmic ray. Subsequently, all spectra were extracted with a consistent extraction window after subtracting the sky background. We opted for a relatively narrow extraction window of 20 pixels (5.7${}^{\prime\prime}$), which aids in minimizing Poisson noise from sky background and increasing the signal-to-noise (S/N) ratio of the spectra. The sky background was estimated from two adjacent regions ($+7.4^{\prime\prime}\sim+14^{\prime\prime}$ and $-7.4^{\prime\prime}\sim-14^{\prime\prime}$) located on either side of the extraction window.

The red side of ground-based optical spectrum is contaminated by the variable absorptions of the telluric atmosphere. In the observed spectrum of Mrk 1018, there are three telluric absorption windows, which are marked by vertical colored bands in Figure 1. One of these windows (the green band) overlaps with the broad H$\alpha$ line. We first correct the telluric absorption and subsequently calibrate the spectral fluxes.

Following the work of Lu2021b, we processed the observed spectra of Mrk 1018 through a three-step calibration procedure. First, we determined the synthesis stellar template for the comparison star by aligning its observed spectrum with a library of synthesis stellar models (e.g., Husser2013). Next, we modeled the observed spectra of comparison star (in counts) using the chosen stellar template, while masking the telluric absorption bands during this modeling (i.e., the regions marked by color diagram in Figure 1). We then derived the the telluric transmission spectra by dividing the observed spectrum of the comparison star by the modeled spectrum. This allowed us to use the telluric transmission spectra to adjust for the telluric absorption lines in each of the observed spectra. Figure 1 (left y-axis) displays the corrected spectra of telluric absorption (in red) alongside the observed spectra (in black), with the bottom panel representing the spectrum of Mrk 1018 and the top panel showing the spectrum of the comparison star. Finally, similar to the previous step, we generated a wavelength-dependent sensitivity function for each object/comparison star pair by comparing the telluric-corrected spectrum of the comparison star with the stellar template. This sensitivity function was then utilized to calibrate the spectrum of Mrk 1018 (see also Lu2019a; Lu2021a; Lu2022). The flux-calibrated spectra of Mrk 1018 (bottom panel) and its comparison star (top panel) are shown in blue in Figure 1 (right y-axis).

Before our monitoring, there exist nine spectra concurrently covering the broad H$\beta$ and H$\alpha$ lines (Kim2018; Osterbrock1981; Cohen1986; Abazajian2009; Jones2009; Trippe2010; McElroy2016), which enable us to investigate the extensive long-term spectral evolution of Mrk 1018 (see Section LABEL:blrph). We obtained these early spectra from Kim2018 and recalibrated them by assuming the same [O iii]$\lambda 5007$ flux as our spectra. The first spectra of Mrk 1018 was taken in the 1970s, meaning that our spectroscopic baseline spans more than 45 years.

The spectral fitting scheme is commonly adopted in spectral analysis of AGNs to separate the blending components (e.g., Dong2011; Park2012; Guo2014; Barth2015). In line with previous fitting routine (e.g., Lu2022), we conducted the spectral fitting and multicomponent decomposition to extract the spectral characteristics of Mrk 1018 (see Figure 2), including variability and velocity shifts of broad emission lines. To cover the broad H$\beta$ and H$\alpha$ regions while minimizing degrees of freedom, we set the fitting window from 4440 Å to 6800 Å in rest frame for all Mrk 1018 spectra. The fitting includes the following components: (1) A power law ($f_{\lambda}\propto\lambda^{\alpha}$, $\alpha$ is the spectral index) representing the AGN continuum. (2) An iron template from Boroson1992 for the iron multiplets. (3) A stellar template with an age of 11 Gyr and metallicity of Z = 0.05 from Bruzual2003 to account for the host-galaxy starlight. (4) Double Gaussians for each broad Balmer lines (H$\beta$ and H$\alpha$). In practice, we initially attempted to fit the broad Balmer line using double Lorentzians (Lorentzian was employed for the broad Balmer line of Narrow-Line Seyfert 1 galaxies, see Veron2001), but found the fits are inferior compared to those using double Gaussians for the double-peak broad lines. Therefore, we chose to use double Gaussians to model the broad Balmer lines. (5) Two sets of double Gaussians for the [O iii] doublets $\lambda 5007/\lambda$4959. (6) A single Gaussian for each narrow emission line. (7) Four sets of single Gaussian for the [S ii] doublets and the [N ii] doublets. All components were fitted simultaneously across the fitting window using the MPFIT package (Markwardt2009), which employs the Levenberg$-$Marquardt algorithm for $\chi^{2}$-minimization. During the spectral fitting for each spectrum, the flux ratio of the [O iii] doublets was fixed to the theoretical value of 3 (e.g., McGill2008; Hu2015); all the narrow emission lines were tied to the same velocity and shift, while the rest of the fitting parameters were allowed to vary. Various host-galaxy templates from Bruzual2003 were also evaluated and the template with an age of 11 Gyr and a metallicity Z = 0.05 provided a reasonable fit for the spectral index of AGNs continuum ($\sim\lambda^{-1.5}$) and the stellar absorption lines. While several iron templates have been proposed for fitting the optical or ultraviolet spectra of AGNs (e.g., Boroson1992; Veron2004; Kova2010; Park2022), we found that the spectral fitting of Mrk 1018 is not sensitive to the choice of the specific iron template.

Next, we focus on the measurements of spectral characteristics. The broad H$\beta$ and H$\alpha$ lines fluxes for each spectrum are measured from the optimal fitted models and are tabulated in Table 1 along with the uncertainties (including Poisson errors and systematic errors). The fluxes of the AGNs continuum can be measured from the best-fitted power-law component at 5100 Å. During our observation period, the pseudo-continuum is dominated by the starlight from the host galaxy.The mild degeneracy between power-law continuum and host-galaxy component will cause scatter to the obtained AGN continuum at 5100 Å. Therefore, we just use the averaged AGN continuum flux at 5100 Å as a a reference for the intercalibration of different photometric light curves (see Section 2.3). We measure the flux of [O iii]$\lambda 5007$ line from the best-fitted model and report the result in Figure 3, which exhibits a scatter of 6.9%. This level of scatter indicates a sufficient accuracy in our spectral calibration.

In spectral fitting described above, the broad H$\beta$ and H$\alpha$ lines in each observation are well modeled by two Gaussian components, This allows us to extract the characteristics of the broad Balmer lines from the optimal models. After correcting the wavelength shift of broad Balmer lines, which usually arise from varying seeing and mis-centering, using its narrow lines as wavelength reference, we measure the peak separations of two Gaussian components (including broad H$\beta$ and H$\alpha$ lines) based on the differences in peak wavelengths. Additionally, we measure the line widths including full-width at half-maximum (FWHM) and velocity dispersion ($\sigma_{\rm line}$) for broad H$\beta$ and H$\alpha$ lines. To reasonably determine the FWHM of the double-peaked emission line, we identify a blue-side peak and a red-side peak within the line profile, and the FWHM is calculate from the maximum wavelength separation of these peaks at each half-maximum, as described by Peterson2004. In practice, instrument broadening is coupled with varying atmospheric (seeing) broadening, so we estimate the total broadening of the broad emission line for each observation by comparing the width of the [O iii] $\lambda$5007 line with that from the SDSS spectrum, which is then used to derive the broadening-corrected line width. These measurements will be used to investigate the broad-line properties in subsequent analyses.

To obtain the long-term variability of Mrk 1018, we compiled the photometry data from the public archival database of Zwicky Transient Facility (ZTF; Graham2019) and Asteroid Terrestrial-impact Last Alert System (ATLAS; Tonry2018). We use the ZTF $g$- and $r$-band data, along with the ATLAS $c$- and $o$-band data, and combine them using the Python package PyCALI¹¹1https://github.com/LiyrAstroph/PyCALI (LYR2014), which employs a multiplicative factor and an additive factor to intercalibrate different sources of light curves. We adopt ZTF-$g$ as the reference dataset and align other datasets with it.

Mrk 1018 is also covered by the Wide-field Infrared Survey Explorer (WISE; Wright2010) survey. We retrieved the photometric data with the following flags: the best frame image quality score ($qi_{-}fact=1$), no contamination of scattered light from the moon ($moon_{-}masked=0$), the larger South Atlantic Anomaly separation ($saa_{-}sep>=5$), and spurious detections excluded ($cc_{-}flags=0$). When necessary, we convert $W1$ and $W2$ magnitudes into flux densities by adopting zero-point flux densities $F_{\lambda}(W1=0)=8.03\times 10^{-12}~{}{\rm erg~{}s^{-1}cm^{-2}\AA^{-1}}$ and $F_{\lambda}(W2=0)=2.43\times 10^{-12}~{}{\rm erg~{}s^{-1}cm^{-2}\AA^{-1}}$ (Jarrett2013). The resulting fluxes and magnitudes are rebinned every 6 months.