Probing the origin of the two-component structure of broad-line region by reverberation mapping of an extremely variable quasar

ABSTRACT

The physical origins of quasar components, such as the broad-line region (BLR) and dust torus, remain under debate. To gain insights into them, we focused on changing-state quasars (CSQs) which provide a unique perspective through structural changes associated with accretion disc state transitions. We targeted SDSS J125809.31+351943.0, an extremely variable CSQ, to study its central core structure and kinematics. We conducted reverberation mapping with optical spectroscopy to explore the structure of the BLR and estimate the black hole mass. The results from Hβ reverberation mapping indicated a black hole mass of |$10^{9.64^{+0.11}_{-0.20}}\rm {M_\odot }$|⁠. Additionally, we analysed variations in the optical to X-ray spectral indices, α_ox, before and after the state transition, to investigate the accretion disc. These variations in α_ox and the Eddington ratio (from 0.4  per cent to 2.4  per cent) exhibitied behaviour similar to state transitions observed in X-ray binary systems. Spectral analysis of Hβ revealed a predominantly double-peaked profile during dim periods, transitioning to include a single-peaked component as the quasar brightened, suggesting that Hβ contains a mixture of two components. Each of these components has its distinct characteristics: the first is a double-peaked profile that remains stable despite changes in the accretion rate, while the second is a variable single-peaked profile. Using time-lags from reverberation mapping, we estimated the spatial relationships between these BLR components, the accretion disc, and the dust torus. Our results suggest that the BLR consists of two distinct components, each differing in location and origin.

accretion, accretion discs, quasars: emission lines, quasars: individual: SDSS J125809.31+351943.0, quasars: supermassive black holes

1 INTRODUCTION

The broad-line region (BLR) is one of the key components to solve several significant problems in astronomy. First, knowing the mechanism of the black hole’s gaining mass will lead to understanding the formation process of supermassive black holes (SMBHs) and the co-evolution with the host galaxy (Magorrian et al. 1998). The information on the structure and kinematics of the BLR is essential to investigate the black hole evolution because it is generally assumed when astronomers estimate the black hole mass of distant quasars. Secondly, the precise picture of quasars is also valuable for measuring the distances from our galaxy. A controversial issue in observational cosmology is the Hubble tension problem (see Perivolaropoulos & Skara 2022, as a recent review), where the measured Hubble constant (H₀) with cosmic microwave background (Planck Collaboration VI 2020) and with the distance ladder (e.g. Riess et al. 2022) differ by more than 5σ. In this context, distance measurements using the BLR of quasars (e.g. Wang et al. 2020a) are expected to be a clue to solving the problem, as they can estimate H₀ independently of traditional methods. For these reasons, we believe it is important to know about the structure of the BLR.

Broad emission lines (BEL) are major characteristics of Type 1 quasars. A picture based on the unified model (Antonucci 1993; Urry & Padovani 1995), in which dust is distributed in a torus structure around the accretion disc and the BLR, is widely accepted as the general structure of quasars. This simple picture can explain diverse spectra of quasars, Type I and II, only by one parameter of viewing angle. The detailed structures, however, are still under discussion as well as the motion and origin of the BLR (see Gaskell 2009, as a review).

Reverberation mapping (RM) is one of the most promising methods for studying the quasar structure (see Cackett, Bentz & Kara 2021, and references therein). Quasars randomly change luminosity at all wavelengths in diverse timescales and amplitudes. Using that feature, RM is a method to estimate the size of the BLR by measuring the response time between the variation in the continuum emitted from the accretion disc and the variation in the BEL emitted from the BLR. The virial black hole mass can be calculated by the measured size of the BLR and its velocity, which is thought to be the most precise method for measuring the black hole mass of distant quasars. An important result of the RM is the discovery of a strong correlation between luminosity and the BLR radius (L-R relation; Koratkar & Gaskell 1991; Kaspi et al. 2000, 2005), which has been widely used to estimate black hole mass from a single spectrum (e.g. Shen et al. 2011). RM can provide more detailed information by increasing the velocity resolution and observation frequency. As usual, emission lines are often divided into several velocity components to estimate the qualitative structure and kinematic information of the BLR (velocity-resolved RM; Du et al. 2016; De Rosa et al. 2018). By using the information derived from each component, it will be possible to elucidate the structure and kinematics of the BLR. For example, if the red wing of the BEL lags behind the blue one, we can infer that outflow motion is dominant in the BLR. To investigate the properties of the BLR in detail, we adopted RM in this study.

Understanding the properties of the BLR is necessary for accurate estimation of the black hole mass (Peterson & Wandel 1999, 2000), which is a fundamental physical quantity of a quasar, because measurement of the black hole mass by RM depends on the structure, mainly the viewing angle, of the BLR. In RM, the black hole mass is calculated with the following equation;

$$\begin{eqnarray} M_{\rm {BH}} = f \frac{c \tau \Delta V^2}{G}, \end{eqnarray}$$

(1)

where c is the speed of light, G is the gravitational constant, and ΔV is the velocity of the BLR gas. The f is a factor that reflects the structure, kinematics, and inclination of the BLR. The problem is that the f factor is supposed to take different values for each object, because at least inclination is different from object to object, while it is currently assumed as the same value. Furthermore, in many cases, the radius of the BLR is estimated from a single spectroscopic observation using the empirical L-R relation. We express concern that discussions premised on such a series of assumptions may potentially lead to misleading interpretations. To solve these problems, it is necessary to create models that are closer to reality and check them against observations, but there are many observational uncertainties. Therefore, the purpose of this study is to gain insight into the structure of the BLR.

To understand the structure and kinematics of the BLR in detail, consideration of the origin of the ionized gas in the BLR is important. There are two major hypotheses for the origin of the BLR gas. One is from the accretion discs as the disc wind (e.g. Czerny & Hryniewicz 2011) and the other is the inflow (including tidally disrupted clumps) from the dust torus (e.g. Gaskell 2009). On one hand, if the gas is supplied by the disc wind from the accretion disc, it is expected that the amount of disc wind changes with the mass accretion rate of the accretion disc, resulting in a variation in the BEL intensity (Elitzur, Ho & Trump 2014). On the other hand, some studies have suggested that tidally disrupted dust torus gas is responsible for the asymmetric component of the BEL (Wang et al. 2017). As the observational results are different from object to object, the origin of the BLR is still under discussion. For example, the velocity-resolved RM results revealed that some objects exhibited outflow features (e.g. Lu et al. 2019, 2021) while some exhibited inflow features (e.g. Gaskell 1988; Koratkar & Gaskell 1989; Denney et al. 2010; Grier et al. 2013a; Du et al. 2016; Zhang et al. 2019; Lu et al. 2021; Bao et al. 2022), or both (e.g. Du et al. 2018; Brotherton et al. 2020; Lu et al. 2021). A possible scenario is that multiple origins may be mixed. The key is to sort out what could exist as an origin and to what extent it is dominant.

