A Candidate Tidal Disruption Event in a Quasar at z = 2.359 from Abundance Ratio Variability

Figure 1 shows the two archival spectra of SDSS J1204+3518. The SDSS spectrum was taken on 05-08-2005, whereas the BOSS spectrum was taken on 03-01-2011, i.e., 1.7 years (in the quasar's rest frame) after the SDSS spectrum. The flux intensity of the N iii] λ1750 emission decayed significantly, whereas C iii] λ1909 stayed constant within uncertainties (see Figure 2 and Table 1 for details). We do not detect significant blueshift or redshift in N iii] λ1750.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

We do not detect significant N iv] λ1486 or He ii λ1640 in either epoch (see Section 5.1 for discussion). The intensity of N v λ1240 seems to have decayed (Appendix A, Figure 8), though, it is highly uncertain due to the blending with Lyα. We cannot robustly disentangle broad N v λ1240 absorption, if any, from emission due to blending and the intrinsically broad widths of the lines.

The narrow absorption lines seen in the spectra include four intervening C iv λλ1548,1551 doublet systems and two associated systems (a C iv λλ1548,1551 doublet and a N v λλ1239,1243 doublet) that have consistent redshifts (Appendix A, Figure 9). There are no clear and significant variations in the strengths or velocities of the narrow absorption lines.

The significant N/C decrease observed in SDSS J1204+3518 based on the N iii] λ1750/C iii] λ1909 emission-line intensity ratio cannot be easily explained by extreme metallicities. The timescales for the stellar population evolution necessary for metallicity changes (in quasar host galaxy stellar core or in the outer self-gravitating part of the accretion disk) are much longer than a few years.

Significant changes in the ionizing spectrum are unlikely given the similar ratios seen in the two-epoch spectra for other BLR emission lines such as the ionizing species of carbon (e.g., the C iii] λ1909/C iv λ1549 intensity ratio is an indicator of the ionization level; Shields 1976; Baldwin et al. 2003b, Table 1). Furthermore, the N iii] λ1750/C iii] λ1909 intensity ratio is sensitive to the N/C abundance ratio and is rather insensitive to the ionization parameter and the slope of the ionizing continuum (Shields 1976; Osmer 1980; Baldwin et al. 2003b; Yang et al. 2017).

Below we proceed with the hypothesis that the rapid N/C variability seen in SDSS J1204+3518 was caused by a candidate TDE of a star by the SMBH that powers the quasar.

The optical (rest-frame UV) light curves of SDSS J1204+3518 are qualitatively consistent with the TDE hypothesis, although the baseline falls short of a proof. Figure 3 shows the heterogenous set of archival data spanning photometric epochs from 2002 to 2013, encompassing the two spectroscopic epochs. There is evidence of a decaying component in the flux on top of a background of stochastic variability (at the level of ∼10%) typically seen in optical quasars. As a toy model, we fit the light-curve data using the combination of a constant background plus a decaying component as f ∝ t^α (see Sections 3.2 and 3.3 for details) where α = −5/3 is expected for TDEs for stars according to conventional TDE theory (Rees 1988). A fit to the V-band data assuming α = −5/3 and a constant background of 0.09 mJy yields a peak-luminosity date of 51,860 ± 240 (MJD), which is ∼1.3 years (rest-frame) before the "nitrogen-high" state caught in the SDSS spectrum. The fit is highly uncertain since the peak-luminosity date is unknown. A model with α = −4 (expected for TDEs of partially disrupted stars; Guillochon & Ramirez-Ruiz 2013) fits the data equally well.

Zoom In Zoom Out Reset image size

Our toy models are for the purpose of illustration only and are not meant to be proof of a TDE. Recent observations have shown that the conventional t^−5/3 profile provides a poor fit to the light curves of most low-redshift TDEs (Arcavi et al. 2014; Holoien et al. 2014, 2016b, 2016a; Brown et al. 2016, 2017; Gezari et al. 2017). A range of power-law profiles seem to fit various decline rates at different times after disruption.

While the available optical light-curve data are broadly consistent with our baseline model, the null hypothesis (i.e., purely stochastic quasar variability) cannot be easily ruled out given the limited temporal coverage and the poor photometric accuracy for most of the light-curve data. Nevertheless, the large variability amplitude observed (≳30%) is unusual for purely stochastic quasar variability given its estimated high Eddington rate (Section 3.5), since high Eddington-rate quasars are unlikely to show large variability (Rumbaugh et al. 2017). Section 5.2 discusses the brightness of the flare showing that it can be accommodated under the TDE hypothesis.

Archival X-ray observations of SDSS J1204+3518 are broadly consistent with the TDE hypothesis. The X-ray flux decayed by ∼40% over rest-frame ∼0.3 years (within 1 rest-frame year after the assumed peak luminosity date; Section 3.6), in agreement with our toy model at least qualitatively. Furthermore, An archival XMM-Newton observation taken on MJD = 52820 suggests an extremely soft X-ray spectrum ("photometric" X-ray photon index Γ_X = 4.1 estimated from the slope between the luminosities at 1 and 5 keV; Lusso & Risaliti 2016). This is similar to X-ray flares in low-redshift TDEs (Γ_X ≳ 4; e.g., Bade et al. 1996; Komossa & Greiner 1999; Lin et al. 2015; Auchettl et al. 2017, 2018), but is significantly different from typical optical quasars with Γ_X ≈ 1.9 (Young et al. 2009). We have analyzed the distribution of Γ_X using archival XMM-Newton observations for a control sample of ordinary quasars (Lusso & Risaliti 2016) that have similar redshift and luminosity to SDSS J1204+3518. SDSS J1204+3518 is an ∼4σ outlier in the distribution of Γ_X. An archival Chandra observation taken on MJD = 52476 also suggests a soft X-ray spectrum, though, the counts were too few for a robust measurement. Section 3.6 presents details on the X-ray data and analysis.

