Reverberation mapping of narrow-line Seyfert 1 galaxy I Zwicky 1: black hole mass

We report results of the first reverberation mapping campaign of I Zwicky 1 during $2014$-$2016$, which showed unambiguous reverberations of the broad H$\beta$ line emission to the varying optical continuum. From analysis using several methods, we obtain a reverberation lag of $\tau_{\rm H\beta}=37.2^{+4.5}_{-4.9}\,$ days. Taking a virial factor of $f_{_{\rm BLR}}=1$, we find a black hole mass of $M_{\bullet}=9.30_{-1.38}^{+1.26}\times 10^6 M_{\odot}$ from the mean spectra. The accretion rate is estimated to be $203.9_{-65.8}^{+61.0}\,L_{\rm Edd}c^{-2}$, suggesting a super-Eddington accretor, where $L_{\rm Edd}$ is the Eddington luminosity and $c$ is the speed of light. By decomposing {\it Hubble Space Telescope} images, we find that the stellar mass of the bulge of its host galaxy is $\log (M_{\rm bulge}/M_{\odot}) = \rm 10.92\pm 0.07$. This leads to a black hole to bulge mass ratio of $\sim 10^{-4}$, which is significantly smaller than that of classical bulges and elliptical galaxies. After subtracting the host contamination from the observed luminosity, we find that I Zw 1 follows the empirical $R_{\rm BLR}\propto L_{5100}^{1/2}$ relation.


INTRODUCTION
Narrow-line Seyfert 1 galaxies (NLS1s) are thought to be a special subclass of active galactic nuclei (AGN). Compared to the broad-line Seyfert 1 galaxies, NLS1s have: (1) narrower Balmer lines (FWHM Hβ 2000 km s −1 , by definition, where FWHM Hβ is the full width at half maximum of the broad Hβ emission line), (2) smaller intensity ratio of [O III] λ5007, to Hβ line ([O III]/Hβ < 3), (3) stronger optical Fe II multiplet emissions, and (4) usually steeper soft X-ray spectra and more rapid X-ray variability (Osterbrock & Pogge 1985;Goodrich 1989;Boroson & Green 1992;Boller et al. 1996;Sulentic et al. 2000;Wang et al. 2004). These distinctive properties of NLS1s can be explained by less massive black holes accreting with higher mass accretion rates (Boroson & Green 1992;Boller et al. 1996;Wang & Netzer 2003;Mathur & Grupe 2005;Grupe 2004;Peterson et al. 2004). Black holes in NLS1s are therefore undergoing fast growth through super-Eddington accretion (e.g., Kawaguchi et al. 2004;Wang & Zhang 2007). In the high-z universe, such fast growth was suggested to be a possible way of forming supermassive black holes (Volonteri & Rees 2005;Wang et al. 2006;Mortlock et al. 2011;Bamnados et al. 2018). Moreover, these super-Eddington accreting massive black holes (SEAMBHs) have saturated luminosity predicted by slim accretion disk model (Abramowicz et al. 1988;Wang & Zhou 1999) and could be used as a new kind of standard candle to study the expansion history of the high-z Universe Marziani & Sulentic 2014;Cai et al. 2018;Marziani et al. 2019) since they are quit common from low-z to high-z Universe (Du et al. 2016b;Negrete et al. 2018;Martínez-Aldama et al 2018). Reliably measuring black hole mass of local super-Eddington AGNs greatly helps elucidate these issues.
The reverberation mapping (RM) technique is a powerful tool for measuring the mass of black hole in the centre of AGN (Peterson 1993;Peterson et al. 2004;Peterson 2014). The widely accepted scenario is that gas around the central back hole is photoionized by the continuum emissions from the accretion disk and emits the observed broad emission lines. This is known as broad-line regions (BLR). Because of the light-travel time from the central black hole to the BLR, the variation in the flux of emission lines will delay that of the continuum. RM-campaigns monitor the flux variations in the continuum and the broad emission lines to measure the lags between them. Thus the observed lags are regarded as a measurement of the radius of the BLR. Assuming that the motion of BLR gas is governed by the gravity of the black hole, the profiles of these lines would provide kinematic information of the BLR and the viral mass of the black hole can be estimated by where f BLR is the virial factor, G is the gravitational constant, R BLR = c × τ Hβ is the emissivity-weighted size of the BLR emitting the broad Hβ line, with c the light speed and V FWHM the velocity of the BLR clouds inferred from the width of the broad Hβ line. Recently, the mass of the black hole in 3C273 has been determined by another novel technique using GRAVITY on VLTI (Very Large Telescope Interferometer) (Sturm et al. 2018), and the result is in good agreement with that given by a 10-year RM-campaign .
As one of the most famous NLS1s, I Zwicky 1 (PG 0050+124 or Mrk 1502, z = 0.061; hereafter I Zw 1) is selected as one of candidates in our large RM-campaign focusing on AGNs with SEAMBHs . I Zw 1 is a nearby and bright (m V = 14.06 ± 0.05 from Slavcheva-Mihova & Mihov 2011b) prototypical NLS1 (Schmidt & Green 1983;Osterbrock & Pogge 1985;Halpern & Oke 1987) with FWHM Hβ 1200 km s −1 , large relative strength of optical Fe II emission lines (R Fe = 1.47 from Boroson & Green 1992, where R Fe = F Fe /F Hβ ), weak intensity of [O III] line and a steep X-ray (2 − 10 keV) photon index (Γ = 2.15 ± 0.02 from Gallo et al. 2007). Its Fe II spectrum is commonly used as a template to fit the optical and ultraviolet spectra of AGNs and quasars (Boroson & Green 1992;Vestergaard & Wilkes 2001;Véron-Cetty et al. 2004). Similar to other NLS1s, Hubble Space Telescope (HST) images show clearly that the host galaxy of I Zw 1 is a spiral galaxy with asymmetric knotty spiral arms (Slavcheva-Mihova & Mihov 2011a) consisting of young stellar populations and ongoing nuclear star formation activity (Hutchings & Crampton 1990;Scharwächter et al. 2007). It would be highly instructive to measure the black hole mass and accretion rate of I Zw 1 in order to understand the relation of its central engine to its host galaxy. As far as we know, only Giannuzzo et al. (1998) obtained seven epochs with a mean cadence of about 170 days to monitor optical variability of I Zw 1. This is not enough for RM measurements.
In this paper, we report results of the first RM campaign for I Zw 1. In §2, we describe the observations and data reduction. In §3, we analyze the light curves to measure the Hβ time lags for black hole mass and fit the HST images to separate the host component to test the empirical R BLR -L 5100 relation. Brief discussions are provided in §4. A summary is given in the last section. Throughout this work, we adopt a cosmology with H 0 = 67 km s −1 Mpc −1 , Ω Λ = 0.68 and Ω m = 0.32 (Ade et al. 2014) .