We suggested in Nagoshi & Iwamuro (2022) that an ideal target for investigating the origin of the BLR is changing-state (look) quasars¹ (CSQs). We focus on the variability of BEL and continuum luminosity, which show significant variations (⁠|$\gt 30~{{\ \rm per\ cent}}$|⁠) in the BELs and in mid-infrared (mid-IR) luminosity (LaMassa et al. 2015; MacLeod et al. 2016, 2019; Ruan et al. 2016; Gezari et al. 2017; Noda & Done 2018; Wang, Xu & Wei 2018; Graham et al. 2020; Ross et al. 2020; Nagoshi et al. 2021; Wada et al. 2021; Green et al. 2022). This significant variation represents the structural transformation process, which is an ideal feature for investigating each component of the unified model. In particular, if we can observe the emergence process of the BLR, we can place restrictions on the origin of the BLR gas. Although the variation process of accretion disc also provides essential information, the observational studies are still not enough for quasars because of their long variation time-scale. Therefore, the purpose of this study is to apply RM to one of the CSQs (Dexter et al. 2019; Li et al. 2022) together with X-ray observations to gain insight into the physical origin of the BLR and the physical process of drastic variation in the accretion disc.

In this paper, we study the BLR structure of SDSS J125809.31+351943.0 (hereafter J1258; z = 0.30939), which was discovered by Nagoshi et al. (2021) as the object that experienced the most significant state transition in history (optical brightness variation of about 4 mag in about 30 yr). The significant brightness variation is expected to have a large signal-to-noise ratio of the varied component in the spectrum, making it a suitable target for RM.

Each section is organized as follows. Section 2 describes our observations and the archival data we used. Section 3 describes the results of the RM analysis and the comparison of SED before and after the extreme variation. Section 4 discusses the black hole mass, the cause of the variation, and the quasar structure based on the above results. The cosmological parameters are assumed to be H₀ = 70 |$\rm kms^{-1}Mpc^{-1}$|⁠, Ω_M = 0.3, and Ω_Λ = 0.7 based on the lambda-cold dark matter cosmology throughout this paper. All errors that appear in this paper represent 1σ errors.

2 OBSERVATIONS AND DATA REDUCTION

This section explains the data we used. First, we describe the optical/mid-IR photometry data, which mainly correspond to the continuum emission from accretion disc/dust torus, respectively. Secondly, the spectroscopic data mainly used to derive the light curve of Hβ are explained. Finally, we describe the X-ray observations used to investigate the accretion state of the inner part of the accretion disc.

2.1 Optical photometry

We included archival data of photometric observations, because the cadence of spectral observations in this study is insufficient to measure time-lag. The catalogue data used were the Catalina Real-time Transient Survey (CRTS; Drake et al. 2009) and the All Sky Survey Automated Survey for Super-Nova (ASAS-SN; Kochanek et al. 2017). CRTS is unfiltered photometric data, while ASAS-SN observations mainly use g-band and V-band filters. Observations between different filters show consistent patterns of light curves, but a certain amount of offset exists. Therefore, an offset was applied to match the mean values in the overlapping period of the two light curves. Detailed analysis is described in Nagoshi et al. (2021). The essential information to measure the time-lag is the patterns of brightness variation. Therefore, we treat the combined light curve as a continuum light curve of J1258 in this analysis, though it contains values of different filters.

2.2 Mid-infrared (WISE) photometry

The Wide-field Infrared Survey Explorer (WISE; Wright et al. 2010), which is a mid-IR satellite launched in 2009 December and observed the entire sky in four bands (3.4, 4.6, 12, and 22 |$\rm {\mu m}$|⁠), also observed J1258. WISE operations were temporarily terminated in 2011 February, but observations resumed i Decembern 2013 as the Near-Earth Object Wide-field Infrared Survey Explorer (NEOWISE) project (Mainzer et al. 2014). NEOWISE has two bands of 3.4 and 4.6 |$\rm {\mu m}$| and observed the entire sky, visiting each area in the sky about 10 times every six months.

As the information of J1258, we referred to the single-exposure profile-fit magnitude within 3 arcmin from RA = 12:58:09.31 and Dec. = +35:19:43.03. From these values, we selected those with no contamination and confusion flags (‘cc|$\_$|flags’ = 0000). In the present analysis, we used the combined data with each bin accumulated for six months period, sufficient for the variation time-scale of our interest.

2.3 Spectroscopy 1: Sloan Digital Sky Survey

J1258 has been spectroscopically observed twice by the Sloan Digital Sky Survey (SDSS), in 2006 and 2016 April. The SDSS (York et al. 2000) is based on a photometric survey using a 2.5-m wide-field telescope (Gunn et al. 2006) equipped with 30 |$\rm {2k \times 2k}$| CCDs and following spectroscopic observations. Between the two observations, the spectroscopic instrument was upgraded, the wavelength range was extended and the thickness of the fibre on the SDSS plate was changed from 3 to 2 arcsec with the upgrade of the spectrograph.

2.4 Spectroscopy 2: LAMOST

J1258 was observed on 2017 April, by the Large Sky Area Multi-Object Fiber Spectroscopic Telescope (LAMOST, Cui et al. 2012; Zhao et al. 2012), a Schmidt telescope with a primary mirror size of |$6.67 \rm {m} \times 6.05 \rm {m}$|⁠. The wavelength ranges from 370 to 900 nm, and the wavelength resolution is about 1000. The data are reduced with LAMOST pipelines (Luo et al. 2015).

2.5 Spectroscopy 3: MALLS

We performed two spectroscopic observations in 2018 December and 2019 May with the slit spectrograph called ‘MALLS’ (Medium And Low-dispersion Long-slit Spectrograph) installed in the 2-m telescope in Nishi Harima Astronomical Observatory in Japan. The grating was 150 l mm⁻¹, giving a spectral resolution of ∼600, with the 1.2 arcsec slit width. We executed five 1200-s exposures for two nights. The observed data were reduced with the standard processing of slit spectroscopy with iraf (dark subtraction, flat correction, matching sky subtraction, wavelength correction, and flux correction using the spectrophotometric standard).

2.6 Spectroscopy 4: KOOLS-IFU

We performed spectroscopic monitoring observation in 2020 using the low-resolution spectroscopic data from KOOLS-IFU (Kyoto Okayama Optical Low-dispersion Spectrograph with optical-fibre Integral Field Unit; Matsubayashi et al. 2019) installed in the 3.8-m Seimei Telescope (Kurita et al. 2020) at Okayama Observatory of Kyoto University. The data are reduced with KOOLS-IFU pipeline.²

2.7 X-ray

The target was observed twice using the ROSAT (Truemper 1982) Position Sensitive Proportional Counter detector in December 1991 (ObsID = RP600164N00) and December 1992 (RP701178N00) with net exposures of 16.2 and 3.9 ks, respectively. The ROSAT data were processed with Standard Analysis Software System. The source spectra were extracted from a circular region with a radius of 150 arcsec centred on the source by using xselect. The background spectra were selected from a partial annulus region with inner and outer radii of 225 and 450 arcsec by excluding the region of a nearby X-ray source (1RXS J125740.2+351334; Voges et al. 1999). We used the standard response matrix (‘pspcb_gain2_256.rsp’) and generated the ancillary response files using pcarf command in heasoft v6.25.

We also used the data of Swift/X-Ray Telescope (XRT; Burrows et al. 2005) in 2021 February (ObsID = 14079001) with a net exposure of 1.9 ks, covering the ∼0.3–10 keV band. The data reduction was performed with the xrtpipeline v0.13.4. The source spectrum was taken from a circular region of 45 arcsec radius, while the background was from an annulus between 90 and 150 arcsec radii.

3 DATA ANALYSIS AND RESULTS

In this section, we explain the analysis of optical spectra to investigate the variation of |$\rm {H\beta }$|⁠. Then, we measured the time-lag of |$\rm {H\beta }$|⁠, W1, and W2 to the continuum.