SDSS J1204+3518 is contained in the SDSS DR7 quasar catalog (Schneider et al. 2010; Shen et al. 2011). It has two spectra available from the SDSS DR13 data archive. The first epoch spectrum⁷ (Plate = 2089, Fiber ID = 328, MJD = 53498) was from the SDSS-I/II survey (York et al. 2000) and the second epoch⁸ (Plate = 4610, Fiber ID = 652, MJD = 55621) was taken by the BOSS spectrograph within the SDSS-III survey (Dawson et al. 2013). The SDSS spectrum covers the wavelength range of 3800–9200 Å with a spectral resolution of R = 1850–2200, whereas the BOSS spectrum covers 3650–10400 Å with a similar spectral resolution.

To measure the intensity ratio N iii] λ1750/C iii] λ1909 ("NCR" for short, defined in Equation (2)) as an indicator of the N/C abundance ratio, we fit spectral models to the observed spectra following the procedures as described in detail in Shen & Liu (2012). The model is a linear combination of a power-law continuum, a pseudo-continuum constructed from Fe II emission templates, and single or multiple Gaussians for the emission lines. As the uncertainties in the continuum model may induce subtle effects on measurements of the weak emission lines, we first perform a global fit to the emission-line free region to better quantify the continuum components. After subtracting the continuum we then fit multiple Gaussian models to the emission lines around the N iii] λ1750 and C iii] λ1909 region locally. The C iv λ1549 region is also fit locally but separately to measure the intensity ratio C iii] λ1909/C iv λ1549 ("CCR" for short, defined in Equation (3)) as an indicator of the ionization parameter or the hardness of the ionizing continuum. We also fit the Mg ii λ2800 region for virial black hole mass estimates from the BOSS spectrum.

More specifically, we model the N iii] λ1750 line with a single Gaussian whose center and width are free parameters. We fit the wavelength range rest-frame 1700–1970 Å. We model the C iii] λ1909 line using up to two Gaussians. We adopt two Gaussians for the Si iii] λ1892 and Al iii λ1857 lines which are often blended with the blue side of C iii] λ1909. To break degeneracy in decomposing the C iii] λ1909 complex, we tie the centers of the two Gaussians for C iii] λ1909 (i.e., the C iii] λ1909 emission-line profile is constrained to be symmetric); we also tie the velocity offsets of Si iii] λ1892 and Al iii λ1857 (relative to C iii] λ1909) to the laboratory values. We also adopt three additional Gaussians for the other possibly detectable emission lines in the region (N iv] λ1718, Fe II UV 191 at 1786.7 Å, and Si ii λ1816; Baldwin et al. 2003b). For the C iv λ1549 line, we fit the rest-frame wavelength range 1500–1700 Å. We model the C iv λ1549 line using two Gaussians. We adopt four additional Gaussians for the narrow and broad components of He ii λ1640 and O iii] λ1663. For the Mg ii λ2800 line (covered only in the BOSS spectrum but not in the SDSS spectrum), we fit the wavelength range rest-frame 2700–2790 Å. We model the Mg ii λ2800 line using a combination of up to two Gaussians for the broad component and one Gaussian for the narrow component. We impose an upper limit of 1200 km s⁻¹ for the FWHM of the narrow lines. As an example, Figure 4 shows our spectral decomposition modeling for SDSS J1204+3518. Figure 8 shows the spectral fit around the N v λ1240 region, which is highly uncertain because of blending. Figure 9 shows the narrow absorption line systems seen in the spectra of SDSS J1204+3518.

Zoom In Zoom Out Reset image size

SDSS J1204+3518 has available light-curve data from the Catalina Real-time Transient Survey (CRTS; Drake et al. 2009, V band), the Palomar Transient Factory (PTF; Law et al. 2009, R band), and the SDSS (u, g, r, i, and z bands). The SDSS measurements have the smallest photometric uncertainties, but there is only one photometric epoch and two subsequent spectroscopic epochs. We estimate the synthetic flux density in the corresponding SDSS filter calculated from convolving the SDSS and BOSS spectra with the filter throughput curves. We adopt the CRTS V-band magnitudes (converted to flux measurements) for fitting the light-curve model because they cover the longest time baseline that encompasses the two spectroscopic epochs. The PTF data do not provide more temporal coverage but serve as a double check for systematics. To mitigate the large photometric uncertainties associated with the CRTS and PTF data, we focus on the yearly inverse-variance-weighted-mean values in the analysis. We provide all the available photometry data in the literature in Table 2. We have also checked the available LINEAR (Sesar et al. 2011) data (re-calibrated to SDSS r band) that extends the temporal coverage to earlier epochs, but SDSS J1204+3518 is close to the detection limit of the survey and its photometric uncertainties are too uncertain to be included in the analysis.