Observations
Our observations were performed from November 2014 to February 2016 using the Lijiang 2.4 m telescope of Yunnan Observatories, Chinese Academy of Sciences. A versatile instrument named Yunnan Faint Object Spectrograph and Camera (YFOSC) was adopted since it can switch from spectroscopy to photometry within one second . During the spectroscopic observations, we used Grism 14, which covers the wavelength range of 3800−7200Å, and a long slit with a fixed width of 2. 5 to minimize flux loss caused by the different seeing conditions. The final spectral resolution was roughly 500 km s −1 . The spectral resolution was estimated by comparing the width of the [O III] line of the Lijiang spectrum with the one obtained from the Sloan Digital Sky Survey ?). For accurate flux calibration, we adopted the method described by Maoz et al. (1990) and Kaspi et al. (2000), putting I Zw 1 and a nearby comparison star simultaneously in the long slit during the exposure. The comparison star (R.A. J2000 = 00:53:46.92, Dec. J2000 = + 12:39:56.8) is a stable G star without intrinsic changes, as confirmed from the photometric light curve shown in Appendix A. The position angle between the comparison star and I Zw 1 is 120 • , and the angular distance is 3.4 . Flux changes caused by unstable weather can be corrected by the comparison star. To mitigate the influence of cosmic rays, we used two consecutive exposures with typical exposure time of 600 s each. Between the two exposures, we took an image to inspect whether I Zw 1 and the comparison star still sit in the middle of the slit. If not, we adjust the slit for realignment. Every night we also took two consecutive exposures of a spectrophotometric standard star. In order to minimise the impact of atmospheric differential refraction, we took spectra only when the air mass was 1.5 (median air mass 1.3 for our data). In addition, we also made V -band photometric observations to verify the accuracy of the spectral calibration. Three consecutive 90 s exposures were taken to reduce the impact of cosmic rays. We obtained 66 epochs of spectroscopic observations and 74 epochs of photometric observations for I Zw 1 during our campaign. The median cadence is 4 days.