3.1 Variation of |$\rm {H\beta }$|

To obtain accurate light curves, we calibrated the spectral flux assuming [O iii]5007 luminosity did not change during all observations. The region where [O iii]5007 is emitted extends over ∼100 pc, so it is often assumed that the intensity remains unchanged over a few decades of observations. We fitted the spectra to determine their [O iii]5007 flux with pyqsofit (Guo, Shen & Wang 2018). The components applied in the fitting are a power-law continuum, three Gaussian components for the broad |$\rm H\beta$|⁠, and one Gaussian each for the narrow |$\rm {H\beta }$|⁠, [O iii]4959/5007. First, the continuum is fitted with power-law using the window of 4450–4500 and 5150–5250 Å in the rest frame, and then the continuum component is subtracted from the spectrum (Fig. 1). We use these continuum-subtracted spectra to fit the emission lines (Fig. 2). The velocity offsets and the line width for [O iii]4959/5007, and the narrow |$\rm {H\beta }$| are set to be equal, and the flux ratio was fixed to 1:3:0.07. Here, the ratio of the narrow component of |$\rm {H\beta }$| was obtained by fitting the first SDSS spectrum, MJD 53848, with free parameters. We fixed the ratio of the narrow emission lines to this value in the fitting of the other spectra. After the fitting, we normalize all flux density using the flux of [O iii]5007 calculated with the Gaussian fitting. The total |$\rm H\beta$| flux is measured by simply integrating the subtracted spectra in 4611–5111 Å after the flux normalization. Each measurement error was estimated using the Monte Carlo method, which consists of 100 fitting calculations repeated with random noise based on the error of the spectrum. The spectra we used is summarized in Appendix A (Fig. A1 and Fig. A2). The value of the light curve of |$\rm H\beta$| is summarized in Table A1 and Fig. 3.

Figure 1.

Example of power-law continuum subtraction. The upper portion displays the original spectrum and the fitted continuum component. The dashed lines indicate the wavelength windows used for the continuum fitting. The lower spectrum line shows the spectrum after subtracting the continuum component. The left panel displays the entire wavelength range, while the right panel focuses on the region around |$\rm {H\beta }$|⁠.

Figure 2.

Example of emission-line decomposition. The solid line presents the emission line spectrum with the continuum subtracted. The dashed line represents the spectrum after subtracting the narrow emission-line components, isolating the broad |$\rm {H\beta }$|⁠. The spectrum of the narrow-line components is illustrated with a dashdotted line. We have assumed that the narrow-line components remain constant throughout the observation period.

Figure 3.

The light curves of continuum (top), |$\rm {H\beta }$| (the second), WISE W1 (the third), and WISE W2 (bottom) of J1258.

We also calculated the mean spectrum and root-mean-square residual spectrum (RMS spectrum) to visualize the typical shape of the emission line and variable component (Peterson et al. 2004). Because of the non-uniform time periods between observations (1 to 3633 d), the calculation is weighted by the length of the observation period as a percentage of the light curve. The definition formula is as follows.

$$\begin{eqnarray} f_{\rm {mean}} &=& \frac{1}{2(t_N-t_1)} \Bigl (f_1 (t_{2}-t_{1}) + f_N (t_{N}-t_{N-1})\\ &+& \sum _{i = 2}^{N-1} f_i (t_{i+1}-t_{i-1}) \Bigr), \end{eqnarray}$$

(2)

$$\begin{eqnarray} f_{\rm {RMS}} &=& \Biggl (\frac{1}{2(t_N-t_1)} \Bigl ((f_1-f_{\rm {mean}})^2 (t_{2}-t_{1})\\ &+& (f_N-f_{\rm {mean}})^2 (t_{N}-t_{N-1}) \\ &+& \sum _{i = 2}^{N-1} (f_i-f_{\rm {mean}})^2 (t_{i+1}-t_{i-1}) \Bigr) \Biggr)^{\frac{1}{2}} . \end{eqnarray}$$

(3)

Here, f represents each spectrum, t represents the date of each observation, and the numbers in the subscripts indicate the order of each spectrum by date of observation.

Here, it is important to note that the RMS spectrum is affected by the observational errors present in each individual spectrum, in addition to the actual fluctuating component of the BLR. Because the measurement error is difficult to estimate accurately in practice, we assumed in the simplest case that the error is constant within the wavelength range of the BLR. In other words, the flux of the RMS (≡ f_RMS(λ)) was assumed to be expressed in the following equation using the actual variable component (≡ f_variable(λ)) and observational error (≡ ϵ).

$$\begin{eqnarray} f_{\rm RMS}(\lambda) = (f_{\rm variable}(\lambda)^{2} + \epsilon ^{2})^{\frac{1}{2}} . \end{eqnarray}$$

(4)

To estimate the constant ϵ, we fitted the RMS spectrum with a Gaussian and ϵ with equation (4) (see the middle panel of Fig. 4). This procedure of removing the measurement error from the RMS spectrum has a significant impact on the calculation of the standard deviation (σ) of the emission line.

Figure 4.

Mean (top), RMS (middle), and the variable component (bottom; baseline subtracted RMS spectrum) spectra of J1258. The middle panel depicts the fitting result for the estimation of ϵ with equation (4). The variable component spectrum shows clearly the broad and symmetric |$\rm {H\beta }$|⁠. The narrow [O iii]4959/5007 lines and asymmetric double-peaked component of broad |$\rm {H}\beta$| appear only in the mean spectrum, except for weak residuals.

These line spectra (mean, RMS, and variable) are all derived from line spectra with the local continuum already subtracted. The FWHM (full width at half-maximum) and σ values have been calculated within the wavelength range of 4600–5120 Å. For the mean spectrum, this calculation was performed after subtracting the narrow-line components. The results of measuring the σ and FWHM of |$\rm {H\beta }$| for the mean spectrum, the RMS spectrum, and the variable component are summarized in Table 1.

Table 1.

Line width of the mean, RMS, and base reduced RMS spectra. The value of the mean spectrum is measured with narrow component subtracted spectrum.

	Mean − narrow	RMS	Variable
FWHM (km s−1)	10268 ± 408	6462 ± 656	5242 ± 556
σ (km s−1)	4727 ± 13	4119 ± 79	3429 ± 157

The mean, RMS, and the variable component spectra for J1258 are shown in Fig. 4. The mean spectrum represents the typical shape of the |$\rm {H\beta }$| and |$\rm {[O\,{\small III}]4959/5007}$|⁠. One of the characteristics of the mean spectrum is asymmetry (the bluer side is stronger), but the variable component spectrum looks symmetric. The RMS spectrum shows that [O iii]4959/5007 profiles are little variable, which indicates that the flux normalization using [O iii]5007 is working well.

3.2 Time-lag measurement

To estimate the radius of the BLR and the dust torus, we calculate the time-lag from the continuum light curve for |$\rm {H\beta }$|⁠, WISE W1, and WISE W2, respectively. We used the results using two different software algorithms to obtain time-lags: interpolated cross-correlation function (ICCF, Peterson 1998; Sun, Grier & Peterson 2018), and Just Another Vehicle for Estimating Lags In Nuclei (JAVELIN, Zu, Kochanek & Peterson 2010, 2011; Zu & Weinberg 2013).