The V-band corresponds to the rest-frame UV that best traces the mass accretion rate in TDEs. We fit the light curve with C + f(t) where C is a constant and f(t) is a power law of the form

$$\begin{eqnarray}&&f{(t)={f}_{0}\times (t+{t}_{0})}^{\alpha }.\end{eqnarray} \tag{ 1 }$$

A fit to the nine data points assuming α = −5/3 and a fixed constant background of C = 0.09 mJy (estimated using the latest CRTS epochs) yields a best-fit model of t₀ = −51,860 ± 240 (MJD) and f₀ = 7.5 ± 1.5 Jy (χ²/ν = 4.3 for ν = seven degrees of freedom). Assuming α = −5/3 and a positive constant background yields a model of C = 0.00 mJy, t₀ = −39,500 ± 1500 (MJD), and f₀ = (1.0 ± 0.2) × 10³ Jy (χ²/ν = 1.4 for ν = six degrees of freedom). While the latter model is statistically improved, we prefer the former model since the background quasar emission should be nonzero. The implied shorter evolutionary timescale is also more consistent with the variability timescale suggested by the two-epoch spectra. Allowing α to vary instead (assuming a fixed constant background of C = 0.09 mJy) does not provide a statistically significant improvement to the fit (yielding χ²/ν = 4.0 for ν = six degrees of freedom for a model of α = −3.97 ± 0.05).

Assuming a constant background quasar emission of 0.09 mJy, we estimate that the V-band brightness of the candidate TDE is ∼0.04 mJy at MJD ∼ 53,000. At the redshift of z = 2.359, the implied absolute V-band magnitude is M_V = −23.4 mag, assuming a Friedmann–Robertson–Walker cosmology with Ω_m = 0.3, Ω_Λ = 0.7, and h = 0.7. A K-correction of 3.1 mag has been applied assuming that the unknown SED follows a black body with T = 2 × 10⁴ K which implies a rest-frame V-band apparent magnitude of ∼2.4 μJy at MJD ∼ 53,000. Assuming that the V-band bolometric correction is in the range of 1.8–60 (appropriate for black body temperatures of 10⁴ – 5 × 10⁴ K as observed in known TDEs; Tadhunter et al. 2017), the bolometric luminosity of the candidate TDE is estimated to be in the range of 9.5 × 10⁴⁴–3.2 × 10⁴⁶ erg s⁻¹ at MJD ∼ 53,000. The total energy released during the observed part of the candidate TDE flare (using a simple trapezium rule integration of the V-band light curve) was estimated as ∼1 × 10⁵² – 4 × 10⁵³ erg. The implied total mass accreted by the black hole (from the candidate TDE only) is estimated as ∼0.02–0.5M_⊙, assuming an accretion disk efficiency of η = 0.42 for a Kerr black hole (Leloudas et al. 2016). This is a lower limit because the peak luminosity is not covered by the available light-curve data.

We estimate the black hole mass using the single-epoch estimator assuming virialized motion in the broad-line region clouds (Shen 2013). Spectral fit to the SDSS spectrum suggests a C iv λ1549 -based virial mass of M_• = 10^9.6±0.1 M_⊙ (statistical error only) using the calibrations of Vestergaard & Peterson (2006). The BOSS spectrum suggests a C iv λ1549-based virial mass of M_• = 10^9.5±0.1 M_⊙, or a Mg ii λ2800-based virial mass estimate of M_• = 10^9.3±0.1 M_⊙, using the calibrations of (Vestergaard & Osmer 2009) for Mg ii λ2800. We adopt the BOSS estimates as our baseline values because they are more representative of the nitrogen-low (i.e., more steady accretion) state of the black hole. Mg ii λ2800-based masses are generally considered more reliable than C iv λ1549-based masses (e.g., Shen 2013), given the larger scatter between C iv λ1549 and Hβ masses for high-redshift quasars (e.g., Shen & Liu 2012). C iv λ1549 is more subject to nonvirial motion such as outflows. We estimate the Eddington ratio as L_bol/L_Edd, where the Eddington luminosity is L_Edd = M_•/M_⊙ × 10³⁸ erg s⁻¹ and the bolometric luminosity L_bol is calculated from L₁₃₅₀ using the bolometric correction BC₁₃₅₀ = 3.81 from the composite quasar SED of Richards et al. (2006).

There are three X-ray observations of the field of SDSS J1204+3518 available in the public data archives. The ROSAT All Sky Survey (RASS; Voges et al. 1999, 2000) scan with the Position Sensitive Proportional Counter (Pfeffermann et al. 1987) from 1990 November yielded a 3σ upper limit of 0.02 counts s⁻¹, corresponding to an unabsorbed X-ray flux limit of F_{0.1–2.4 keV} < 2.9 × 10⁻¹³ erg cm⁻² s⁻¹ (assuming a Γ_X = 2 power-law spectrum corrected for the Galactic column density⁹ of 1.7 × 10²⁰ cm⁻²; Dickey & Lockman 1990).