Data reduction
All the spectroscopic and photometric data were reduced following standard methods using IRAF v2.16. For spectroscopic data, the extraction aperture was 8. 5. Standard neon and helium lamps were used for wavelength calibration, and the comparison star was used for flux calibration, as described by Du et al. (2014). Here we briefly describe the method of flux calibration. (1) We used spectrophotometric standard stars to calibrate the spectra of the comparison star. (2) A fiducial spectrum of the comparison star was made from its flux-calibrated spectra taken over several nights under good weather conditions. The uncertainty of the absolute flux calibration was ∼ 10% (Du et al. 2016a). (3) By comparing the observed spectrum of the comparison star to the calibrated fiducial spectrum, we obtained a sensitivity function for each exposure. (4) This sensitivity function was then used to calibrate the spectrum of I Zw 1. (5) We combined the separate exposures taken each night into a single-epoch spectrum for that night. This method can achieve an accuracy of ∼ 3% in the calibrated spectra , which is especially necessary for targets such as I Zw 1 and other SEAMBHs, whose [O III]λ5007 line tends to be too weak to be used for relative flux calibration, as conventionally practiced for RM (Foltz et al. 1981;Peterson et al. 1982).
For generating the photometric light curves of I Zw 1, six stars in the same fields were used to perform differential photometry. The radius of the aperture for I Zw 1 was 5. 1, and that for the background was 8. 5 − 17. 0. Specifically, for each exposure, we obtained the instrumental magnitudes of I Zw 1 and the six stars by IRAF. We calculated the differential magnitude between I Zw 1 and the mean of the six stars, and the uncertainties in all the instrumental magnitudes were propagated to the uncertainty in this differential magnitude. Then we averaged the differential magnitudes of the three exposures in the same night as the value of that individual night, and calculated a statistical uncertainty from the uncertainty of each exposure by error propagation. In addition, we calculated scatters between the differential magnitudes of the three exposures as the systematic uncertainty. The two uncertainties are added in quadrature as the final uncertainty of the differential magnitude in each night.