ICCF calculates the cross-correlation function with two linearly interpolated light curves. The time-lag and its error are determined by the bootstrapping method; it calculates the cross-correlation function 10 000 times with light curves of different values within error bars, resulting in the distribution of the correlation peak (i.e. the FR/RSS; see Peterson et al. 1998).

On the other hand, JAVELIN assumes a primary light curve to follow the dumped random walk (DRW) process (e.g. Kelly, Bechtold & Siemiginowska 2009) and convolves it with a top-hat transfer function (TF) to make a responding light curve, searching with the Markov chain Monte Carlo method for the most likely values for the following five parameters: the amplitude and time-scale of the DRW process, height and width of the top-hat TF, and response time delay.

The results of time-lag measurements are summarized in Table 2 and Fig. 5.

Figure 5.

The cross-correlation function and the centroid distribution to the continuum with |$\rm {H\beta }$| (top), WISE W1 (middle), and WISE W2 (bottom) of J1258.

Table 2.

Rest-frame reverberation lags based on ICCF and JAVELIN analysis.

Time-lag target	τrest(ICCF)	τrest(JAVELIN)
|$\rm {H\beta }$|	|$439^{+68}_{-119}$|	|$461^{+60}_{-61}$|
W1	|$805^{+170}_{-137}$|	|$739^{10}_{-291}$|
W2	|$998^{+80}_{-98}$|	|$904^{41}_{-29}$|

4 DISCUSSION

In this section, we first estimate the black hole mass. Secondly, we discuss the mechanism of the extreme brightening event of J1258. Finally, based on a review of the information obtained in this study, we discuss the structure of the BLR, and the origin of the BLR.

4.1 Black hole mass estimation

In order to discuss the structure of the BLR and the dust torus, we first estimate the mass of the central black hole based on the time-lag information of the BLR. The black hole mass is calculated using equation (1), including the object-specific parameters f, τ, and ΔV, which require certain assumptions for each. In this section, we discuss the values of these parameters.

Although the factor f is inherently different for each object because of the influence of BLR structure and viewing angle, previous studies have often applied the same value to different targets as a typical value for many objects. To determine the value of f specifically for each object, a flexible model fitting with respect to the structure of the BLR would be required. However, in this paper, we consider this to be a subject for future work, and for the time being, adopt the value of f = 4.31 ± 1.05 (Grier et al. 2013b) from previous research.

The parameter τ represents the typical radius of the BLR, which is estimated in this paper using two different algorithms, ICCF and JAVELIN. As shown in Table 2, there is no significant discrepancy between the results obtained by the two algorithms. Consequently, accepting the results of either algorithm will not significantly affect the following discussion. However, given reports that JAVELIN tends to estimate time-lags in a way that fits into observational gaps (Shen et al. 2023), here we adopt the results from the ICCF, which involves fewer assumptions.

In this study, we propose a procedure for estimating the ϵ parameter, which contributes to determining a physically plausible velocity width. In equation (1), since τ represents the response time of the fluctuating component of the emission line, it is considered appropriate for the velocity ΔV of the emission line to be measured in the variable component of the spectrum. As the RMS spectrum is affected by substantial observational errors, we contend that it is physically more appropriate to estimate the velocity using the variable component, the error-subtracted RMS spectrum, for the calculation of the black hole mass.

In summary, this paper refers the f value from Grier et al. (2013b), and adopts the ICCF estimate for |$\tau (439^{+68}_{-119} \ {\rm d})$|⁠, and σ(3429 ± 157 km s⁻¹) for the variable component of ΔV. It is important to note that regardless of the method employed, there is no difference in the order of magnitude, and no significant impact on the discussion of BLR and dust torus structure in the subsequent sections. In the following analysis, we use |$10^{9.64^{+0.11}_{-0.20}}\rm {M_\odot }$| for the black hole mass calculated using equation (1).

4.2 Mechanism of the extreme variability

In general, the leading cause of the extreme variation accompanied by CSQs is considered to be intrinsic variability in the mass accretion rate of the disc for the following reasons. First, the variability time-scale of BELs and continuum is shorter than the orbital time-scale of dust torus (>∼100 yr). Second, the reddening in optical (e.g. MacLeod et al. 2016) and the X-ray hydrogen column density (Husemann et al. 2016) did not change before and after the state transition. Third, the brightness in the mid-IR changed in response to continuum emission (Stern et al. 2018; Nagoshi et al. 2021). Finally, polarization was not detected in spectropolarimetry (Hutsemékers et al. 2017, 2019). The above results support that CSQs are not explained by any shielding theories, but are consistent with the intrinsic variation of mass accretion rate.

In Nagoshi et al. (2021), we proposed a hypothesis that the J1258’s brightening event was caused by a state transition in the accretion disc, based on its time-scale. To reconsider the hypothesis, we performed spectroscopic monitoring and X-ray observations to estimate the optical-to-X-ray spectral indices (α_ox) in the on/off state each.

4.2.1 α_ox and Eddington ratio

We analysed the X-ray spectra of the ROSAT and Swift/XRT data with xspec v12.10.1 (Arnaud 1996). We consider Galactic absorption by the phabs model and refer the hydrogen column density fixed at the value of Willingale et al. (2013). Abundances were assumed at solar values. Cash (1979) statistics were utilized to fit ROSAT (0.1–2.4 keV) and Swift (0.3–10 keV) spectra, and the spectra were binned to have one count per bin. The spectra were well reproduced by a single power-law model (zpowerlw) with Galactic absorption, as shown in Fig. 6. The best-fitting parameters of observed-frame 2 keV (⁠|$\nu F_{\rm 2\, keV}$|⁠) and 2–10 keV (⁠|$F_{\rm 2-10\, keV}$|⁠) fluxes, and rest-frame 2 keV (⁠|$\nu L_{\rm 2\, keV}$|⁠) and 2–10 keV (⁠|$L_{\rm 2-10\, keV}$|⁠) luminosities,³ all of which are corrected for Galactic absorption, are presented in Table 3. The photon indices are ∼1.7–1.9, consistent with a typical value in active galactic nuclei (AGNs, e.g. Ueda et al. 2014; Ricci et al. 2017). We find that the rest-frame 2 keV luminosities increase from ∼1 × 10⁴⁴ erg s⁻¹ in 1991 and 1992 to ∼4 × 10⁴⁴ erg s⁻¹ in 2021.

Figure 6.

The unfolded X-ray spectra and best-fitting models of the observations with ROSAT in 1991 (black), 1992 (red), and Swift/XRT in 2021 (green). The bottom panel shows the ratios between the data and models.

Table 3.

Best-fitting X-ray spectral parameters, observed-frame absorption-corrected fluxes, rest-frame luminosities, and reduced C statistics.

Parameters	ROSAT	ROSAT	Swift
	(1991)	(1992)	(2021)
Photon index	|$1.82^{+0.06}_{-0.05}$|	|$1.94^{+0.07}_{-0.08}$|	|$1.73^{+0.10}_{-0.10}$|
|$\nu F_{\rm 2\, keV}$| (10−13 erg s−1 cm−2)	|$3.25^{+0.29}_{-0.28}$|	|$4.44^{+0.54}_{-0.49}$|	|$14.07^{+1.18}_{-1.13}$|
|$F_{\rm 2-10\, keV}$| (10−13 erg s−1 cm−2)	|$6.03^{+0.54}_{-0.51}$|	|$7.52^{+0.91}_{-0.83}$|	|$28.27^{+2.37}_{-2.27}$|
|$\nu L_{\rm 2\, keV}$| (1044 erg s−1)	|$0.96^{+0.07}_{-0.07}$|	|$1.36^{+0.13}_{-0.13}$|	|$4.07^{+0.29}_{-0.29}$|
|$L_{\rm 2-10\, keV}$| (1044 erg s−1)	|$1.78^{+0.14}_{-0.13}$|	|$2.29^{+0.24}_{-0.22}$|	|$8.16^{+0.59}_{-0.57}$|
C/dof	166.8/174	147.2/140	147.8/192