A Chandra observation of HS 1202+3538 (observation ID = 3070) on 2002 July 21 (MJD = 52476) serendipitously caught SDSS J1204+3518 at an off-axis angle of 10.7 arcmin. The observed-frame soft (0.5–2 keV) and hard (2–8 keV) counts in the 6.7 ks (4.5 ks corrected for vignetting) exposure were 16 and 9, respectively (Gibson et al. 2008). We derive an unabsorbed X-ray flux (corrected for the Galactic column density) of F_{0.5–10 keV} = 2.6 × 10⁻¹⁴ erg cm⁻² s⁻¹ for Γ_X = 4, or F_{0.5–10 keV} = 4.4 × 10⁻¹⁴ erg cm⁻² s⁻¹ for Γ_X = 2. The X-ray counts are too few for a robust measurement of Γ_X. Nevertheless, the hardness ratio HR ≡ (H−S)/(H+S) = −0.28 suggests a soft X-ray spectrum, similar to the values observed in the X-ray flares of TDEs not long after peak luminosity (Auchettl et al. 2017).

Finally, SDSS J1204+3518 is contained in the 3XMM serendipitous source catalog DR5 (Rosen et al. 2016) with DETID = 101487424010032. The XMM-Newton European Photon Imaging Camera detected SDSS J1204+3518 on 2003 May 23 and 2003 June 30. The longer exposure taken on 2003 June 30 (observation ID = 0148742401) covered SDSS J1204+3518 with a "live" exposure time of 20.1 ks at an off-axis angle of 10.0 arcmin. The X-ray photon index estimated from the slope between the luminosities at 1 and 5 keV is 4.09 with an X-ray S/N = 4.7 (Lusso & Risaliti 2016). This "photometric" X-ray photon index is not as reliable as the spectroscopic tracer, but provides a reasonable estimate of the X-ray spectral hardness. The derived unabsorbed X-ray flux (averaged over the two XMM-Newton epochs) is F_{0.5–10 keV} =1.5 × 10⁻¹⁴ erg cm⁻² s⁻¹ for Γ_X = 4, or F_{0.5–10 keV} = 1.9 ×10⁻¹⁴ erg cm⁻² s⁻¹ for Γ_X = 2. For the Chandra and XMM-Newton detections, we adopt the X-ray estimates from assuming Γ_X = 4 as our baseline values (shown in Figure 3), although the qualitative decaying trend still holds for Γ_X = 2.

SDSS J1204+3518 was covered by the FIRST survey footprint but was undetected with a 3σ upper limit of ${f}_{6\,\mathrm{cm}}^{\mathrm{obs}}\lt 0.381$ mJy. Assuming that the radio flux follows a power law (i.e., ${f}_{\nu }\propto {\nu }^{\alpha }$), this translates into ${f}_{6\,\mathrm{cm}}^{\mathrm{rest}}\lt 0.70$ mJy for a spectral index α = −0.5 (Jiang et al. 2007), or <1.0 mJy assuming α = −0.8 (Gibson et al. 2008). Combined with the f₂₅₀₀ measurement from the optical spectrum, the implied limit on the radio loudness parameter (Kellermann et al. 1989), i.e., R ≡ f_{6 cm}/f₂₅₀₀, is <4.9 (<7.8) using the SDSS (BOSS) spectrum assuming α = −0.5, or <7.2 (<11) assuming α = −0.8, which is marginally inconsistent with the radio-loud criterion R > 10.

We start with a sample of 311 N-rich quasars by combining two N-rich quasar catalogs in the literature (Jiang et al. 2008; Batra & Baldwin 2014). The first catalog (Jiang et al. 2008) contains 293 quasars at 1.7 < z < 4.0 with strong N iv] λ1486 or N iii] λ1750 emission lines (rest-frame EW > 3 Å). The second catalog contains 43 quasars (of which 25 overlap with the first catalog) at 2.29 < z < 3.61 that have the strongest N iv] λ1486 and N iii] λ1750 lines in addition to strong N v λ1240 lines. Both catalogs were selected from the fourth edition of the SDSS Quasar Catalog (Schneider et al. 2007) based on the SDSS fifth data release (Adelman-McCarthy et al. 2007).

To study spectroscopic variability for these N-rich quasars, we analyze all the available high-quality spectrum pairs (with median S/N > 10 pixel⁻¹ over rest-frame 1700–2000 Å) in the SDSS thirteenth data release (SDSS Collaboration et al. 2016) for a sample of 82 unique quasars. We visually examine the ratio spectrum (similar to that shown in the lower panel of Figure 2 but over the entire relevant spectral range) for each spectrum pair; SDSS J1204+3518 was identified as our best candidate for having significant N/C variability (i.e., the ratio spectrum changed significantly over the N iii] λ1750 region but stayed constant over the C iii] λ1909 region). There are 193 high-quality spectrum pairs with rest-frame time separations >1 year (for 78 of the 82 unique N-rich quasars; the other 4 only have high S/N repeat spectra separated by <1 year). This "long-term" group serves as our parent sample for estimating the fraction of N-rich quasars that show significant N/C variability over ≳1 year timescales. To calibrate measurement uncertainties, we also define a control sample of a "short-term" group containing 355 high-quality spectrum pairs for N-rich quasars with rest-frame time separations <1 month. The median time baseline $\langle {\rm{\Delta }}t\rangle$ in the "long-term" ("short-term") group is 2.7 (0.034) years (rest frame).