Light curves
We use the following method to determine the 5100Å flux density and the Hβ flux. The 5100Å flux is the median flux between 5075 and 5125Å in the rest frame. The Hβ flux is obtained by integrating the flux between 4810−4910Å from the continuum-subtracted spectrum, while the continuum underlying the Hβ line is determined by a linear interpolation between two continuum bands (4740-4790Å and 5075-5125Å). The continuum windows are selected to minimize the contamination from other emission lines, such as [O III], Fe II, and He II. In Figure 1, we mark the window for the Hβ line by the red band and the windows for the continuum by the grey bands. The measurement errors of the light curves come from both Poisson noise and systematic uncertainty. The systematic uncertainty arises from poor weather conditions, slit positioning, telescope tracking, and choice of continuum window. It can be estimated from the standard deviation of the residuals after subtracting a median-filtered, smoothed light curve. Table 1  To improve the sampling of the continuum light curves, we used the V -band photometry data from ASAS-SN (All-Sky Automated Survey for SuperNovae) 1 . ASAS-SN is a long-term project aiming to monitor the entire visible sky to a depth of V < 17 mag with a cadence of 2-3 days (Shappee et al. 2014;Kochanek et al. 2017 We merge the V -band photometry data from the Lijiang and ASAS-SN into the 5100Å light curves by applying a multiplicative scale factor and an additive flux adjustment, which are determined by a Markov-chain Monte Carlo (MCMC) implementation (Li et al. 2014;Du et al. 2015). The combined continuum light curves are shown in Figure 7i, which contain 136 epochs from the 2014-2015 season and 149 epochs from the 2015-2016 season. By combining ASAS-SN data, the continuum light curve shows an obvious structure for 2014-2015 season (shown in Figure 7i), which corresponds to that in the light curve of Hβ and the seasonal gap is shortened to 94 days. However, as shown in Figure 7i, the 5100Å light curves increase the scatter of the continuum light curve, so we merge only the V -band photometry data from ASAS-SN into the photometry data from Lijiang to obtain the final photometry continuum light curve shown in Figure 2i. The final photometric continuum light curve contains 109 and 110 epochs from the 2014-2015 and 2015-2016 seasons, respectively. We adopt the photometric continuum light curve to obtain the Hβ lags in this paper. However, for comparison, we also show the results for the 5100Å and the combined continuum light curves in Appendix B.

Variability characteristics
We use the methods described by Rodríguez-Pascual et al. (1997) to calculate light curve characteristics of the continuum and the Hβ emission line. The variability characteristics are given by where and the error of F var is given by (Edelson et al. 2002), where N is the total number of data, F i is the flux of the i-th observation, and ∆ i is the uncertainty of F i . We use F 5100 var , F combine var , F V var and F Hβ var to denote the amplitude of the variation in the light curves of the 5100Å continuum, the combined continuum, the photometry continuum and broad Hβ line. Table 2 lists F 5100 var , F combine var , F V var and F Hβ var for the 2014-2015 data, 2015-2016 data, and 2014-2016 data, respectively. We find that F 5100 var , F combine var , F V var and F Hβ var of I Zw 1 are clearly much smaller than that of other PG quasars (see Table 5 in Peterson et al. 1998 and Table 5 in Kaspi et al. 2000), by a factor of a few.
The mean and the rms spectra of I Zw 1 are calculated bȳ We show the results for the 2014-2015 season in Figure 1. The [O III] line disappears in the rms spectrum, which implies that the spectral calibration is quite good. We also plot the mean of the error of the individual spectra (see the red line in Figure 1), which can be taken as the rms spectrum for the variance caused by noise only (Hu et al. 2016).

Hβ lags
We employ three methods to calculate Hβ lags: the interpolation cross-correlation function (ICCF; Gaskell & Sparke 1986;Peterson et al. 2004), JAVELIN (Zu et al. 2011), and MICA 2 (Li et al. 2016). The ICCF method calculates the cross-correlation function (CCF) between the light curves of the continuum and Hβ fluxes. The time lag is determined by measuring either the location τ peak of the CCF peak (r max ) or the centroid τ cent of the points around the peak above a threshold r ≥ 0.8 r max . The associated uncertainties are estimated by the Monte Carlo flux randomization/random subset sampling (FR/RSS) method of Peterson et al. (1998) and are given at 63.8% confidence levels. Both JAVELIN and MICA are forward-modelling methods that use the damped random walk model (e.g., Kelly et al. 2009;Zu et al. 2011;Li et al. 2013) Table 3 summarises the time lags determined by the above three methods. They are consistent to each other, within uncertainties.
We detect significant time lags with small uncertainties using the photometric light curves. The Hβ time lags for the 2015-2016 data have large uncertainties (Figures 2g and 2h) because of the small variability amplitude and large flux dispersion of the light curves (see Figures 2e and 2f). In principle, JAVELIN and MICA can reconstruct the light curves during the large seasonal gaps, such that, for data for the whole campaign, the lags obtained by JAVELIN and MICA should be more robust than that derived by ICCF. We also find that the lags obtained using JAVELIN and MICA (Figure 2l) are more consistent with that obtained by ICCF for the 2014-2015 data (Figure 2c). From inspection of the CCF peak value r max , we adopt τ cent = 37.2 +4.5 −4.9 days from the ICCF analysis of the photometric continuum and Hβ light curves from 2014−2015, and we use this value to calculate the black hole mass below.

Host galaxy
We use images obtained with the HST to study the host galaxy, in particular its overall morphology, bulge-to-total ratio, and bulge stellar mass. We also need to determine the degree to which the spectroscopically measured optical continuum is contaminated by host galaxy emission. The highest quality HST observations currently available are those obtained with the WFC3 camera on 2 November 2013 (GO-12903, PI: Luis C. Ho). I Zw 1 was observed for 300 s with the F438W filter (395-468 nm) in the UVIS channel and for 147 s with the F105W filter (901-1204 nm) in the IR channel. An additional short (40 s) F438W exposure was taken to warrant against saturation of the nucleus in the long exposure. These two filters were purposefully chosen to mimic B and I in the rest frame of I Zw 1, at the same time avoiding contamination from strong emission lines. The large wavelength separation between the two filters also offers the most leverage for constraining the stellar population (Section 5.2). To better sample the point-spread function (PSF), the long UVIS observation was taken with a three-point linear dither pattern, while the IR observation was taken with the four-point boxy dither pattern. To avoid overheads due to buffer dump, we employed the UVIS2-M1K1C-SUB 1k×1k subarray for the UVIS channel and the IRSUB512 subarray for the IR channel, which resulted in a restricted field-of-view of 40 ×40 and 67 ×67 , respectively.
Because of the sparseness of field stars in the vicinity of I Zw 1, the observations were conducted in GYRO mode, which, unfortunately, led to considerable degradation of the PSF of the dither-combined image generated from the standard data reduction pipeline. Instead, we use the DrizzlePac task AstroDrizzle (v1.1.16) to correct the geometric distortion, align the sub-exposures, perform sky subtraction, remove cosmic rays, and, finally combine the different exposures. The core of the AGN was severely saturated in the long F438W exposure, and it was replaced with an appropriately scaled version of the 40 s short exposure. Because the PSF of HST is undersampled, we broaden both our science and PSF images by convolving them with a Gaussian kernel so that the PSF can reach Nyquist sampling . This largely removes the sub-pixel mismatch of the core of the PSF and the smearing of the nucleus due to image drift.
The F105W image ( Figure 3a) reveals a nearly face-on spiral galaxy with two dominant spiral arms. A bar-like structure may also be present. The F438W image (Figure 3b) is considerably shallower, but the host is still clearly detected. To extract quantitative measurements of the bulge, we use the program GALFIT (v3.0.5; Peng et al. 2002Peng et al. , 2010 to fit two-dimensional surface brightness distributions to the HST images. A crucial ingredient is the PSF, which will have a strong effect on the bright and active nucleus. Unfortunately, no suitably bright star is available to be used as the PSF within the limited field-of-view of the subarray WFC3 images. Instead, we generated a high-S/N stacked empirical PSF by combining a large number (24 for F105W and 57 for F438W) of bright, isolated, unsaturated stars observed during the course of other WFC3 programs. Extensive tests, consisting of fits to isolated bright stars, indicate that our stacked empirical PSF is far superior to synthetic PSFs generated from the TinyTim program (Krist & Hook 1999), and it has higher S/N than the PSFs of individual stars. The reduced χ 2 of the fits are ∼3 times larger for the TinyTim synthetic PSF. Comparison of empirical PSFs observed from different programs indicate that the WFC3 PSF does not vary significantly with time (< 10%).
We concentrate first on obtaining the best global fit on the deeper F105W image, whose red wavelength is also more sensitive to the host. After much experimentation, we adopt a model with three components: a point source (represented by the PSF) for the nucleus, a bulge parameterized as a Sérsic (1968) function with index n, and an exponential disk (equivalent to a Sérsic function with n = 1), with coordinate rotation turned on to fit the spiral arms. We could not obtain a robust solution that considers the faint, bar-like feature, and in the end we did not treat it. We use the m = 1 Fourier mode of the disk, which is sensitive to lopsidedness, to gauge the degree of global asymmetry of the galaxy (see, e.g., Kim et al. , 2017). The best model ( Figure  3a) reveals a bulge with n = 1.73 ± 0.06, formally but barely below the conventional threshold for pseudo-bulges (n < 2; Fisher & Drory 2008) and an overall B/T = 0.52 ± 0.04. The disk is moderately asymmetric, with Fourier amplitude a 1 = 0.11 ± 0.01. As the host is considerably weaker in the F438W image, we fit it keeping the structural parameters fixed to the values obtained from the F105W model, solving only for the magnitude.
The uncertainties of the decomposition are dominated by PSF mismatch in the nucleus. We estimate the effect of the PSF by generating variants of the empirical PSF by combining different subsets of stars, and then repeating the fit. The final error budget is the quadrature sum of these two contributions.

The R − L relation
A proper evaluation of the R − L relation must consider the influence of host galaxy contamination on the luminosity (Kaspi et al. 2000;Bentz et al. 2013). We use the decomposition of the HST images of I Zw 1 to estimate the contribution of the host galaxy within the spectral extraction aperture. The host flux at 5100Å is transformed from the F438W magnitude with the IRAF task synphot, assuming a stellar population template with an age of 5.0 Gyr (Section 4.2), after correcting for redshift (z = 0.0611) and Galactic reddening [E(B − V ) = 0.057 mag; Schlafly & Finkbeiner 2011] using the extinction curve of Cardelli et al. (1989). Table 4 lists the observed total, host galaxy, and AGN fluxes at 5100Å. With an AGN luminosity of L 5100 = 3.19 × 10 44 erg s −1 , I Zw 1 follows the empirical R BLR − L 5100 relation (Figure 4a).  (Table 4). Adopting τ Hβ = 37.2 days and f BLR = 1, Equation (1) yields M • = 9.30 +1.26 −1.38 × 10 6 M . It is known that the observed kinematics of the BLR are generally influenced by inclination effects, which will lead to uncertainties of the virial factor (Krolik 2001;Collin et al. 2006).
In principle, f BLR can be calibrated using the M • − σ relation of inactive galaxies, and the current best estimates of f BLR are 1.3 ± 0.4 for classical bulges and 0.5 ± 0.2 for pseudo bulges, when f BLR is calibrated using FWHM based on mean spectra (Ho & Kim 2014). The host galaxy of I Zw 1 likely contains a pseudo bulge (Section 5.2). However, the large scatter of the M • − σ relation of pseudo bulges (Kormendy & Ho 2013) introduces significant uncertainty into the estimate of f BLR . For consistency with our previous work on SEAMBHs ?), we adopt f BLR = 1.0 (Netzer & Marziani 2010;Woo et al. 2013). We also estimate the black hole mass from the FWHM of the rms spectrum (V rms FWHM = FWHM rms Hβ = 606 +28 −28 km s −1 ) using f rms BLR = 1.12 (Woo et al. 2015) and obtain M rms • = 2.99 +0.46 −0.48 × 10 6 M . Meanwhile, polarized spectra of type 1 AGNs can provide invaluable information for estimating black hole masses. Considering electron scattering in the equatorial plane (Smith et al. 2005), a polarized spectrum is equivalent to a spectrum seen by an observer on the mid-plane. Thus, polarized spectra can yield more reliable estimates of the black hole mass (Afanasiev & Popovic 2015;Baldi et al. 2016;Songsheng & Wang 2018). For equatorial scattering of a BLR with a flattened geometry, where the virial factor f BLR ≈ 0.24 (with small uncertainty) and V FWHM is the FWHM of the broad-line profile in the polarized spectrum (Figure 3 in Songsheng & Wang 2018). Fortunately, a polarized spectrum of I Zw 1 for the Hα emission line has been taken using the 6 m telescope of the Special Astrophysical Observatory of the Russian Academy of Science (L. Popović 2018, private communications) 3 . We find FWHM Hα = 1983 km s −1 from the polarized spectrum. The relation between the widths of Hβ and Hα in the polarized spectrum is unclear, but a reasonable assumption is that the same relation is followed as in the total spectrum if 1) both lines are scattered by the same population of electrons, and 2) the regions emitting Hβ and Hα lines are much smaller than the electron scattering region. So we use the relation FWHM Hβ = 1.07 × 10 3 FWHM Hα /10 3 km s −1 1.03 from Greene & Ho (2005), and obtained V FWHM = FWHM Hβ = 2165 +189 −182 km s −1 . This value needs to be tested observationally in the future. Equation (6) then yields M • = 8.18 +1.74 −1.75 × 10 6 M , which agrees remarkably well with the mass estimates based on the mean spectrum.