We compare the Eddington ratio and α_ox of before and after the state transition with the similar method of Jin et al. (2021). α_ox is a spectral index between optical and X-ray characterizing the state of the accretion disc. We use the following definition to calculate α_ox (Tananbaum et al. 1979),

$$\begin{eqnarray} \alpha _{\rm {ox}} = - \frac{{\rm log}(\lambda L_{2500\, {\mathring{\rm A}}})-{\rm log}(\nu L_{\rm {2\, keV}})}{{\rm log}\nu _{2500\, {\mathring{\rm A}}}-{\rm log}\nu _{\rm {2\, keV}}} + 1 \end{eqnarray}$$

(5)

in which the continuum luminosity at 2500 Å is estimated by extrapolation of the fitted power-law component. For the spectral fitting, the power-law component, Fe ii template (Boroson & Green 1992), and the Balmer continuum component were applied to 4450–4500 and 5150–5250 Å avoiding narrow emission lines and BELs. The errors were estimated by Monte Carlo simulation. The bolometric luminosity is estimated using the correction formula in Lusso et al. (2010).

$$\begin{eqnarray} {\rm log}(L_{\rm {bol}}) = {\rm log}(L_{2-10\, \rm {keV}}) + 15.61 - 1.853\alpha _{\rm {ox}} + 1.226\alpha _{\rm {ox}}^2 . \end{eqnarray}$$

(6)

The obtained α_ox, bolometric luminosity, and Eddington’s ratio before and after the state transitions are summarized in Table 4.

Table 4.

Table summarizing the parameters representing the accretion disc state before and after state transitions. Each parameter is defined as follows: |${\rm log}(\lambda L_{2500\, {\mathring{\rm A}}})$| represents the λL at 2500 Å, L_bol is the estimated bolometric luminosity using equation (6), the Eddington ratio is calculated from L_bol and the black hole mass, and α_ox is calculated using equation (5).

Parameters	Off state	On state
log(λL2500Å)	44.289 ± 0.010	45.337 ± 0.001
αox	1.02 ± 0.07	1.29 ± 0.05
logLbol	45.3 ± 0.22	46.1 ± 0.11
Eddington ratio	0.004 ± 0.009	0.024 ± 0.051

We used optical spectra and X-ray data taken during off and on states each. Since the optical light curve in Fig. 7 shows that the state transition happened in 2007 (MJD around 54000) to 2015 (MJD around 57000), we consider the status before 2007 as off state, while the status after 2015 as on state.

Figure 7.

The optical light curve with periods of optical spectra and X-ray observation. The black dots represent the optical light curve, and the vertical lines drawn with dotted lines represent the time of spectral or X-ray observations. In addition to the CRTS and ASAS-SN light curves, data from The Guide Star Catalogue (Lasker et al. 2008) are plotted. We interpret the data before MJD = 54 000 are considered as off state, and the data after MJD = 57 000 as on state.

Jin et al. (2021) compared the α_ox of 10 previously discovered CSQs with the simulated results from Sobolewska, Siemiginowska & Gierliński (2011) based on X-ray observations of the X-ray binary GRO J1655 − 40 extended to AGNs. They confirmed that it is consistent with the prediction of the state transition model. Jin et al. (2021) found the negative correlation between Eddington ratio and α_ox for the objects with low (< 0.01) Eddington ratio, while Lusso et al. (2010) reported the positive correlation for the objects with high (> 0.01) Eddington ratio.

4.2.2 Accretion rate variation

The measured Eddington ratio and α_ox of J1258 in the off and on states each are consistent with the state transition model suggested in the previous research (Noda & Done 2018). We show the comparison of the Eddington ratio and α_ox on Fig. 8 together with the best-fitting results of Jin et al. (2021) and Lusso et al. (2010).

Figure 8.

Comparison of the measured α_ox and λ_Edd for J1258 in its on and off states. the dashed line and surrounding region indicate the fitting results and associated errors for the off state derived from quasars and CSQs, as referenced by Jin et al. (2021). In contrast, the solid line and surrounding the region represent the α_ox − λ_Edd relationship and its associated errors for the on state, based on findings from the XMM the Cosmic Evolution Survey (Lusso et al. 2010).

One of the possible causes of the state transition is thermal disc instability (Noda & Done 2018). According to the disc instability model, the accretion disc forms a geometrically thin disc at high Eddington ratios (Shakura & Sunyaev 1973) characterized by soft excess in X-ray, whose accretion state is called the high/soft state. During the low Eddington ratio, the inner region of the thin disc is evaporating (Esin, McClintock & Narayan 1997), resulting in a hot and radiatively inefficient accretion flow (Shapiro, Lightman & Eardley 1976), possibly advection dominated (Narayan & Yi 1994). In this state, the X-ray spectrum is dominated by Compton scattering of hard X-rays originating from the accretion flow, called the low/hard state. The spectral index from optical to soft X-ray component, α_ox, is used to characterize these spectral states. The correlation between α_ox and the Eddington ratio is predicted theoretically to change according to certain rules (negatively correlated to the Eddington ratio during low/hard state, while positively correlated during high/soft state) during the state transition, and has been confirmed observationally in X-ray binaries (e.g. Maccarone & Coppi 2003; Debnath, Chakrabarti & Nandi 2010).

Another possibility is that the state transition was triggered by external factors such as the tidal disruption of a star (Merloni et al. 2015; Ricci et al. 2020). The outburst caused by tidal disruption events (TDEs) is expected to show a decline in flux with t^−5/3 (van Velzen et al. 2020). In the case of J1258, however, it shows brightening again from around MJD = 58 000 (Fig. 7), which is inconsistent with the expected patterns of the light curve of TDEs.

One other possibility is that J1258 is a pair of SMBHs in a binary system. Asymmetric BELs are often referred to as one of the characteristics of binary SMBHs (Eracleous et al. 2012; Runnoe et al. 2015). That is a situation where two SMBHs each form an AGN with their own BLRs. Given this situation, if a state transition occurs in one of the AGNs there is no need for occurring the state transition in the other AGN. Furthermore, if the interaction of the two SMBHs efficiently extracted angular momentum, it could be considered a trigger for the large-scale brightening event (Wang & Bon 2020). To confirm this hypothesis, it would be helpful to test whether the velocity shift in the asymmetric component of the BELs changes according to the orbital motion. However, it is not easy to calculate the velocity shift of BELs independently of varying intensities because the quality of the obtained spectra in this study is not high enough. If a velocity shift were to present, it would be valid to compare the spectra again after sufficient time has passed, which could be verified in future studies.

We conclude here the main cause of the J1258 brightening event is a variation of the intrinsic accretion power caused by the state transition of the accretion disc, as we mentioned in Nagoshi et al. (2021).

4.3 Structure of the BLR