Figure 5 (upper left panel) shows the distribution of the fractional variability in the N iii] λ1750/C iii] λ1909 emission-line intensity ratio ("NCR" for short). We define the fractional variability as

$$\begin{eqnarray}&&{\rm{\Delta }}\mathrm{NCR}/\mathrm{NCR}\equiv ({\rm{N}}/{{\rm{C}}}_{2}-{\rm{N}}/{{\rm{C}}}_{1})/({\rm{N}}/{{\rm{C}}}_{1}),\end{eqnarray} \tag{ 2 }$$

where N/C₁ is from the earlier epoch and N/C₂ is from the later epoch. Negative values mean a decrease in N/C over time. Monte Carlo simulations suggest that the observed distribution of the short-term NCR variability is consistent with being purely induced by statistical noise. The dispersion in ΔNCR/NCR for the long-term pairs is significantly larger than that in the short-term pairs (with the standard deviation (SD) being 1.2 for the long-term sample compared to 0.14 for the short-term sample; Table 3). This much larger scatter is most likely caused by systematic uncertainties associated with continuum modeling, which is larger in the long-term pairs due to increased quasar variability (e.g., MacLeod et al. 2010; Morganson et al. 2014).

Zoom In Zoom Out Reset image size

Table 3.  Emission-line Intensity Ratio Variability Statistics

Sample	Short-term (<1 mon)		Long-term (>1 years)
Measurements	($\langle {\rm{\Delta }}t\rangle =0.4$ mon)		($\langle {\rm{\Delta }}t\rangle =2.7$ years)
Statistics	Median	SD	Median	SD	Pnull
ΔNCR/NCR	−0.01 ± 0.01	0.14	−0.12 ± 0.08	1.2	10−21
ΔCCR/CCR	0.00 ± 0.01	0.18	0.06 ± 0.06	0.79	10−8
ΔNCR/σΔNCR	−0.2 ± 0.1	1.4	−1.0 ± 0.2	2.8	10−10
ΔCCR/σΔCCR	0.0 ± 0.2	3.1	0.5 ± 0.2	2.2	0.05

Note. NCR stands for the N iii] λ1750/C iii] λ1909 intensity ratio and CCR stands for the C iii] λ1909/C iv λ1549 intensity ratio. The difference is defined as the later epoch minus the earlier epoch so that negative values mean a decrease in the intensity ratios over time. The long-term (short-term) sample contains 193 (355) pairs of high-quality spectra for N-rich quasars. SD: standard deviation. P_null: null probability from the Kolmogorov–Smirnov test of the long-term (>1 years) sample compared against the short-term (<1 month) sample. See Section 4.2 for details.

Download table as:  ASCIITypeset image

Figure 5 also shows the effective "S/N" of the NCR and CCR variation (lower panels), i.e., emission-line intensity ratio variation normalized by measurement uncertainties (accounting for both statistical and systematic errors). The dispersion in ΔNCR/σ_ΔNCR for the long-term pairs is more similar to that in the short-term pairs after the normalization (with SD being 2.8 for the long-term sample compared to 1.4 for the short-term sample), that accounts for differences in the systematic uncertainties. The long-term sample shows a 5σ detection that the median value is nonzero (−1.0 ± 0.2, where 0.2 is the 1σ error in the median, estimated as SD$/\,\sqrt{N}$, where N is the number of data points; Table 3). Furthermore, its probability distribution is significantly different from that expected from pure measurement noise as characterized by the short-term pairs. Kolmogorov–Smirnov test shows that the probability that the two distributions are drawn from the same sample P_null is 10⁻¹⁰ (Table 3). Compared to the control sample of short-term pairs, the long-term pairs show a net negative ΔNCR/σ_ΔNCR, providing statistical evidence that NCR decreases with time in N-rich quasars. SDSS J1204+3518 stands out among those that show the most significant decreases in the NCR. We have checked the other positive and negative value tails seen in the ΔNCR/NCR distribution, which are largely caused by systematic uncertainties from our emission-line measurements and do not show significant variability in the flux-ratio spectrum.

Figure 5 (upper right panel) also shows the C iii] λ1909/C iv λ1549 emission-line intensity ratio ("CCR" for short) variation. Similarly to the NCR variation, it is defined as

$$\begin{eqnarray}&&{\rm{\Delta }}\mathrm{CCR}/\mathrm{CCR}\equiv ({\rm{C}}3/{\rm{C}}{4}_{2}-{\rm{C}}3/{\rm{C}}{4}_{1})/({\rm{C}}3/{\rm{C}}{4}_{1}),\end{eqnarray} \tag{ 3 }$$