Accretion rates
Photon trapping causes the radiated luminosity of slim disks to saturate (Abramowicz et al. 1988;Wang & Zhou 1999;Mineshige et al. 2000). Under these conditions, the normal Eddington ratio L bol /L Edd is very insensitive to the accretion rate and mainly depends on the black hole mass, where L Edd is the Eddington luminosity and L bol is the bolometric luminosity. This renders L bol /L Edd unsuitable to indicate the accretion rate of SEAMBHs. As derived from the self-similar solution of slim disk, the photon trapping radius is R trap ≈ 72 (Ṁ /80)R g , where R g = GM • /c 2 , and the radius emitting optical (5100Å ) photons is R 5100 ≈ 4.3 × 10 3 m −1/2 7 R g (Wang & Zhou 1999), which is much larger than R trap . So optical photons can escape freely, and hence accretion rates can be reliably estimated from the formalism of the standard accretion disk model (Shakura & Sunyeav 1973). The dimensionless accretion rate is defined asṀ =Ṁ • /L Edd c −2 , withṀ • is the mass accretion rate. Given the 5100 A luminosity and the black hole mass,Ṁ = 20.1( 44 / cos i) 3/2 m −2 7 , where 44 = L 5100 /10 44 erg s −1 , m 7 = M • /10 7 M , and i is the inclination angle between the line-of-sight and the accretion disk yieldingṀ = 203.9 +61.0 −65.8 , for cosi = 0.75. As discussed in ?, cosi = 0.75 is a reasonable mean value for type 1 AGNs. If we take ∆ log cos i 0.1, then the uncertainty onṀ due to i will be ∆ logṀ = 1.5 ∆ log cos i 0.15. The inclination also has insignificant effect on the observed width of the broad emission lines V FWHM . If we know the Keplerian velocity V K and the height H BLR of the flattened BLR, by the zeroth-order approximation, Li et al. (2013) and Pancoast et al. (2014) suggest H BLR /R ∼ 1, thus sini ≤ H BLR /R, which implies the inclination has insignificant influence on V FWHM , hence M • and thenṀ . The high dimensionless accretion rate indicates that I Zw 1 is a SEAMBH. Note that the accretion rate could be even higher by a factor of ∼ 10 than the present, if the FWHM obtained from rms spectrum are used in the calculation.

