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

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 disk 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$\beta$ reverberation mapping indicated a black hole mass of $10^{9.64^{+0.11}_{-0.20}}\rm{M_\odot}$. Additionally, we analyzed variations in the optical to X-ray spectral indices, $\alpha_{\rm{ox}}$, before and after the state transition, to investigate the accretion disk. These variations in $\alpha_{\rm{ox}}$ and the Eddington ratio (from 0.4 \% to 2.4 \%) exhibitied behavior similar to state transitions observed in X-ray binary systems. Spectral analysis of H$\beta$ revealed a predominantly double-peaked profile during dim periods, transitioning to include a single-peaked component as the quasar brightened, suggesting that H$\beta$ 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 disk, and the dust torus. Our results suggest that the BLR consists of two distinct components, each differing in location and origin.


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 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 con-★ E-mail: shumpei@kusastro.kyoto-u.ac.jp stant ( 0 ) with cosmic microwave background (Planck Collaboration et al. 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  0 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 disk 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 et al. 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 disk 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. 2000Kaspi et al. , 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; where  is the speed of light,  is the gravitational constant, and Δ is the velocity of the BLR gas.The  is a factor that reflects the structure, kinematics, and inclination of the BLR.The problem is that the  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 disks as the disk-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 disk wind from the accretion disk, it is expected that the amount of disk wind changes with the mass accretion rate of the accretion disk, resulting in a variation in the BEL intensity (Elitzur et al. 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. 2019Lu et al. , 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) Quasars1 (CSQs).We focus on the variability of BEL and continuum luminosity, which show significant variations (> 30%) in the BELs and in mid-infrared luminosity (LaMassa et al. 2015;MacLeod et al. 2016MacLeod et al. , 2019;;Ruan et al. 2016;Gezari et al. 2017;Wang et al. 2018;Noda & Done 2018;Ross et al. 2020;Graham 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 disk also provides essential information, the observational studies are still not enough for quasars because of their long variation timescale.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 disk.
In this paper, we study the BLR structure of SDSS J125809.31+351943.0 (hereafter J1258;  = 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 magnitudes in about 30 years).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  0 = 70 kms −1 Mpc −1 , Ω  = 0.3 and Ω Λ = 0.7 based on the ΛCDM cosmology throughout this paper.All errors that appear in this paper represent 1 errors.

OBSERVATIONS AND DATA REDUCTION
This section explains the data we used.Firstly, we describe the optical/mid-infrared photometry data, which mainly correspond to the continuum emission from accretion disk/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 disk.

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 catalog data used were the Catalina Realtime 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 -band and -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.

Mid Infrared (WISE) photometry
The Wide-field Infrared Survey Explorer (WISE; Wright et al. 2010), which is a mid-infrared satellite launched in December 2009 and observed the entire sky in four bands (3.4, 4.6, 12, and 22 𝜇m), also observed J1258.WISE operations were temporarily terminated in February 2011, but observations resumed in December 2013 as the NEOWISE project (Mainzer et al. 2014).NEOWISE has two bands of 3.4 and 4.6 m and observed the entire sky, visiting each area in the sky about ten times every six months.
As the information of J1258, we referred to the single-exposure profile-fit magnitude within 3" from RA = 12:58:09.31and 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 timescale of our interest.

Spectroscopy 1: Sloan Digital Sky Survey
J1258 has been spectroscopically observed twice by the SDSS, in April 2006 and April 2016.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 2k × 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 fiber on the SDSS plate was changed from 3 ′′ to 2 ′′ with the upgrade of the spectrograph.

Spectroscopy 2: LAMOST
J1258 was observed on April 2017, 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.67m × 6.05m.The wavelength ranges from 370 nm to 900 nm, and the wavelength resolution is about 1000.The data is reduced with LAMOST pipelines (Luo et al. 2015).

Spectroscopy 3: MALLS
We performed two spectroscopic observations in December 2018, and May 2019 with the slit spectrograph called "MALLS" (Medium And Low-dispersion Long-slit Spectrograph) installed in the 2meter telescope in Nishi Harima Astronomical Observatory in Japan.The grating was 150 l/mm, giving a spectral resolution of ∼ 600, with the 1. ′′ 2 slit width.We executed five 1200-second 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).

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-fiber 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 pipeline2 .

X-ray
The target was observed twice using the ROSAT (Truemper 1982) Position Sensitive Proportional Counter (PSPC) detector in 1991 December (ObsID = RP600164N00) and 1992 December (RP701178N00) with net exposures of 16.2 ks and 3.9 ks, respectively.The ROSAT data were processed with Standard Analysis Software System (SASS).The source spectra were extracted from a circular region with a radius of 150 ′′ centered 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 ′′ 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 ′′ radius, while the background was from an annulus between 90 ′′ and 150 ′′ radii.

DATA ANALYSIS AND RESULTS
In this section, we explain the analysis of optical spectra to investigate the variation of H.Then, we measured the time-lag of H, 1, and 2 to the continuum.

Variation of H𝛽
To obtain accurate light curves, we calibrated the spectral flux assuming [OIII]5007 luminosity did not change during all observations.The region where [OIII]5007 is emitted extends over ∼ 100 pc, so it's often assumed that the intensity remains unchanged over a few decades of observations.We fitted the spectra to determine their [OIII]5007 flux with PyQSOfit (Guo et al. 2018).The components applied in the fitting are a power-law continuum, three Gaussian components for the broad H, and one Gaussian each for the narrow H, [OIII]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.We use these continuum-subtracted spectra to fit the emission lines.The velocity offsets and the line width for [OIII]4959/5007, and the narrow H are set to be equal, and the flux ratio was fixed to 1:3:0.07.Here, the ratio of the narrow component of H 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 [OIII]5007 calculated with the Gaussian fitting.The total H 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 value of the light curve of H is summarized in the Appendix (Table A1 and Figure 3).
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 day to 3633 days), the calculation is weighted by the length of the observation period as a percentage of the light curve.The definition formula is as follows. (2) (3) Here,  represents each spectrum,  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 (≡  RMS ()) was assumed to be expressed in the following equation using the actual variable component (≡  variable ()) and observational error (≡ ).
To estimate the constant , we fitted the RMS spectrum with a Gaussian and  with Equation 4 (see the middle panel of Figure 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.
These line spectra (mean, RMS, and variable) are all derived from line spectra with the local continuum already subtracted.The FWHM and  values have been calculated within the wavelength range of 4600 Å to 5120 Å.For the mean spectrum, this calculation was performed after subtracting the narrow line components.The results of measuring the  and FWHM of H for the mean spectrum, the RMS spectrum, and the variable component are summarized in Table 1.The mean, RMS, and the variable component spectra for J1258 are shown in Figure 4.The mean spectrum represents the typical shape of the H and [OIII]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 [OIII]4959/5007 profiles are little variable, which indicates that the flux normalization using [OIII]5007 is working well.

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 H, WISE 1, and WISE 2, respectively.We used the results using two different software algorithms to obtain time-lags: interpolated cross-correlation function (ICCF; Sun et al. 2018;Peterson 1998), and JAVELIN (Zu et al. 2010(Zu et al. , 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 crosscorrelation 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 et al. 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 Figure 5.

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.

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  , , and Δ, which require certain assumptions for each.In this section, we discuss the values of these parameters.
Although the factor  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  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  = 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 Δ 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  value from Grier et al. (2013b), and adopts the ICCF estimate for (439 +68 −119 days), and (3429 ± 157 km/s) for the variable component of Δ.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 M ⊙ for the black hole mass calculated using Equation 1.

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 disk for the following reasons.First, the variability timescale of BELs and continuum is shorter than the orbital timescale 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(Hutsemékers et al. , 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 disk, based on its timescale.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.

𝛼 ox and Eddington ratio
We analyzed 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 statistics (Cash 1979) 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 Figure 6.The best-fitting parameters of observed-frame 2 keV ( 2 keV ) and 2-10 keV ( 2−10 keV ) fluxes, and rest-frame 2 keV ( 2 keV ) and 2-10 keV ( 2−10 keV ) luminosities3 , 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 AGNs (e.g., Ueda et al. 2014;Ricci et al. 2017).We find that the rest-frame 2 keV luminosities increase from ∼ 1 × 10 44 erg/s in 1991 and 1992 to ∼ 4 × 10 44 erg/s in 2021.
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 disk.We use the following definition to calculate  ox (Tananbaum et al. 1979), 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, FeII 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).The obtained  ox , bolometric luminosity, and Eddington's ratio before and after the state transitions are summarized in Table 4.We used optical spectra and X-ray data taken during Off and On states each.Since the optical light curve in Figure 7 shows that Table 4. Table summarizing the parameters representing the accretion disk state before and after state transitions.Each parameter is defined as follows: log( 2500 Å ) represents the  at 2500 Å,  bol is the estimated bolometric luminosity using Equation 6, the Eddington ratio is calculated from  bol and the black hole mass, and  ox is calculated using Equation 5.

Parameters
Off  The optical light curve with periods of optical spectra and Xray 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 = 54000 are considered as Off state, and the data after MJD = 57000 as On state.
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.
Jin et al. ( 2021) compared the  ox of 10 previously discovered CSQs with the simulated results from Sobolewska et al. (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.

Accretion rate variation
The measured Eddington ratio and  ox of J1258 in the Off and On state 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 Figure 8 together with the best-fit results of Jin et al. (2021) and Lusso et al. (2010).
One of the possible causes of the state transition is thermal disk instability (Noda & Done 2018).According to the disk instability model, the accretion disk forms a geometrically thin disk 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 disk is evaporating (Esin et al. 1997), resulting in a hot and radiatively inefficient accretion flow (Shapiro et al. 1976) 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 et al. 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  −5/3 (van Velzen et al. 2020).In the case of J1258, however, it shows brightening again from around MJD = 58000 (Figure 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 disk, as we mentioned in Nagoshi et al. (2021).

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 Figure 4 show the mean, the RMS, and the variable component of H, respectively.
Examining Figure 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- ure 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 Figure 9 represents the mean spectrum.The second spectrum from the top is the mean spectrum that the narrow lines (the narrow components in H and [OIII]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 Figure 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).
Indeed, since this object was observed during a relatively faint period, the double-peaked shape can be seen in real data.Figure 10 compares the H shape by subtracting the continuum and narrow emission lines (H and [OIII]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.
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 Figure 9), respectively.Here we name these components Region X and Region Y, respectively: Region X: Emitting a relatively slow (∼3,000 km/s), singlepeaked line that responds to the state transition of the accretion disk.
Region Y: Emitting a relatively fast (∼7,000 km/s), doublepeaked line, which is not likely affected by the state transition of the accretion disk.
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 center (Section 3).
The kinematics responsible for the double-peaked emission    4).This observation suggests that these emissions might occur independently of the accretion disk'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 disk'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 H 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 pc ∼ 1400 R g (≡ GM c 2 ), which can be explained by the rotation of the disk, given that the black hole mass is 10 9.64 +0.11 −0.20 M ⊙ and the typical speed is about 8,000 km/s (the peak velocity of the doublepeaked 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 years.Therefore, detecting a velocity shift over a 20-year observation period would be difficult.The distance to the dust torus is estimated by the response time of 1 to be 0.68 pc (805 light days).These locations are shown in Figure 11.
Assuming the structure of Figure 11, it is also possible to explain why Region Y was not affected by the state transition of the accretion disk.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 disk being larger than that of hard X-rays from the hot corona ( ∼ 10 9 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 disk 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 Figure 11 with model calculation.We adopt a model of Chen et al. (1989) here.This model can calculate the BEL spectrum, including relativistic effects, for a BLR distributed as a thin Keplerian disk.This model is defined by five parameters; the viewing angle , the inner edge of the BLR  1 , the outer edge of the BLR  2 , the emissivity index , and the local broadening parameter   .The model assumes that the surface intensity distribution of the disk is determined by Equation 7, which is a function of  1 and  2 and , and that locally emission lines are emitted with a velocity component of Keplerian rotation and a velocity dispersion of   .The surface intensity  () is defined as follows: where,  0 is a peak value of the emissivity and  is radius from the center of the disk.We show here as an example of parameters that can reproduce the qualitative features of Region Y, the results calculated with  = 47 • , 1 = 150  (  is Schwarzschild radius; 150  = 0.04 pc),  2 = 900  (= 0.25 pc),  = 0.9, and   = 1700 km/s are shown in Figure 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 disk-shaped BLR (Region Y) with a thin hole inside the BLR radius (at Region X) estimated by the time delay, as shown in Figure 11.Additionally, utilizing the parameters presented in Figure 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.Figure 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.

Origin of the two components of the BLR
Based on the structure shown in Figure 11, we discuss the origins of the gas in the BLR from two major theories: one is that the gas is supplied from disk winds from accretion disks (e.g., Gaskell 1982;Nicastro 2000;Nicastro et al. 2003;Czerny & Hryniewicz 2011) and the other is that it is derived from accreted components from the dust torus (e.g., Gaskell 2009Gaskell , 2010;;Hopkins et   Goosmann 2013; Wang et al. 2017).The disk-wind origin theory is based on a model in which outflow generated from the accretion disk according to an increase of the mass accretion rate winds up the gas, which is ionized by the emission from the inner accretion disk.This model suggests that no BLR is produced below a certain threshold because the disk 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 disk.
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 disk winds from the accretion disk 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 disk 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 disk 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 H.As a result, the measured time-lag is 4.1 +3.5 −3.1 days (Figure 14).If the gas in Region X was come from disk winds, the outflow reach Region X in 0.37 pc from the center in 20 years, which means the velocity is about 1% 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 timelags 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 doesn't have significant radial motion.
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 timelag of 1 is about twice the radius of Region X.The time-lag of 1 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 1 and 2 bands (see Figure A3), the predicted increase/decrease pattern and the actual observation are different, especially in the 2 band (the JAVELIN light curve shows a decline around MJD=58500 days, while the actual observation in the 2 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 years.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 timescale (< gas circling timescale: ∼ 200 years).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 disk 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).

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 M ⊙ .• The brightening event can be explained by the state transition of the accretion disk.
• The location of the BLR near the central core is inferred from the RM results (Figure 11).There are two components in the BLR, each located at 0.29 pc and 0.37 pc from the center.
• 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. .Spectra around the H, 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 1 .
Figure 1.Example of power-law continuum subtraction.The original spectrum is shown in gray, and the fitted continuum component is represented in red.The black dashed lines indicate the wavelength windows used for the continuum fitting, while the black solid 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 H.

Figure 2 .Figure 3 .
Figure2.Example of emission line decomposition.The black line presents the emission line spectrum with the continuum subtracted.The red line represents the spectrum after subtracting the narrow emission line components, isolating the broad H.The spectrum of the narrow line components is illustrated with a green dashed line.We have assumed that the Narrow Line components remain constant throughout the observation period.

Figure 4 .
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 H.The narrow [OIII]4959/5007 lines and asymmetric double peaked component of broad H appear only in the mean spectrum, except for weak residuals.

Figure 5 .
Figure 5.The Cross-Correlation Function and the centroid distribution to the continuum with H (top), WISE 1 (middle), and WISE 2 (bottom) of J1258.

Figure 6 .
Figure6.The unfolded X-ray spectra and best-fit 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.

Figure 7 .
Figure7.The optical light curve with periods of optical spectra and Xray 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 = 54000 are considered as Off state, and the data after MJD = 57000 as On state.

Figure 9 .
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.

Figure 10 .
Figure 10.Narrow lines subtracted H from SDSS.The black spectrum was taken on 2006/04/23, and the red spectrum was taken on 2016/04/03.

Figure 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 H 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 1 time-lag.

Figure 12 .
Figure 12.Example of reproduction of Region Y component spectrum by a model ofChen et al. (1989).The solid black line is the mean spectrum, which is decomposed into a narrow emission line (blue dash-dot line), a constant multiple of the variable component (blue dot line), and Region Y component (solid red line).The green dotted line shows an example of fitting the Region Y component with the model ofChen et al. (1989).

FluxFigure 13 .Figure 14 .
Figure13.Fitting results of the H emission line using two components: the variable component as Region X and the model ofChen et al. (1989)) as Region Y.Each panel, from top to bottom, displays the H and [OIII] 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.

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

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

Table 3 .
Best-fit X-ray spectral parameters, observed-frame absorptioncorrected fluxes, rest-frame luminosities, and reduced C statistics (Lusso et al. 2010on Comparison of the measured  ox and  Edd for J1258 in its On and Off states.The red dashed line and surrounding red 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 blue solid line and surrounding the blue region represent the  ox - Edd relationship and its associated errors for the On state, based on findings from the XMM-COSMOS survey(Lusso et al. 2010).
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.