where C3/C4₁ is from the earlier epoch and C3/C4₂ is from the later epoch. The long-term CCR variation, on the other hand, is more similar to the short-term CCR variation. While the dispersion is still larger for the long-term CCR variation in terms of ΔCCR/CCR (SD being 0.79 for the long-term compared to 0.18 for the short-term sample; Table 3), their difference is smaller than that in the dispersion of ΔNCR/NCR (a factor of ∼4 difference rather than ∼9). This is understandable considering that stronger lines (i.e., C iv λ1549) are less affected by systematic uncertainties in the continuum modeling. In terms of ΔCCR/σ_ΔCCR, i.e., after normalizing the measurement uncertainties, the median value of the long-term sample is consistent with being zero within 2.5σ (0.5 ± 0.2; Table 3); the 2.5σ positive deviation from zero is most likely explained by a noise-induced bias, as we discuss below. Furthermore, the long-term and short-term distributions are statistically identical (with P_null being 0.05; Table 3). SDSS J1204+3518 does not show significant CCR variation. Table 3 summarizes the statistical properties of the NCR and CCR variability measurements.

Figure 6 examines whether there is any general trend between ΔNCR/NCR and ΔCCR/CCR. Although ΔCCR/CCR appears to be increasing as ΔNCR/NCR decreases in the long-term sample (shown with cyan dots), this apparent trend is also seen in the short-term sample (shown with black crosses), which is most likely caused by the correlation between the NCR and CCR terms due to measurement errors. The apparent trend arises because NCR and CCR are not independent from each other since the C iii] λ1909 emission-line flux intensity is included in both terms. Considering measurement errors, a stronger C iii] λ1909 (caused by noise) will lead to a smaller NCR and a larger CCR, resulting in an apparent anti-correlation between the NCR and CCR. This also likely explains the 2.5σ positive deviation from zero in the median value of ΔCCR/σ_ΔCCR for the long-term sample, in which systematic uncertainty from continuum modeling adds to the error, resulting in more bias in the long-term sample than in the short-term sample because of increased quasar variability.

Zoom In Zoom Out Reset image size

Figure 7 examines the potential time dependence of the intensity ratio variability. We show both the individual data points (small dots) and the median values in a given bin of time separation (large open squares with error bars denoting uncertainties in the median values) to help assess the general trend. While there is some hint of the variability increasing as a function of time in support of the TDE hypothesis, we cannot draw a firm conclusion in view of the limited dynamic range in time separation and the large scatter in the intensity ratio variability measurements.

Zoom In Zoom Out Reset image size

The absence of evidence for strong variability in the N iv] λ1486 and He ii λ1640 lines is puzzling, but does not rule out the TDE hypothesis. In the only three known TDEs with available UV spectra (Cenko et al. 2016; Yang et al. 2017; Brown et al. 2018), N iv] λ1486 (He ii λ1640) is not always detected or relatively weak in some epochs, whereas N iii] λ1750 is always detected. So empirically it is possible to see N iii] λ1750 but not N iv] λ1486 (He ii λ1640) emission in TDEs, which may be related to the different ionizing potentials of the lines. In addition, a practical complication is that N iv] λ1486 is on the extended blue wing of C iv λ1549 for SDSS J1204+3518 and may be buried underneath, especially if it is weak. For He ii λ1640, another possibility is that if the disrupted star were about to ascend the giant branch, most helium would be sequestered in the dense helium core, which would not be disrupted by the more massive black hole (unlike in low-redshift TDEs around smaller black holes). Finally, our selection is based on the N iii] λ1750/C iii] λ1909 ratio variability so that we may have missed systems that show stronger N iv] λ1486 and/or He ii λ1640 variability.

While bearing some resemblances to N-rich quasar spectra, none of the three TDEs with UV spectra exhibit C iii] λ1909 or Mg ii λ2800 emission features. Their C iv λ1549 is also much weaker than expected. The origins of these discrepancies are still open to debate. They might be due to differences in the physical conditions in the gas and on the shape of the ionizing continuum. For example, the absence of Mg ii λ2800 could result from photoionization by the extremely hot continuum, which may be transient in nature as the continuum temperature eventually cools down (Cenko et al. 2016). The absence of C iii] λ1909 could be due to collisional deexcitation if the gas density is above critical, n = 10^9.5 cm⁻³ (Osterbrock 1989).

The brightness of the candidate TDE of ∼0.04 mJy at MJD ∼ 53,000 implies an absolute V-band magnitude of −23.4 mag at the redshift of the quasar. This is only ∼0.6 mag fainter than the peak luminosity of ASASSN-15lh (Leloudas et al. 2016), the brightest TDE candidate reported to date (see Dong et al. 2016; Godoy-Rivera et al. 2017, however, for an alternative explanation for the nature of the source as a supernova). This high apparent brightness seems difficult to explain under the TDE hypothesis, but it is perhaps understandable for the following reasons. First, the observed V-band samples the rest-frame UV that is close to the peak of the spectral energy distribution (SED) of the candidate TDE, resulting in a large K-correction at the redshift of the quasar. This is dramatically different from local TDEs, in which the V band only samples the Rayleigh–Jeans tail of the SED.

Second, the high luminosity of the candidate TDE could be explained by the large BH mass of the quasar and its likely rapid spin (see Sections 3.5 and 5.3 for details). In addition, TDEs in AGN hosts may have higher radiative efficiencies than those in inactive galactic nuclei potentially due to interaction of the stellar debris with the preexisting accretion disk (Blanchard et al. 2017). However, the estimated efficiency does not strain the Eddington limit even at the projected peak-luminosity date, which is highly uncertain. The total energy released would require the accretion of a mass of ≳0.02–0.5M_⊙ (see Section 3.4 for details).