The R − L relation
There is growing evidence that the empirical R BLR -L 5100 relation does not apply to SEAMBHs. For a given luminosity, SEAMBHs tend to have a shorter Hβ lag, suggesting that the size of the BLR is related not only to luminosity but also to the accretion rate ?, 2018). According to the unified scaling relation for sub-and super-Eddington AGNs suggested by ?, I Zw 1 is predicted to have an Hβ delay of 12.5 days, which is much shorter than our measured value of 37.2 days. Unlike other SEAMBHs, I Zw 1 actually follows the standard R BLR -L 5100 relation and does not show an obvious shortened Hβ lag ( Figure  4a). Defining the deviation ∆R BLR ≡ log(R BLR /R R−L ), where R R−L is given by the empirical R BLR − L 5100 relation, Figure  4b plots ∆R BLR versusṀ . I Zw 1 deviates from the trend defined by most SEAMBHs. So does Mrk 493. The reasons for these outliers are unclear. One possibility is that some lags are too short to detect, given the current observation cadence. Within the SEAMBH framework described by Wang et al. (2014b) and Du et al. (2018), there are two BLRs: a normal one unshadowed, and another shadowed by the inner part of the disk and much closer to the central black hole. The emissivity-weighted gas distribution yields two peaks in the transfer function, corresponding to the unshadowed and shadowed BLRs respectively. Such a two-region BLR scenario gets supported from a recent modelling of the observed RM data in Mrk 142 by Li et al. (2018), which constructed flexible one-region and two-region models that included the spatial distribution and anisotropic emissivity of BLR clouds. They found that the two-region model is preferrable to the one-region model. However, the observed Hβ lag depends on the data quality, especially the cadence. For example, low cadence will smooth the light curves and yields a single broad peak in the CCF, although there are two peaks in the transfer function. To measure the potential shorter lag, campaign with higher cadence than the present are being planned for getting detailed BLRs of I Zw 1 and Mrk 493 as done in Mrk 142.