We discuss the structure and the origin of the BLR using the extreme variation of the BEL of J1258. The spectra in Fig. 4 show the mean, the RMS, and the variable component of |$\rm H\beta$|⁠, respectively. Examining Fig. 4, we can see that the shapes of the mean and variable spectra are very different. To show the differences better, we depict the mean spectrum with subtracted narrow lines in Fig. 9. This figure also shows how the BEL shape would look like if there were greater changes in the variable component. Here, we simply assume that the line shape changes with a constant multiple of the variable component in all velocities. The spectrum on the top of Fig. 9 represents the mean spectrum. The second spectrum from the top is the mean spectrum that the narrow lines (the narrow components in |$\rm {H\beta }$| and [O iii]4959/5007) are subtracted. The third and fourth spectra from the top are 1 or 3 times of the variable component subtracted from the second spectrum (the narrow lines subtracted mean spectrum). We can see that the spectral shape at the bottom of Fig. 9 shows a double-peaked profile. On the other hand, the variable component has a relatively symmetrical, single-peaked shape. These represent that the spectra typically have a double-peaked component, but the single-peaked component changes with continuum variation, indicating that there are two distinct regions in the BLR (Hu et al. 2020).

Figure 9.

Diagram visualizing the double-peaked component on the BEL. The top spectrum is the mean spectrum, the second from the top is the mean spectrum minus the narrow emission line, and the third is the second spectrum minus the variable component. The bottom spectrum is the mean spectrum minus three times the variable component. The baseline of each spectrum has been added to a constant for clarity.

Indeed, since this object was observed during a relatively faint period, the double-peaked shape can be seen in real data. Fig. 10 compares the |$\rm H\beta$| shape by subtracting the continuum and narrow emission lines (⁠|$\rm H\beta$| and [O iii]4959/5007) based on the spectra from SDSS. The spectrum observed in 2006 has a shape with a small peak at the central wavelength as well as peaks on each side. Considering the central component as variable, we can see a double-peaked component less responsive to changes in the continuum luminosity. Using J1258’s extreme luminosity variation, we successfully clarified the existence of components with different characteristics in the BEL.

Figure 10.

Narrow lines subtracted |$\rm H\beta$| from SDSS. The spectrum with dashdotted line was taken on 2006 April 23, and the spectrum with solid line was taken on 2016 April 03.

Based on the above spectral components, the BLR of J1258 is expected to consist of two different regions emitting a single-peaked component (RMS spectrum) and a double-peaked component (the bottom spectrum in Fig. 9), respectively. Here we name these components Region X and Region Y, respectively:

• Region X: emitting a relatively slow (∼3000 km s⁻¹), single-peaked line that responds to the state transition of the accretion disc.

• Region Y: emitting a relatively fast (∼7,000 km s⁻¹), double-peaked line, which is not likely affected by the state transition of the accretion disc.

Region X has the same characteristics as the BLR seen in typical quasars, whose size can be estimated from the measured time-lag of the reverberation observation. The BLR size of Region X was found to be 0.37 pc (439 light days) from the centre (Section 3).

The kinematics responsible for the double-peaked emission lines observed in Region Y could potentially be attributed to fast outflow, fast inflow, or rotating ring-like structures. While strong outflows and inflows are typically influenced by the mass accretion rate, as suggested in studies such as Elitzur et al. (2014), our observations indicate that the double-peaked component of Region Y did not exhibit significant variations during state transitions (see Fig. 4). This observation suggests that these emissions might occur independently of the accretion disc’s activities, implying a more stable and passive kinematic process. Although the lack of response to continuum variation does not entirely rule out the possibility of outflows or inflows, it encourages consideration of a scenario where these processes are steady and unrelated to the accretion disc’s dynamics. While acknowledging that outflow and inflow cannot be completely dismissed, the evidence aligns more consistently with a rotating ring-like structure as the source of the double-peaked broad |$\rm H\beta$| emission, as discussed in previous research (e.g. Eracleous & Halpern 1994, 2003). This interpretation is considered more plausible in our context, given the observed stability of the emission lines, though further investigation is warranted to elucidate these kinematics fully.

With the above points taken into account, it is possible to combine them into a single figure without any geometrical inconsistency. The distance to Region X can be determined from the time-lag, which is about 0.37 pc. The distance to Region Y is about |$0.29\ \rm {pc} \sim 1400\ \rm {R_g} (\equiv \frac{GM}{c^2})$|⁠, which can be explained by the rotation of the disc, given that the black hole mass is |$10^{9.64^{+0.11}_{-0.20}}\rm {M_\odot }$| and the typical speed is about 8000 km s⁻¹ (the peak velocity of the double-peaked component). If the asymmetric shape is caused by a bias in the distribution of ionized gas in Region Y, the period over which the gas makes one cycle with Kepler motion is about 225 yr. Therefore, detecting a velocity shift over a 20-yr observation period would be difficult. The distance to the dust torus is estimated by the response time of W1 to be 0.68 pc (805 light days). These locations are shown in Fig. 11.

Figure 11.

Conceptual view of the central core structure of J1258 in On state inferred from the time-lag and wavelength shift. The distance to Region X and the black hole mass are estimated from the |$\rm {H\beta }$| time-lag. The distance to Region Y is estimated from the determined black hole mass, and the distance to the inner edge of the dust torus is estimated from the WISE W1 time-lag.

Assuming the structure of Fig. 11, it is also possible to explain why Region Y was not affected by the state transition of the accretion disc. The brightening phenomenon due to the state transition from the low/hard state to the high/soft state is caused by the variation of soft X-rays from the inner part of the accretion disc being larger than that of hard X-rays from the hot corona (⁠|$T \sim 10^{9} \ \rm {K}$|⁠) region (Noda & Done 2018). That is, the intensity of the hard X-ray itself does not change much before and after the state transition, which is also confirmed observationally (e.g. Mathur et al. 2018). The Compton scattering component from the coronal region ranges from X-ray to the UV which contributes to the ionization of the Balmer lines. Suppose that Region Y is located right next to the disc so that only the axially extended hot corona region can supply ionizing photons to Region Y. In that case, this is a possible explanation for why the intensity of the emission lines emitted from Region Y is less affected by the state transition.

We confirmed that the double-peaked feature of Region Y can be reproduced by the structure shown in Fig. 11 with model calculation. We adopt a model of Chen, Halpern & Filippenko (1989) here. This model can calculate the BEL spectrum, including relativistic effects, for a BLR distributed as a thin Keplerian disc. This model is defined by five parameters; the viewing angle i, the inner edge of the BLR ξ₁, the outer edge of the BLR ξ₂, the emissivity index q, and the local broadening parameter σ_l. The model assumes that the surface intensity distribution of the disc is determined by equation (7), which is a function of ξ₁ and ξ₂ and q, and that locally emission lines are emitted with a velocity component of Keplerian rotation and a velocity dispersion of σ_l. The surface intensity ϵ(ξ) is defined as follows:

$$\begin{eqnarray} \epsilon (\xi) = \frac{\epsilon _{0}}{4\pi }\xi ^{-q} {\rm \ [ergs \ cm^{-2} \ s^{-1} \ sr^{-1}]}, \end{eqnarray}$$

(7)