Finally, our sample is most likely highly biased to the most luminous events because: (i) by selection we are sensitive only to the high-redshift (i.e., z > 1.7 for the relevant UV lines to be covered in the observed optical spectra) population, and (ii) the SDSS is not a particularly deep survey.

Analysis of the BOSS spectrum suggests that the quasar is powered by an SMBH with a virial mass estimate of

$$\begin{eqnarray}&&\mathrm{log}\,({M}_{\bullet }/{M}_{\odot })=9.3\pm 0.1\pm 0.5\,\mathrm{dex},\end{eqnarray} \tag{ 4 }$$

where 0.1 (0.5) dex is the 1σ statistical (systematic) error (Shen 2013). The implied Eddington ratio is ∼0.005–0.16 (for the candidate TDE flare component only, depending on the unknown SED; it is ∼0.13 for the background quasar emission; see Section 3.5 for details) around MJD of 53,000 (i.e., ∼1 rest-frame years post the peak-luminosity date), consistent with an initial sub-Eddington phase predicted by conventional TDE theory for massive BHs (Rees 1988). Assuming a time lag of

$$\begin{eqnarray}&&{t}_{f}\sim 1.1{({M}_{\bullet }/{10}^{8}{M}_{\odot })}^{1/2}{r}_{\star }^{3/2}{m}_{\star }^{-1}\,\mathrm{years},\end{eqnarray} \tag{ 5 }$$

where r_⋆ ≡ R_⋆/R_⊙ and m_⋆ ≡ M_⋆/M_⊙, between disruption and infall of the tightest-bound debris for a solar-type star, a black hole mass of M_• = 10^9.3 M_⊙ would yield a Newtonian estimate of the disruption date of ∼46,000 (MJD), i.e., ∼5 rest-frame years before the peak-luminosity date. A relativistic estimate of the disruption date would be ∼50,000 (MJD), i.e., ∼2 rest-frame years before the peak-luminosity date (Leloudas et al. 2016).

The gravitational radius increases with black hole mass at a higher rate than the tidal radius R_T does (i.e., R_g ∝ M_• whereas ${R}_{{\rm{T}}}\propto {M}_{\bullet }^{1/3}$). Therefore, stars can be tidally disrupted before being swallowed whole into the horizon only if the black hole is less massive than the Hills mass M_Hills, where

$$\begin{eqnarray}&&{M}_{\mathrm{Hills}}\approx {10}^{8}{M}_{\odot }{r}_{\star }^{3/2}{m}_{\star }^{-1/2}\end{eqnarray} \tag{ 6 }$$

for non-spinning Schwarzschild black holes (Hills 1975). M_Hills can be up to an order of magnitude larger for stars on optimal orbits around rapidly spinning Kerr black holes (Leloudas et al. 2016). The virial mass estimate of SDSS J1204+3518 suggests that a spinning Kerr black hole is required for all allowed masses. A main-sequence supersolar mass star on prograde equatorial orbits can be disrupted by a maximally spinning Kerr black hole with M_• = 10^8.9 M_⊙, or a moderately spinning Kerr black hole with M_• = 10^8.5 M_⊙. If the BH mass is in the higher end in the estimated range, then a giant star is required.

The rapid N/C variability, on the other hand, implies that the evolutionary timescale is not much longer than a few years. Assuming that it is driven by the circularization luminosity (but see, e.g., Jiang et al. 2016, for an alternative explanation) that evolves on the fallback timescale t_f, which is comparable to the viscous timescale in the accretion disk for massive black holes, the rapid evolution suggests that the disrupted is a main-sequence star rather than a giant. Future near-infrared spectroscopy of the Balmer lines and/or reverberation mapping observations are needed to help better constrain the black hole mass for SDSS J1204+3518.

A significant difference between N-rich and normal optical quasars is that N-rich quasars have systematically narrower widths of C iv λ1549 (Jiang et al. 2008; Batra & Baldwin 2014). This implies, under standard assumptions, that N-rich quasars have less massive black holes than the normal quasar population. This is at odds with the more massive host galaxies implied by the high metallicity hypothesis given the galaxy mass–metallicity relation (Tremonti et al. 2004; Erb et al. 2006; Maiolino et al. 2008) and the SMBH-host mass correlation (Ferrarese & Merritt 2000; Gebhardt et al. 2000; Tremaine et al. 2002; Kormendy & Ho 2013; Graham 2016). In contrast, the TDE scenario would provide a natural explanation for the narrower widths of C iv λ1549 seen in N-rich quasars since less massive black holes are more likely to disrupt stars.

N-rich quasars also have a significantly higher radio-loud fraction compared to normal optical quasars (Jiang et al. 2008). The origin of this difference is unknown. If the radio loudness is physically related to the black hole spin (Blandford & Znajek 1977; Sikora et al. 2007), the TDE hypothesis may also explain the high radio-loud fraction in N-rich quasars, since black holes with larger spins are more prone to disrupting stars. SDSS J1204+3518 was undetected by the FIRST 1.4 GHz survey (White et al. 1997) whose upper limit was marginally inconsistent with being radio loud (see Section 3.7 for details).