Black hole and bulge stellar masses
The GALFIT decomposition of the HST images of I Zw 1 yields a bulge magnitude 17.23±0.04 mag in F438W and 14.07±0.07 mag in F105W. We apply K-correction to convert the HST-based magnitudes to rest-frame magnitudes in B and I. We generate a series of template spectra with ages spanning 1 to 12 Gyr, adopting a Bruzual & Charlot (2003) models with solar metallicity, a Chabrier (2003) stellar initial mass function, and an exponentially decreasing star formation history with a star formation timescale of 0.6 Gyr. After accounting for Galactic extinction and redshift, we convolve the spectra with the response functions of the HST filters and generate synthetic F438W and F105W magnitudes. The observed bulge color of F438W−F105W = 3.16 ± 0.08 mag is best matched with a stellar population of age 5.0 +1.0 −0.5 Gyr, resulting in rest-frame M I = −22.70 ± 0.07 mag and B − I = 1.96 ± 0.09 mag. From Table 4 of Bell & de Jong (2001), we derived M bulge = 10 11.18±0.07 M , assuming solar metallicity and a Salpeter (1955) initial mass function. A Chabrier initial mass function gives 45% less mass (Longhetti & Saracco 2009), and therefore M bulge = 10 10.92±0.07 M .
For M • = 10 6.97 M , M • /M bulge ≈ 10 −4 , which is lower by a factor of ∼ 50 than the value inferred from the M • − M bulge relation for classical bulges and elliptical galaxies, but lies within the large scatter and lower ratios found for pseudo bulges (Kormendy & Ho 2013). I Zw 1 very likely hosts a pseudo bulge, in view of its relatively low Sérsic index of n = 1.73 ± 0.06 (Fisher & Drory 2008) and evidence for recent (see discussion above) and ongoing (Scharwächter et al. 2007) nuclear star formation. 3. We derive the stellar for the bulge of the host galaxy from detailed decomposition of HST images. We find a mass ratio of M • /M bulge ≈ 10 −4 , much lower than that for nearby inactive classical bulges and elliptical galaxies.
I Zw 1 is famous for its strong and narrow Fe II lines; however, the current data do not allow us to successfully measure a lag for Fe II. We are in the process of scheduling a new monitoring campaign with improved cadence and homogeneity in observation more favourable for detecting Fe II lags.   NOTE-F obs and F gal are the observed and host galaxy flux densities at (1 + z)5100Å in the observed frame. L 5100 is the mean AGN luminosity λL λ at rest-frame 5100Å after subtracting the host galaxy and correcting for Galactic extinction.     In order to check the stability of the comparison star during the campaign, we obtain its photometric light curve using five stars in the same field to perform differential photometry. The light curve of the comparison star is shown in Figure 5. The standard deviation is 1.9% , demonstrating that the comparison star does not vary significantly.
B. LIGHT CURVES AND Hβ LAGS Figure 6 and 7 show the results of the correlation analysis using the 5100Å and combined continuum light curves. For the 2014-2015 data, the Hβ time lag cannot be constrained well solely using the 5100Å light curves. As shown in Figures 6c and  6d, the CCF is very spiky and the posterior distributions of the time lags obtained by JAVELIN and MICA are not well-defined. This may be caused by the short duration and the large scatter of the continuum light curves. As shown in Tables 5 and 3, the uncertainty of the time lags derived using only the photometry continuum light curves is smaller than that using the combined continuum light curves.
In Figure 8, we show the ASAS-SN light curve of I Zw 1 from 2012 to the end of 2018. Generally, the source was less variable in those 7 years except for the period from 2014-07 to 2016-01. During the whole period, the variability of the continuum flux is F var (%) = 0.63 ± 0.02, which is much smaller than that of the 2014-2016 period. We were lucky in this SEAMBH campaign and captured the largest variations in this period. Figure 8 shows that it is less variable again after the present campaign, which implys that we may need a longer waiting for sharp variations. Higher cadence than 2-3 days is also necessary for more details of the BLR structure. On the other hand, the long-term stability of I Zw 1 supports SEAMBHs as cosmic candles for cosmology Cai et al. 2018). . Light curves and results of correlation analysis using the combined continuum light curves for I Zw 1 (same as Figure 2).