where, ϵ₀ is a peak value of the emissivity and ξ is radius from the centre of the disc. We show here as an example of parameters that can reproduce the qualitative features of Region Y, the results calculated with i = 47^○, ξ1 = 150R_s (R_s is Schwarzschild radius; 150R_s = 0.04 pc), ξ₂ = 900R_s (= 0.25 pc), q = 0.9, and σ_l = 1700 km s⁻¹ are shown in Fig. 12, where the spectral shape obtained by subtracting 2.4 times the variable component from the mean spectrum is well reproduced. This result shows that the observed spectral shape of the BEL can be qualitatively reproduced by considering a disc-shaped BLR (Region Y) with a thin hole inside the BLR radius (at Region X) estimated by the time delay, as shown in Fig. 11.

Figure 12.

Example of reproduction of Region Y component spectrum by a model of Chen et al. (1989). The solid line is the mean spectrum, which is decomposed into a narrow emission line (dashed-dotted line), a constant multiple of the variable component (dotted line), and Region Y component (solid line). The dashed line shows an example of fitting the Region Y component with the model of Chen et al. (1989).

Additionally, utilizing the parameters presented in Fig. 12, we verified that by simply optimizing the magnification factors of the components from Region X (the variable component) and Region Y (based on the model of Chen et al. 1989), the Hβ profile of J1258 for all observed periods can be accurately reproduced. Fig. 13 displays the fitting results for the dimmest, brightest, and intermediate periods, chosen from spectra with high signal-to-noise ratios. This outcome indicates that our model is effectively able to explain the observed phenomena.

Figure 13.

Fitting results of the Hβ emission line using two components: the variable component as Region X and the model of Chen et al. (1989) as Region Y. Each panel, from top to bottom, displays the Hβ and [O iii] emission lines with the continuum subtracted for MJD=53848, 57481, and the mean of 58953, 58963, 58964, 58982, and 58991, respectively. In each panel, the upper section illustrates the fitting results of each component, while the lower section presents the fitting residuals. The black line indicates the observed flux, the blue line denotes the narrow component, the green line corresponds to the Region Y component, the purple line to the Region X component, and the red line represents the sum of all model components. All panels are displayed with the same grid size.

4.4 Origin of the two components of the BLR

Based on the structure shown in Fig. 11, we discuss the origins of the gas in the BLR from two major theories: one is that the gas is supplied from disc winds from accretion discs (e.g. Gaskell 1982; Nicastro 2000; Nicastro, Martocchia & Matt 2003; Czerny & Hryniewicz 2011) and the other is that it is derived from accreted components from the dust torus (e.g. Gaskell 2009, 2010; Hopkins et al. 2012; Gaskell & Goosmann 2013; Wang et al. 2017). The disc-wind origin theory is based on a model in which outflow generated from the accretion disc according to an increase of the mass accretion rate winds up the gas, which is ionized by the emission from the inner accretion disc. This model suggests that no BLR is produced below a certain threshold because the disc wind depends on the mass accretion rate (Elitzur et al. 2014). This dependency is sometimes used as the reason for the appearance/disappearance of the BLR with changes in luminosity in the CSQs (e.g. LaMassa et al. 2015). On the other hand, the dust torus origin model assumes that BLR is continuously connected with the dust torus, where the ionized component becomes the BLR during mass accretion to the accretion disc.

We suggest that the results of this study can be used to validate the model of the origin of the BLR observationally. First, assuming that disc winds from the accretion disc are the origin of the BLR, strong outflow is expected. Features of outflow are detectable by velocity-resolved RM, which predicts that the red side of the BEL follows the blue side. In other words, the disc wind model is not dominant if the velocity-resolved RM does not detect outflow. The independent variation of the BEL from the continuum variation also excludes the disc wind model, because these two should be correlated in this model. We will now discuss the origin of the Region X/Y in terms of two aspects: velocity-resolved RM and response to continuum variation, respectively.

To confirm the origin of Region X, we conducted a simple velocity-resolved RM. We calculated the cross-correlation function for each of the light curves of the red wing (4701.3–4861.3 Å) and the blue wing (4861.3–5021.3 Å) components of |$\rm {H\beta }$|⁠. As a result, the measured time-lag is |$4.1^{+3.5}_{-3.1}$| d (Fig. 14). If the gas in Region X was come from disc winds, the outflow reach Region X in 0.37 pc from the centre in 20 yr, which means the velocity is about 1  per cent of the light speed. When there is strong outflow, the red wing should lag behind the blue wing, because the part far from the observer is moving away (Gaskell 2009). Here, it is unlikely that gas in Region X was supplied by outflow, as no significant time-lags could be measured for the red and blue wings. In addition, the evidence of strong inflow was not observed either, suggesting that gas in Region X does not have significant radial motion.

Figure 14.

ICCF of blue wings relative to red wings of |$\rm {H\beta }$|⁠, and the distribution of ICCF centre.

Because no significant radial motion was detected in Region X, we can consider that the gas existed near the torus before being ionized by the increase in central luminosity (see also Noda et al. 2023). The radius of the ionizing region is thought to be proportional to the central luminosity to the power of 1/2; the BLR’s ‘breathing’ phenomena have been reported for several objects (e.g. Wang et al. 2020b). Indeed, long-term optical/near-infrared observation reports that the inner radius of the torus of a CSQ, Mrk 590, decreases rapidly following a reduction in the AGN luminosity (Kokubo & Minezaki 2020). The dust sublimation radius measured by the time-lag of W1 is about twice the radius of Region X. The time-lag of W1 is mainly dominated by the time-lag when J1258 turns to dim, which may reflect the sublimation radius after it brightens. It is plausible to assume that the dust sublimation radius increased by a factor of 2 during the RM observation period since the continuum luminosity increased 4–5 times during the observation period. Therefore, we conclude that the origin of the gas in Region X can be said to be the sublimated torus, which was ionized passively with the luminosity variation.

We note, however, that the time-lag measurements in this study may underestimate the inner edge radius of the torus. Examining the JAVELIN light curves for the estimated time-lag between the continuum and the W1 and W2 bands (see Fig. A3), the predicted increase/decrease pattern and the actual observation are different, especially in the W2 band (the JAVELIN light curve shows a decline around MJD=58500 d, while the actual observation in the W2 band continues to increase). This potential underestimation of the time-lag may result from the limited observation period. However, even if the actual inner radius of the dust torus measured from the time-lag is 2–3 times larger than that in this paper, the inner radius of the dust torus can be expected to have also increased 6–7 times since the object has brightened approximately 40 times over 40 yr. Therefore, the gas comprising the BLR observed today is likely the former dust torus.

Region Y is unlikely to have been formed by outflow, as it is not affected by variation in central luminosity. The asymmetric double-peak component suggests that the gas distribution is biased, formed within a relatively short time-scale (< gas circling time-scale: ∼200 yr). We hypothesize that the dust clumps have been tidally disrupted to create an inhomogeneous gas distribution over a short period (Wang et al. 2017). In other words, Region Y was formed independently of the state of the accretion disc and thus is less affected by the state transitions.

In summary, the BLR gas in J1258 has two very different positional regions, both of which are likely sourced from the dust torus. Although both originate from the dust torus, we conclude that they have different mechanisms of their origin. One is that a large variation due to state transition changes the radius of being ionized to extend the BLR (Region X). The other is the ionization of components that enter inside the sublimation radius by accretion of the dust torus (Region Y).

5 CONCLUSIONS

In this research, we conducted optical and mid-IR RM and X-ray follow-up observation of an extremely variable CSQ; J1258. This is the first result of RM of the most extremely variable quasar. Our conclusions are as follows.

• Black hole mass of J1258 is |$10^{9.64^{+0.11}_{-0.20}}\rm {M_\odot }$|⁠.