SDSS J1204+3518 represents the first case of a quasar with significant N/C variability over yearly timescales. Our analysis on a large sample of N-rich quasars shows statistical evidence that N/C decays over time, while the ionization level remains unchanged (Section 4). However, similar to the case of Q0353−383 (Osmer 1980; Baldwin et al. 2003b), the typical decay rate estimated in the average SDSS N-rich quasars (median ΔNCR/NCR = −12 ± 8% in ∼2.7 rest-frame years, or ∼−4 ± 3% yr⁻¹; Table 3) is much less than that seen in SDSS J1204+3518 (ΔNCR/NCR = −86 ± 14% in 1.7 rest-frame years, or approximately −50 ± 8% yr⁻¹; Table 1). This could be explained in the TDE scenario if most N-rich quasars were caused by TDEs of giants rather than main-sequence stars, since TDEs of giants would evolve more slowly (MacLeod et al. 2012). We cannot rule out the possibility that only some but not all of the N-rich quasars are caused by TDEs. Testing whether most N-rich quasars are due to TDEs of giants would require high S/N spectroscopic follow-up observations on much longer timescales (i.e., decades and even centuries) rather than being probed by current surveys.

While the lifetimes of giant stars are much shorter than those of main-sequence stars (and therefore TDEs of giants should be much less frequent), it is still possible for the N-rich population to be dominated by TDEs of giants considering the disruption constraint given the large black hole masses of high-redshift quasars (i.e., main-sequence stars are more likely to be swallowed whole rather than being disrupted unless the black hole mass is small enough and/or the spin is high), in contrast to the demographics seen in local, inactive galaxies with smaller black holes.

We estimate the implied TDE rate in quasars assuming that (i) SDSS J1204+3518 is a TDE of a main-sequence star, and (ii) nearly all of the rest of N-rich quasars are TDEs of giants. The total order-of-magnitude rate is estimated as

$$\begin{eqnarray}r & = & {r}_{\mathrm{MS}}+{r}_{{\rm{G}}}={f}_{{\rm{N}}-\mathrm{rich}}\times \left(\displaystyle \frac{{f}_{\mathrm{MS}}}{2\,\mathrm{years}}+\displaystyle \frac{{f}_{{\rm{G}}}}{{10}^{3}\,\mathrm{years}}\right)\\ & \sim & {10}^{-4}\,{\mathrm{yr}}^{-1}\,{\mathrm{galaxy}}^{-1},\end{eqnarray} \tag{ 7 }$$

where f_N-rich ≲ 10⁻², ${f}_{\mathrm{MS}}\sim \tfrac{1}{78}$, and $0\lt {f}_{{\rm{G}}}\lesssim \tfrac{77}{78}$ (${f}_{{\rm{G}}}\sim \tfrac{77}{78}$ assuming all N-rich quasars are TDEs). The evolution timescale (which goes into the denominator of the rate estimate) is much longer for TDEs of giants (MacLeod et al. 2012, here assumed to be ∼10³ years; assuming it is ∼10² years instead would make r_G comparable to r_MS but still not an order-of-magnitude larger) so that their contribution to the estimated rate is small (or similar if f_G is close to unity) compared to that by SDSS J1204+3518 alone.

Under the TDE scenario, our systematic search implies a much higher TDE rate in quasars than in normal galaxies. The fraction of N-rich quasars is ≲1%, and we found one definitive case SDSS J1204+3518 out of a parent sample of 78 N-rich quasars in total within a time window of ∼2 years (rest-frame). Thus the implied TDE rate is ∼10⁻⁴ yr⁻¹ galaxy⁻¹. This is ∼100 times higher than the expected TDE rate at M_• ∼ 10⁹ M_⊙, which is ∼10⁻⁶ yr⁻¹ galaxy⁻¹ for supersolar mass stars (Kochanek 2016b). Similarly, the observed TDE rate in post-starburst galaxies is ∼100 times higher than normal galaxies (Arcavi et al. 2014). This is perhaps not a coincidence since both quasars and post-starburst activities may be associated with recent merger events that would greatly increase the number of stars in low angular momentum orbits, fueling TDEs.

Our findings provide motivation for future research programs on N-rich quasars, which have been enigmatic and hard-to-understand objects until now. Despite their rarity, we show that N-rich quasars may be important links for understanding how SMBHs disrupt stars.

Examples of future research programs include:

• 1.  
  Finding additional cases of TDEs of stars in quasars to establish their frequency of incidence and learn more about their astrophysics.
• 2.  
  Learning more about the physics of the encounters of stars with SMBHs.
• 3.  
  Provide improved knowledge of the evolutionary nature and distributions of stars near SMBHs. For example, timescales for disruption of giant stars are longer than for main-sequence stars.

The abundance ratio variability offers a potentially new method for identifying TDEs, complementary to the traditional method based on X-ray and/or UV/optical flux variability. In particular, it may open a new window of discovering TDEs at significantly higher redshifts (z > 2) than previous work. In comparison, the majority of known TDEs are at low redshifts (z < 0.3), and the current redshift record holder is Swift J2058+05 at z = 1.1853¹⁰ (Cenko et al. 2012). This has implications for the current efforts to use TDEs to study supermassive black holes. Its potential may be better realized with the upcoming large-scale time-domain spectroscopic surveys such as the SDSS-V (Kollmeier et al. 2017) and MSE (McConnachie et al. 2016) projects.