• The brightening event can be explained by the state transition of the accretion disc.

• The location of the BLR near the central core is inferred from the RM results (Fig. 11). There are two components in the BLR, each located at 0.29 and 0.37 pc from the centre.

• The origin of the ionized gas in Region X in the BLR is likely to be the sublimated dust torus. The origin of the ionized gas in Region Y is likely to be accretion from dust torus.

ACKNOWLEDGEMENTS

We express our gratitude to the anonymous reviewer for their insightful comments and valuable suggestions, which have significantly contributed to enhancing the quality of this paper.

This work is financially supported by JSPS KAKENHI grant nos 22J13428 (SN); 22K20391 and 23K13154 (SY). SY is grateful for support from RIKEN Special Postdoctoral Researcher Program. MO was supported by JSPS Grants-in-Aids for Scientific Research(C) (JP19K03914 and 22K03675). This work is partly supported by JSPS KAKENHI grant nos JP18H05439 and JP20K14521.

We would like to express our gratitude to Dr Miho Kawabata for support during the observations.

Funding for the Sloan Digital Sky Survey-IV has been provided by the Alfred P. Sloan Foundation, the U.S. Department of Energy Office of Science, and the Participating Institutions.

SDSS-IV acknowledges support and resources from the Center for High Performance Computing at the University of Utah. The SDSS website is www.sdss4.org.

SDSS-IV is managed by the Astrophysical Research Consortium for the Participating Institutions of the SDSS Collaboration including the Brazilian Participation Group, the Carnegie Institution for Science, Carnegie Mellon University, Center for Astrophysics| Harvard & Smithsonian, the Chilean Participation Group, the French Participation Group, Instituto de Astrofísica de Canarias, The Johns Hopkins University, Kavli Institute for the Physics and Mathematics of the Universe (IPMU)/University of Tokyo, the Korean Participation Group, Lawrence Berkeley National Laboratory, Leibniz Institut für Astrophysik Potsdam (AIP), Max-Planck-Institut für Astronomie (MPIA Heidelberg), Max-Planck-Institut für Astrophysik (MPA Garching), Max-Planck-Institut für Extraterrestrische Physik (MPE), National Astronomical Observatories of China, New Mexico State University, New York University, University of Notre Dame, Observatário Nacional/MCTI, The Ohio State University, Pennsylvania State University, Shanghai Astronomical Observatory, United Kingdom Participation Group, Universidad Nacional Autónoma de México, University of Arizona, University of Colorado Boulder, University of Oxford, University of Portsmouth, University of Utah, University of Virginia, University of Washington, University of Wisconsin, Vanderbilt University, and Yale University.

Guoshoujing Telescope (the Large Sky Area Multi-Object Fiber Spectroscopic Telescope, LAMOST) is a National Major Scientific Project built by the Chinese Academy of Sciences. Funding for the project has been provided by the National Development and Reform Commission. LAMOST is operated and managed by the National Astronomical Observatories, Chinese Academy of Sciences.

This work made use of data supplied by the UK Swift Science Data Centre at the University of Leicester.

The CSS survey is funded by the National Aeronautics and Space Administration under grant no. NNG05GF22G issued through the Science Mission Directorate Near-Earth Objects Observations Program. The CRTS survey is supported by the U.S. National Science Foundation under grants AST-0909182.

DATA AVAILABILITY

The data underlying this article will be shared on reasonable request to the corresponding author.

Footnotes

1

Strictly speaking, the definitions of changing-look quasars and changing-state quasars are different (see a review of Ricci & Trakhtenbrot 2023). However, for the purpose of this paper, they are treated similarly.

2

We used the pipeline here http://www.kusastro.kyoto-u.ac.jp/~iwamuro/KOOLS/index.html

3

The uncertainties of |$\nu F_{\rm 2\, keV}$| and |$\nu L_{\rm 2\, keV}$| are obtained by the cflux and clumin models in xspec. The uncertainties of the 2–10 keV fluxes and luminosities are assumed to the log-scale uncertainties at the 2 keV band.

Appendix A:

In Fig. A1, we show all spectra used in this study. In Table A1, we summarize the information of each spectrum and the result of the spectral fitting. In Fig. A3, we show the estimated light curves produced by JAVELIN.

Figure A1.

Spectra used in this study. For better visibility, constants have been added to the baseline for each observation. The colours correspond to the date of observation.

Figure A2.

Spectra around the |$\rm H\beta$|⁠, with the continuum of each spectrum subtracted. For better visibility, constants have been added to the baseline for each observation. The colours correspond to the date of observation.

Figure A3.

The light curves of continuum (top), |$\rm {H\beta }$| (the second), WISE W1 (the third), and WISE W2 (bottom) of J1258.

Table A1.

List of the spectroscopic observations we use here.

Observation date	Instrument	Exposure time	Standard star	λL5100Å	|$\rm {H\beta }$|
MJD (MM-DD-YYYY)		second × number		1044 erg s−1	1042|$\rm {erg\,s^{-1}}$|
53848 (04-23-2006)	SDSS	3000 × 3	−	2.41 ± 0.02	7.4 ± 1.5
57481 (04-03-2016)	SDSS	3600 × 4	−	10.37 ± 0.02	19.7 ± 2.7
57846 (04-03-2017)	LAMOST	5400 × 3	−	8.69 ± 0.02	20.6 ± 2.2
58470 (12-18-2018)	MALLS	1200 × 5	HZ44	6.67 ± 0.09	11.6 ± 6.1
58606 (05-03-2019)	MALLS	1200 × 5	BD+33d2642	7.63 ± 0.03	12.3 ± 2.6
58607 (05-04-2019)	MALLS	1200 × 4	BD+33d2642	7.22 ± 0.05	12.1 ± 3.9
58612 (05-09-2019)	MALLS	1200 × 4	HZ44	7.44 ± 0.05	12.4 ± 3.3
58899 (02-20-2020)	KOOLS-IFU	600 × 10	HZ44	8.06 ± 0.03	12.6 ± 2.3
58913 (03-05-2020)	KOOLS-IFU	600 × 6	HZ44	7.66 ± 0.05	10.2 ± 3.6
58927 (03-19-2020)	KOOLS-IFU	600 × 8	HZ44	9.60 ± 0.03	13.4 ± 2.4
58940 (04-01-2020)	KOOLS-IFU	600 × 8	HZ44	9.87 ± 0.03	13.2 ± 3.7
58953 (04-14-2020)	KOOLS-IFU	600 × 8	HZ44	11.05 ± 0.03	14.6 ± 1.7
58963 (04-24-2020)	KOOLS-IFU	600 × 7	HZ44	11.29 ± 0.03	14.9 ± 2.6
58964 (04-25-2020)	KOOLS-IFU	600 × 6	HZ44	12.70 ± 0.03	16.4 ± 3.1
58967 (04-28-2020)	KOOLS-IFU	600 × 8	HZ44	11.30 ± 0.03	16.0 ± 2.2
58982 (05-13-2020)	KOOLS-IFU	600 × 8	HZ44	11.21 ± 0.03	16.6 ± 2.8
58989 (05-20-2020)	KOOLS-IFU	600 × 4	HZ44	11.77 ± 0.05	15.2 ± 3.5
58991 (05-22-2020)	KOOLS-IFU	600 × 7	HZ44	12.09 ± 0.03	14.5 ± 2.5