Unveiling the Small-scale Jets in the Rapidly Growing Supermassive Black Hole IZw1

Accretion of black holes at near-Eddington or super-Eddington rates represents the most powerful episode driving black hole growth, potentially occurring across various types of objects. However, the physics governing accretion and jet–disk coupling in such states remains unclear, primarily due to the difficulty in detecting associated jets, which may emit extremely weakly or exhibit episodic behavior. Only a few near/super-Eddington systems have demonstrated radio activity, and it remains uncertain whether jets exist and what their properties are in super-Eddington active galactic nuclei (AGNs) and ultraluminous X-ray sources. This uncertainty stems mainly from the complex radio emission mix, which includes contributions from jets, star formation activity, photoionized gas, accretion disk wind, and coronal activity. In this work, we conducted high-resolution, very long baseline interferometry observations to investigate jets in the highly accreting narrow-line Seyfert I system I Zw 1. Our observations successfully revealed small-scale jets (with a linear size of ∼45 pc) at both 1.5 and 5 GHz, based on the high radio brightness temperature, radio morphology, and spectral index distribution. Additionally, the parsec-scale jet observed in I Zw 1 displays a knotted morphology reminiscent of other sources accreting at similar rates. In summary, the high accretion rates and jet properties observed in the AGN I Zw 1 may support the AGN/X-ray binary analogy in this extreme state.

The Eddington ratio1 is a key indicator of black hole accretion and ejection states, in both stellar-mass black holes (SBHs) and supermassive black holes (SMBHs) (e.g.Fender et al. 2004;Falcke et al. 2004;Körding et al. 2006), which is generally < 1 under the assumption of spherical accretion.The accretion flows and associated ejection processes with low and moderate Eddington ratios generally can be described with Advection Dominated Accretion Flows (ADAFs, Narayan & Yi 1994;Esin et al. 1997) and standard accretion disk (or Shakura-Sunyaev disk, SSD, Shakura & Sunyaev 1973), respectively, and corresponding revisions (e.g.Fender et al. 2004;Done et al. 2007;Yuan & Narayan 2014).
However, super-critical accretion (with the super-Eddington accretion rate for which the black hole radiates above the Eddington luminosity) is viable in both observations and physics and potentially plays an essential role in feeding the black hole growth in the early Universe (see Yang et al. 2020, and references therein).Furthermore, super-Eddington accretion of the first-generation SMBHs may have a deep impact on regulating the (host) galaxy evolution and the epoch of reionization through feedback processes.As accretion increases to near or super-Eddington rates, the standard disc geometry cannot be maintained and the accretion flow will inevitably evolve into a 'slim disc' (Vierdayanti et al. 2013).The corresponding state is sometimes called the 'ultraluminous state' (Gladstone et al. 2009).
Regardless of the importance and the viability of super-Eddington accretion, our understanding of the accretion and ejection processes in this accretion state remains limited, which is primarily due to that only a few systems can temporarily trigger super-Eddington accretion (e.g.Greiner et al. 2001;Dai et al. 2018), and even fewer systems can maintain long-lived super-Eddington accretion (e.g.Middleton et al. 2021;King et al. 2023).Even worse, the mechanism for sustaining super-Eddington accretion in those sources which persistently accrete at super-Eddington rates is unclear.
It is also widely accepted that supermassive and stellar-mass black holes have similarities in accretion physics, i.e., Active galactic nuclei (AGNs) and XRBs have similar accretion state transitions and associated ejection processes.However, it is still unclear whether the AGN/XRB analogy holds in the 'ultraluminous state' and whether the geometry of the disc-corona system and jet-disc coupling are similar.Here, our interest is the connection between the short-lived canonical 'very high state' (universally found in XRBs) and the longstanding super-Eddington accretion in, for example, the microquasar SS 433, and ULXs, and which parameters are driving the long-lived super-Eddington accretion.As the time scale of state transition is proportional to black hole mass (Svoboda et al. 2017;Yang et al. 2020), a 'very high state' in SMBHs (e.g.M BH = 10 7 M ), would last 10 6 times longer than in 10 M stellar-mass black holes found in XRBs.Therefore, the study of near/super-Eddington AGNs provides an opportunity to understand the ejection process in a quasi-steady 'very high state' and may shed light on the physics to sustain such a near/super-Eddington accretion.For this reason, we present very long baseline interferometry (VLBI) observations of an AGN, I Zw 1, which radiates close to or above the Eddington limit.
I Zw 1 is one of the closest quasars located at a redshift of z = 0.0589 (Ho & Kim 2009) and is regarded as an archetypal narrow-line Seyfert 1 galaxy (NLS1) based on its optical properties (Schmidt & Green 1983;Pogge 2000).The black hole mass of I Zw 1 was estimated to be M BH = 9.3 × 10 6 M from optical reverberation mapping (Huang et al. 2019).The bolometric luminosity estimated from the spectral fitting is log L bol = 45.50 − 45.68 erg s −1 (Martínez-Paredes et al. 2017), which exceeds its Eddington luminosity with an Eddington ratio of λ Edd = 2.77−4.20.Another work obtained a higher black hole mass of M BH = 2.8 × 10 7 M using X-ray reverberation (Wilkins et al. 2021) and estimated the Eddington ratio of I Zw 1 as unity (or 0.3) based on the optical monochromatic luminosity (or the X-ray luminosity).However, the authors note that the luminosity could be underestimated due to photons being trapped in the disk.If we take the bolometric luminosity of log L bol = 45.50 − 45.68 erg s −1 (Martínez-Paredes et al. 2017), which is thought to be more accurate than the estimation from the single-band luminosity, and use the larger black hole mass measurement of M BH = 2.8 × 10 7 M (Wilkins et al. 2021), then the Eddington ratio would be 0.92 − 1.40.We should bear in mind that when the bolometric luminosity approaches and exceeds the Eddington luminosity, the actual mass accretion rate would be significantly higher than expected from the observable luminosity assuming a typical radiative efficiency (η < 1, Bian & Zhao 2003) because of the photon trapping effect (Mineshige et al. 2000).The radiative efficiency of I Zw 1 was estimated, as one of the Palomar-Green (PG) quasars PG 0050+124, to be log η = −2.21 or −1.18 ± 0.04 (based on the mass estimates with the broad emission line widths and the M − σ * correlation, respectively, see Davis & Laor 2011).With such a high Eddington ratio and low radiative efficiency, therefore, the SMBH in I Zw 1 must be growing with a mass accretion rate notably higher than the Eddington limit, i.e. the super-Eddington accretion rate.
I Zw 1 is also an extremely radio-quiet AGN with radio loudness parameter R = 0.35 (fn.2 ) (Yang et al. 2020).Radio emission from radio-quiet AGNs is complex and remains a subject of debate (Panessa et al. 2019).On the other hand, the presence of jets in near or super-Eddington systems and how it is launched are also questions that need to be explored.X-ray observations of I Zw 1 indicate that the X-ray corona exhibits some structure and part of it may be collimated and ejected (Gallo et al. 2007;Wilkins et al. 2017).I Zw 1 is an extreme example of a nearby highly accreting and radioquiet AGN, providing an ideal laboratory for studying outflow activities3 including jets and winds, with high spatial resolutions.
In this work, we report the Very Long Baseline Array (VLBA) and European VLBI Network (EVN) plus the enhanced Multi-Element Remote-Linked Interferometer Network (e-MERLIN) observations of the nuclear region in I Zw 1 and we also analyze the archival data from VLA and MERLIN.Our paper is organized as follows, Section 2 details the multi-band observations, data reduction, and analysis of the target I Zw 1, while Section 3 presents the results and discussions.Finally, we provide our conclusions in Section 4. Throughout this work, we adopt the standard ΛCDM cosmology with H 0 = 71 km s −1 Mpc −1 , Ω Λ = 0.73, Ω m = 0.27, and the corresponding physical scale is 1.125 pc mas −1 in I Zw 1.

VLBI observation and data reduction
We observed I Zw 1 on 2018 September 23 with 10 antennas of the Very Long Baseline Array (VLBA) and on 2020 November 17 with 19 antennas of the European VLBI Network (EVN) plus the enhanced Multi-Element Remote-Linked Interferometer Network (e-MERLIN).The VLBA observation was carried out at L-band (1.548 GHz or 1.5 GHz for short, the project code BY145), and the EVN+e-MERLIN observation was conducted at C-band (4.926 GHz or 5 GHz for short, the project code EY037), respectively.The total VLBA observing time is 2 h with a data recording rate of 2 Gbps, and the total time of the EVN+e-MERLIN observation is 8 h with a data recording rate of 4 Gbps.Both observations were performed in the phase-referencing mode, using J0056+1341 (R.A.: 00 h 56 m 14.816010 s ±0.000013 s , Dec.: +13 • 41 15.75506 ±0.00044 ) as the phase reference calibrator.
We calibrated the VLBI data in the Astronomical Image Processing System (aips), a software package developed by the National Radio Astronomy Observatory (NRAO) of U.S. (Greisen 2003), following the standard procedure.A-prior amplitude calibration was performed using the system temperatures and the antenna gain curves provided by each station.The earth orientation parameters were obtained and corrected using the measurements from the U.S. Naval Observatory database and the ionospheric dispersive delays were corrected based on a map of the total electron content provided by the GPS satellite observations4 .The opacity and parallactic angles were also corrected based on the auxiliary files attached to the data.The delay of the visibility phase and the telescope bandpass were calibrated using the bright radio source 3C 454.3.Next, we performed a global fringe-fitting on the phase-referencing calibrator, J0056+1341, by assuming a point source model to solve miscellaneous phase delays.
The phase calibrator J0056+1341 shows a core-jet structure that extends upto ∼100 mas to the north (see supplementary Figure 1 in Appendix B for their 1.5 and 5 GHz images).We performed self-calibration to the phase calibrator and obtained its CLEAN model, which was then used as the input model to re-solve the phases in aips.This operation can eliminate phase reference errors due to the jet structure.Finally, both phase and amplitude solutions obtained from the phase calibrator were applied to the target I Zw 1.The calibrated uv-data was exported to difmap (Shepherd 1997) for deconvolution.Based on a signal-to-noise ratio of ∼ 30 and ∼ 14 in 1.5 and 5 GHz in the residual map, respectively, we decide to not perform self-calibration to the target source.
We performed different deconvolution algorithms in difmap to produce radio maps, i.e.CLEAN and Gaussian model-fit.It is important to note that the solutions for visibilities are not necessarily unique when complex structures are to be handled.In order to compare the goodness of the deconvolution results, we summarised statistical parameters in Table 1.Obviously, the statistical parameters of the Gaussian model-fit are close to the CLEAN, while the CLEAN is better than the Gaussian model-fit based on the χ 2 r -value (especially at All the images are produced with natural weight and the map reference is at the Gaia position.The restoring beams are displayed as grey ellipses in the lower-left/right corner of each panel, which are 19.9 × 14.9, 11.5 × 4.67, 9.27 × 5.95 and 3.22 × 1.14 mas from panel (a) to (d), respectively.The contours are at 3σ × (−1, 1, 1.41, 2, 2.83, . . .), where positive contours are white and black solid and negative ones are red dashed.From panel (a) to (d), the image peak is 1.66, 0.88, 0.24, and 0.10 mJy/beam, respectively, and the rms noise is 0.06, 0.025, 0.008 and 0.007 mJy/beam, respectively.The circles with crosses inside indicate the corresponding Gaussian models.In panels (b) and (c), the red stars mark the Gaia optical positions.The uncertainty in the Gaia position is ∆α = 0.18 and ∆δ = 0.17 mas, which includes an astrometric excess noise error of 0.14 mas.At the redshift of I Zw 1, 1 mas corresponding to 1.125 pc.
5 GHz), which is reasonable as the emission regions are generally more complex than the representation of a few Gaussian components.The resulting CLEAN images and Gaussian model-fitted components are shown in Figure 1.The root-mean-square (rms) noise in the residual map (Column 3 of Table 1) is larger than the off-source rms noise in the map (Column 4 of Table 1), indicating that some diffuse emission cannot be recovered in the full resolution map.The off-source rms repre-sent the background noise fluctuation, whereas the rms of the residual map incorporate residual flux densities.Therefore, the loss of flux density in the images can be characterized as σr−σrms σrms (Column 5 of Table 1).Simultaneously, we also produced a uv-tapered map to reduce the weight of the long-baseline visibilities (and thus also reduce the resolution) in an attempt to recover the weak and extended emission, see panels (a) and (c) in Figure 1.

Archived VLA data
We retrieved the raw visibility data of I Zw 1 observed by the Very Large Array (VLA) from the NRAO data archive5 , including historical VLA and the newly observed Karl G. Jansky VLA (JVLA) data.Although some data have been published (see Table 2), to ensure consistency in the data reduction, we performed a manual calibration for all available datasets using the Common Astronomy Software Application (CASA v5.1.1, McMullin et al. 2007).Our data reduction followed the standard routines described in the CASA cookbook.We adopted the 'Perley-Butler 2017' flux density standard to set the overall flux density scale for the primary flux calibrator and then bootstrapped the secondary flux density calibrators and the target.For the historical VLA datasets, we determined the gain solutions using a nearby phase calibrator and transferred them to the target I Zw 1.For the JVLA datasets, we also determined antenna delay and bandpass by fringe-fitting the visibilities.For the data observed after 1998, we performed an ionosphere correction using the data obtained from the CDDIS archive.Deconvolution, self-calibration, and model-fit were performed in difmap.The final images were created using natural weight.Due to the good uv-coverage, simple emission structure, and high signalto-noise ratio (SNR>9), the VLA data allows for selfcalibration using a well-established model.For data with lower SNR, we used three times the image noise as the upper limit for the flux density.

Astrometry of the VLBI data
We measure the uncertainties of the astrometric measurements from three main origins: (1) Positional uncertainties of phase-referencing calibrators.In phasereferencing observations, the coordinates of the target are referenced to a close calibrator.The calibrator J0056+1341 was selected from the catalog rfc 2022a in Astrogeo6 with precise position with accuracies ∆α = 0.20 mas in right ascension and ∆δ = 0.44 mas in declination; (2) Astrometric accuracy of phase-referenced observations.Primarily concerning the station coordinate, Earth orientation, and troposphere parameter uncertainties, which can be measured through the formula and data from Pradel et al. (2006).This portion contribute position errors ∆α ≈ 0.26 mas and ∆δ ≈ 0.50 mas in VLBA 1.5 GHz observation and ∆α ≈ 0.27 mas and ∆δ ≈ 0.47 mas in EVN+e-MERLIN 5 GHz observation; (3) Thermal error due to the random noise (e.g.Thompson et al. 1986;Rioja et al. 2017).This uncertainty can be characterized as σ t ∼ θ B /(2×SNR), where θ B is the full-width at half maximum of the restoring beam and SNR is the signal-to-noise ratio.In this work, we take this value from difmap.
During the self-calibration process, the absolute coordinate position of the phase-referencing calibrator is lost and the brightest feature of the image is shifted to the phase center of the map.In general, due to the frequency-dependent shift in the peak of the optically thick component and the slightly different distribution of the radio emission at different resolutions, the brightest component may not be the same component from one observation to the next.This would induce a systematic offset between two images.The alignment between the images of two frequency observations can be done using an optically thin component as a reference since its position is less affected by the frequency-dependent opacity effect (e.g.Marr et al. 2001;Kovalev et al. 2008;Sokolovsky et al. 2011;Fromm et al. 2013).
In our VLBI observations, we used J0056+1341 as the phase-referencing calibrator, which has a flat-spectrum radio core.We obtained VLBI C (4.34 GHz) and X-band (7.62 GHz) data from Astrogeo.As neither the C and X-band data are self-calibrated, the core shift at C-band can be directly estimated as ∆R.A. ∼ 0.69 mas, ∆DEC.∼ −0.29 mas relative to the X-band data (see supplementary Table 3 in Appendix B).J0056+1341 also has a significant offset between C and L bands (see supplementary Figure 1 in Appendix B) in our observations, and the offset of the VLBA L-band image is determined using the optically thin component J1 of the jet, estimated as ∆R.A. = 1.742 mas and ∆DEC.= −3.228mas (see supplementary Table 3 in Appendix   ordinates for the target are also listed in supplementary Table 3 in Appendix B.

Radio Spectrum
To obtain the radio spectral index, we first checked the variability of I Zw 1. Figure 2 shows the radio flux density versus the observing epoch.The largest variability we identify is ∼8%, which is from the VLA Aarray observations at the C-band between epochs 1983 and 1995.Since there is no evidence for extreme variability on a time scales of ∼30 years, we plot the radio flux density versus the frequency in Figure 3 using all archival data.The least-square fitting gives an overall radio spectral index of −0.89 ± 0.10.From Figure 3, we can see the radio flux density changes with the collection area.The radio spectral index between 1.4 and 5 GHz, using the datasets with a similar resolution (i.e., 1.3 ∼ 1.5 and 4.3 ∼ 5.3 arcsec), is −0.69 ± 0.08 and −0.61 ± 0.02, respectively.These yields flatter spectral indices than the overall fit.Given the total flux density in Section 2.5, the spectral index between VLBA 1.5 and EVN+e-MERLIN 5 GHz is −1.06±0.13,consistent with the overall radio spectral index.
To obtain the spectral index distribution for the highresolution data observed with the VLBA and EVN+e-MERLIN, we created a spectral index map following the procedure described in Hovatta et al. (2014).Here, we  2. The blue dashed line is the model-fitting result with a power-law spectrum using all the data points presented here.The power-law slope (spectral index) is −0.89 and the blue belt shows the 95% confidence interval (0.10).The green and red dashed lines show the power-law fitting between 1.4 and 5 GHz datasets with similar size scales, i.e. 1.3∼1.5 (arcsec, red) and 4.3∼5.3(arcsec, green), respectively.The green and red belts indicate their 95% confidence intervals.
used the uv-tapered image by 0.5 at 20 Mλ at 5 GHz and restored it to match the 1.5 GHz map.Similarly, the alignment of two images is through the brightest optically thin component E1.The spectral index was calculated pixel-by-pixel between the 1.5 and 5 GHz total intensity maps.For a given frequency, pixels with an intensity less than 3σ rms were removed.The spectral index map between 1.5 and 5 GHz is shown in the left panel of Figure 4.Both Gaia position and component C are in the flat spectral region (α > −0.5, see also Figure 5.) We study the radio spectral index distribution along the jet trajectory by estimating a ridge line of the jet.We define the jet ridge as the line that connects the peaks of one-dimensional Gaussian profiles fitted to the profile (slices) of the jet brightness drawn orthogonally to the jet direction (see Vega-García et al. 2020).To obtain the ridge line, we performed a fitting to the tapered and restored (with circular beam) 5 GHz EVN+e-MERLIN image.The first slice starts from the Gaia position and a position angle of 220 • (regarding to the East direction anticlockwisely) was initially adopted for the jet direction.The step between individual slices was set to be smaller than the beam size (a measure taken in order to ensure the continuity of the brightness profiles  We estimate flux density uncertainties following the instructions described by Fomalont (1999).In this work, the integrated flux densities S i were extracted from Gaussian model-fit in difmap, where a standard deviation in model-fit was estimated for each component and considered as the fitting noise error.Additionally, we assign the standard 5 and 10% errors originating from amplitude calibration of VLBA (see VLBA Observational Status Summary 2018B8 ) and EVN (e.g.Radcliffe et al. 2018) data, respectively.
The radio brightness temperature was estimated using the formula (Ulvestad et al. 2005): where S i is the integrated flux density of each Gaussian model component in units of mJy (column 5 of Table 4); φ min and φ maj is minor-and major-axis of the Gaussian model (i.e.FWHM) or the restored beam in milli-arcsec; ν is the observing frequency in GHz (column 2 of Table 4), and z is the redshift.The resolution limit θ lim for Gaussian components can be estimated using the formula (Lobanov 2005), where θ B is the FWHM of the synthesized beam, SNR is the signal-to-noise ratio and β = 2 for the natural weight.If the fitted component size is lower than the corresponding resolution limit, then θ lim was used instead as the component size, and the component was identified as unresolved (or the radio-emitting region cannot be constrained).The resolution limits for each component are listed in Column 7 of Table 4 and the estimated 1.5 GHz and 5 GHz radio brightness temperatures are listed in Column 8 of Table 4. Since the measured component size is only the upper limit, the radio brightness temperature should be considered as the lower limit.We also estimated a total radio flux density from a uv-tapered image, which is 2.636±0.282and 0.765±0.085mJy for 1.5 and 5 GHz, respectively, and the corresponding source angular sizes are ∼50 and ∼40 mas, respectively.Panel (b) of Figure 1 shows the 1.5 GHz VLBA image, displaying a quasi-continuous emission structure elongated along the east-west direction with an extent of ∼45 parsec (pc).Panel (d) of Figure 1 shows a higherresolution (the beam FWHM is 3.22 × 1.14 mas) image obtained from the 5 GHz EVN+e-MERLIN observation.The bright components in the 1.5 GHz image are resolved into a series of knots in the 5 GHz image.Most of these components (except for components S and W1) have brightness temperatures > 10 7 K (Table 4), and the whole radio emitting structure has an overall steep spectrum (Figure 3), which favors a jet origin and is unlikely from star-forming activities (Condon et al. 1991) and thermal free-free radiation of the hot molecular disc surrounding AGNs (Gallimore et al. 1997).Furthermore, based on the identification of optical and radio cores (see below), the bilateral radio structures in the 5 GHz EVN image (panel d of Figure 1) are consistent with the assembling of approaching and receding jets.Correspondingly, the brightness asymmetry between the two branches of the bilateral structures indicates the Doppler boosting effect of relativistic jets.
The jet base of an AGN typically has a flat radio spectrum due to the synchrotron self-absorption in the optically-thick region (α > −0.5, fn. 9 ).The radio core of I Zw 1 is likely close to component C because it is located in a relatively flat spectral region (see Figure 4), has a high radio brightness temperature of log (T B /K) > 7 in both 1.5 and 5 GHz (see Table 4), and is roughly associated with the Gaia position.The 5 GHz image further resolves component C and measures a more accurate position of its peak, which slightly offsets from the Gaia position (see panel d of Figure 1 and Figure 5).Actually, it was already shown that there is a significant offset between VLBI and Gaia positions in AGNs (Petrov & Kovalev 2017;Kovalev et al. 2017;Plavin et al. 2019).Interestingly, the offset between component C and the Gaia position is consistent with the observations of Seyfert I galaxies (i.e. the upstream offsets of Gaia positions correspond to VLBI positions, see Plavin et al. 2019).Therefore, the VLBI/Gaia offset in I Zw 1 can be attributed to the dominance of the accretion disk in the optical band (the Gaia position), while the VLBI observations alternatively track the emissions of jets (Plavin et al. 2019).
Based on both the radio and optical cores, there is an obvious spectral steepening along the eastern jets, while it is less prominent in the western jet (see the right panel of Figure 4).This property resembles canonical jets in blazars and can be explained with radiative losses of the synchrotron-emitting plasma or with the evolution of the high-energy cutoff in the electron spectrum (Hovatta et al. 2014).

Is the ejection process in I Zw 1 episodic?
Component C is most likely the closest component to the radio core, however, the spectrum of C itself is too steep (as a reference, the spectral index at the peak position of component C is roughly −0.4 from the spectral index map, see Figure 5) to be from the nucleus of I Zw 1. Taking the integrated flux density of C at 1.5 and 5 GHz (see supplementary Table 4 in Appendix B), we can estimate the 1.5 − 5 GHz spectral index of component C as −0.88 ± 0.11.Recalling the offset between the VLBI position of component C and the Gaia position of the optical nucleus, the natural explanation of component C is nascent ejecta.
Interestingly, the positions of local brightness enhancements (the knots E1 and W1) are symmetric (see Figure 6) and components E1 and C seem to be embedded on a continuous background of the jet stream and mimic discrete blobs/knots.On the origin of the knots, several works envisage them as the shocks traveling along the jets (e.g.Hovatta et al. 2014), while the other works interpret them as discrete blobs from episodic ejection (e.g.Shende et al. 2019).We note that the relatively flatter spectra around the Gaussian components (C and E1) in Figure 4 do imply the shock acceleration in the discrete blobs (C and E1) in I Zw 1 (see Hovatta et al. 2014).However, it seems that the heavily knotty jets (comparing with the jets in blazars, Hovatta et al. 2014) in I Zw 1 is inconsistent with the only mild flattening in the spectrum, i.e. the spectral distribution shows only plateaux than clear bumps (see Hovatta et al. 2014).This discrepancy may be indicative of an episodic scenario.Actually, the lacks of a bright radio core and the steep spectrum in the nascent ejecta C both support the episodic ejection scenario.In observations, there is also growing evidence to show that highly accreting AGNs tend to launch episodic jets (Yao et al. 2021;Yang et al. 2022aYang et al. ,b,c, 2023)).

On the physical interpretations of the complex jets
By summing the radio flux density along the jet trajectory, i.e.only excluding the component S, the measured 5 GHz jet luminosity of I Zw 1 is log L 5GHz = 38.547± 0.003 erg s −1 .Taking the X-ray luminosity of log L 2−10keV = 43.65 erg s −1 measured in 2020 (Wilkins et al. 2021), the radio to X-ray luminosity ratio of I Zw 1 is L R /L X = 10 −5.102 , suggesting the nature of the jet is the corona ejection from the accretion disk (Yuan et al. 2009), i.e.L R /L X ∼ 10 −5 (Panessa et al. 2019;Yang et al. 2020), and consistent with the interpretation of the X-ray behavior (Gallo et al. 2007;Wilkins et al. 2017).Models and observational evidence suggest that episodic jets in both AGN and microquasar are ejected from the accretion disc corona (see Shende et al. 2019, and references therein).
The bilateral morphology and the linear size of ∼45 pc (the tapered 5 GHz image in panel c of Figure 1 and Figure 6) are reminiscent of I Zw 1 being a compact symmetric object (CSO, O'Dea & Saikia 2021).However, the lack of a spectral peak at GHz (Figure 3) alternatively implies that it belongs to compact steepspectrum (CSS) radio sources.Assuming minimum energy (approximately equipartition) conditions in synchrotron emission, we can estimate the magnetic field B min of I Zw 1 in units of Gauss (G) through the formula (e.g.Miley 1980;Patil et al. 2022) (G), (3) where S i (in mJy) is the integrated flux density of the source measured at frequency ν (in GHz) and angular size θ (in mas), α is the spectral index, z is the redshift of the source, and r is the comoving distance in Mpc.Here we take p = 0.5, the overall spectral index α = −0.89± 0.10 from ν 1 ∼ 1 GHz to ν 2 ∼ 15 GHz, a filling factor for the relativistic plasma f rl = 1 and a relative contribution of the ions to the energy a = 2.By adopting the standard ΛCDM cosmology and using the cosmology calculator provided by NED 10 , r = 245.7 Mpc in this case.Here we use the integrated flux density of 2.636 ± 0.282 mJy and the angular size of 50 mas from the tapered VLBA 1.5 GHz image, which yield an overall magnetic field B min ≈ 10 −3.6 G.The magnetic field is consistent with the typical value of e.g.CSOs with peaked-spectrum (PS, B ∼ 10 −3 G) and CSS (B ∼ 10 −4 G) radio sources (O'Dea 1998).On the other hand, the total radio power of I Zw 1 is only 10 21.8 W Hz −1 (from the 5 GHz EVN+e-MERLIN observational flux density of jets, see above), which is at least 1 order of magnitude lower than the typical sample of PS and CSS sources (O'Dea & Saikia 2021).However, the radio power of I Zw 1 remains higher than our recently discovered CSO/PS in NGC 4293 (∼ 10 20 W Hz −1 , see Yang et al. 2022c), which hold the current lower limit record of radio power.Simultaneously, the magnetic field at component E1 can be estimated as B min ≈ 10 −2.9 G by using the model-fitting results with either 1.5 GHz VLBA or 5 GHz EVN+e-MERLIN observational result.
As the promising spectral turnover of I Zw 1 must be below 1 GHz, we can estimate the lowest allowed electron lifetime via the formula (O'Dea 1998) where B is the magnetic field in G, B R 4(1 + z) 2 × 10 −6 G is the equivalent magnetic field of the microwave background, and ν p is the break frequency in GHz.For I Zw 1, using values estimated above, i.e.B ≈ 10 −3.6 G and ν p < 1 GHz, we find the electron lifetime is substantially > 10 6.3 years.Again, the estimated electron lifetime of component E1 is > 10 5.2 years.Taking component E1 as a reference (15 mas away from the core) and assuming a typical knot advance speed of 0.1 c for PS/CSS (O'Dea & Saikia 2021), the ejection time scale of component E1 is only ∼ 550 years and smaller than its electron lifetime above, therefore it suggests the knot advance speed of I Zw 1 should be substantially slower than 0.1 c.In stellar mass XRBs, discrete radio blobs are produced in the 'very high state' or 'super-Eddington state' at a time scale spanning from a few days to one year (to the best of our knowledge, see Margon & Anderson 1989;Mirabel & Rodríguez 1994;Hjellming & Rupen 1995;Fender et al. 1999;Tudose et al. 2007;Joseph et al. 2011;Miller-Jones et al. 2012, 2019;Bright et al. 2020).The longer ejection time scale in I Zw 1 seems consistent with its more massive nucleus than XRBs, i.e. the ejection time scale can be scaled through the mass of accretors, which supports the scale-invariant ejection process in both XRBs and AGNs.
Component S in the VLBA 1.5 GHz image is a real structure (panel b of Figure 1) as it is identified in the 1.5 GHz tapered map (panel (a) of Figure 1).The 5 GHz EVN image fails to detect component S possibly due to the loss of large-scale and diffuse emission.The radio emission at component S is likely responsible for the flux density deficit in 5 GHz (see supplementary Figure 3 in Appendix B).Furthermore, the southern bump in the VLA image (Figure 2) mimics component S, which requires further identification.Interestingly, component S clearly deviates from the jet trajectory, because the jet is along the East-West direction and extends up to 0.5 arcsec scale (see Figure 2).Here we interpret the bending of component S as the result of jet-medium interactions.Given that the jet direction (extending up to ∼ 0.5 parsec) is nearly aligned with the kpc-scale molecular disk (Tan et al. 2019;Shangguan et al. 2020) in I Zw 1, it's possible to support a jet-disk interaction.On the other hand, I Zw 1 hosts strong multi-scale wideangle outflows: an ultrafast wind-like outflow with a velocity of > 0.25 c obtained from the fitting of the iron Kline profile (Reeves & Braito 2019), an ionized ultraviolet gas outflow with a velocity of 1870 km s −1 (Laor et al. 1997) and X-ray outflow with velocities of ∼ 2000 km s −1 (Costantini et al. 2007;Silva et al. 2018), and a neutral gas outflow with a velocity of 45 km s −1 (Rupke et al. 2017).Given that the wide-angle outflows tend to be perpendicular to the kpc-scale molecular disk, for example, the neutral gas outflow (see Rupke et al. 2017, for their Figure 13) of I Zw 1 is along the North-South direction, therefore the component S can be a result of jet-wind collision as well.

SUMMARY
In summary, AGNs with near or super-Eddington accretion rates are often discussed as a scaled-up version of the stellar-mass black hole in the 'very high state' or 'ultraluminous state'.The observational evidence of episodic jets in I Zw 1 indicates that the analogy between AGNs and XRBs might also hold in this extreme state.Our observations imply that near or super-Eddington and extremely radio-quiet AGNs can also launch short-lived, small-scale, and weak jets.In the Supplementary section, we present further analysis of the radio emission in I Zw 1, which sheds light on the long-standing question about the origin of radio emission in radio-quiet AGNs.Such a population is important for building our understanding of jet-disc coupling in near or super-Eddington states.There is only one Galactic microquasar (SS 433) that exhibits long time-scale, super-Eddington behavior and quasicontinuous ejection (Fabrika 2004).Only a few XRBs evolve into a near/super-Eddington state and are associated with episodic jets (Revnivtsev et al. 2002;Fabrika 2004;Done et al. 2004;Miller-Jones et al. 2019), but this phase is short-lived.Finally, the radio emis-sion from a handful of (extragalactic) ultraluminous Xray sources is too weak to be detected (Kaaret et al. 2017).The time scale of the super-Eddington state in AGNs is longer than the canonical 'very high state' in XRBs according to the scaling relation, and the radio luminosity of the super-Eddington state in AGNs is higher than that of stellar-mass black holes according to the relation L R /L X = 10 −5 (Panessa et al. 2019;Yang et al. 2020), which is essential for testing the generic model of jet-disc coupling for all near/super-Eddington systems.Our findings here may indicate common features for high-Eddington AGNs, because there are similar knots/discrete radio blobs in jets and corresponding transverse/bending structures in other sources, e.g.NGC 4051 (Giroletti & Panessa 2009) with Eddington ratio 0.2 (Yuan et al. 2021) and Mrk 335 (Yao et al. 2021) with Eddington ratio 0.48-3.11(Yang et al. 2020).Further observational and theoretical studies will be necessary to establish a comprehensive understanding of the outflows of near/super-Eddington AGNs, and the results obtained from this work may contribute to this goal.This work is supported by the National Science Foundation of China (12103076, 11721303, 11991052)  Figure 3 shows I Zw 1 has a power-law radio spectrum over the entire observed frequency range.This indicates the dominance of synchrotron radiation, and it is not significantly affected by the size of collection areas (see supplementary Figure 3 in Appendix B) where the intercept between the L-band and C-band lines can be regarded as the spectral index.The overall spectral index is −0.89 ± 0.10 and there is no significant decrease at higher frequencies, indicating a continuous replenishment of fresh electrons (Morganti 2017).Interestingly, the spectral index at the jet edge farthest from the core is ∼ −1.8 (Figure 4).This is inconsistent with the overall spectral index estimated for the larger area but suggests a non-jet origin because the spectral index decreases along the jet trajectory.The large-scale flux density is dominated by diffuse radio emission with only a fraction coming from the (parsec-scale) core region (only account ∼ 30% and ∼ 47% of radio emission from 60 kpc and ∼ 1.54 kpc scale region at 1.5 GHz, respectively).The distribution of radio flux density can be fitted as S L = (0.85 ± 0.18)r 0.212±0.020and S C = (0.88 ± 0.02)r 0.131±0.005, where S L and S C is the L and C-band flux density in mJy, and r is the angular size in mas (see supplementary Figure 3 in Appendix B).The VLBA 1.5 GHz emission satisfies the flux density versus collection area distribution, while the 5 GHz flux density from our EVN+e-MERLIN observation is underestimated due to the loss of large-scale emissions.
Both star-forming activities and relativistic winds can produce large-scale radio emission (Panessa et al. 2019).Here the star-forming activities are preferentially referred to as supernovae or supernova remnants due to the power-law spectrum.Assuming all of the radio emission is from star-forming activities, we can estimate the star formation rate (SFR) from radio emission by using the SFR-radio relation (formula 3 in Yang et al. 2020).The largest SFR can be obtained from the datasets: NVSS at 1.4 GHz, AE0022 at 1.4 GHz and 4.86 GHz, and AA0048 at 14.94 GHz.These yield a SFR of ∼ 20 M yr −1 , which is similar to other estimates (Molina et al. 2021) of ∼ 26 M yr −1 , suggesting the large scale radio emission can be entirely due to star-forming activities.Whilst the SFR-radio relation is crude and we can not fully rule out the contribution of wind-like outflow, the radio-emitting wind at large scales (a few kiloparsec scales) is negligible.In addition, the radio-emitting wind is still possible at the intermediate scale (tens of parsec scale), as there are no compact supernovae and supernova remnants detected in our VLBI and e-MERLIN observations (e.g.Fenech et al. 2008).2).As I Zw 1 is not resolved in the given observations, the synthesized beams are taken to represent the collection area.The dashed lines and belts show power-law fittings and 95% confidence intervals, respectively, in the L (red) and C-band (green).

Figure 1
Figure 1.1.5 and 5 GHz VLBI images of I Zw 1.All the images are produced with natural weight and the map reference is at the Gaia position.The restoring beams are displayed as grey ellipses in the lower-left/right corner of each panel, which are 19.9 × 14.9, 11.5 × 4.67, 9.27 × 5.95 and 3.22 × 1.14 mas from panel (a) to (d), respectively.The contours are at 3σ × (−1, 1, 1.41, 2, 2.83, . . .), where positive contours are white and black solid and negative ones are red dashed.From panel (a) to (d), the image peak is 1.66, 0.88, 0.24, and 0.10 mJy/beam, respectively, and the rms noise is 0.06, 0.025, 0.008 and 0.007 mJy/beam, respectively.The circles with crosses inside indicate the corresponding Gaussian models.In panels (b) and (c), the red stars mark the Gaia optical positions.The uncertainty in the Gaia position is ∆α = 0.18 and ∆δ = 0.17 mas, which includes an astrometric excess noise error of 0.14 mas.At the redshift of I Zw 1, 1 mas corresponding to 1.125 pc.
B) relative to the EVN+e-MERLIN C-band image.For the target I Zw 1, we use the position of the brightest optically thin (α = −0.86 ± 0.07) component in the tapered EVN+e-MERLIN 5 GHz image (panel c of Figure 1) to align with the VLBA 1.5 GHz image.The peak position of the 1.5 GHz image was moved in difmap to align with the 5 GHz image by ∆R.A. = 1.33±0.38mas, ∆DEC.= −0.77± 0.69 mas, where the positional uncertainties accounts for both 1.5 and 5 GHz astrometric uncertainties of the brightest component E1.The centroid of the optical emission obtained from the second data release 7 of the Gaia mission (Gaia Collaboration et al. 2018a,b) is R.A.= 00 h 53 m 34 s .933288± 0.000012 and DEC.= +12 • 41 35 .93081± 0.00017 (J2000).This includes astrometric excess noise error of 0.14 mas.Co-

Figure 2 .
Figure 2. The radio light curves of I Zw 1 over a time interval of 37 years.The integrated radio flux densities and their uncertainties are taken from Table2, where the data with the same observing band (approximately equal central frequencies) and arrays/sub-arrays are concatenated to show the variability.

Figure 3 .
Figure 3. Wide-band radio spectrum of I Zw 1.The integrated radio flux density measurements of I Zw 1 in five radio frequency bands between 1.4 and 15 GHz are shown, where the flux density and uncertainties are taken from Table2.The blue dashed line is the model-fitting result with a power-law spectrum using all the data points presented here.The power-law slope (spectral index) is −0.89 and the blue belt shows the 95% confidence interval (0.10).The green and red dashed lines show the power-law fitting between 1.4 and 5 GHz datasets with similar size scales, i.e. 1.3∼1.5 (arcsec, red) and 4.3∼5.3(arcsec, green), respectively.The green and red belts indicate their 95% confidence intervals.

Figure 4 .
Figure 4. 1.5-5 GHz spectral index distribution of I Zw 1 on the parsec scale.Left panel: the spectral index map produced by using the naturally weighted clean map at 1.5 and 5 GHz.The region with radio flux density below 3σ was set as blank (white), i.e. the outer region of the red (for 1.5 GHz) and black (for 5 GHz) curves.Radio spectral indices within both red and black curves are reliable.The black dots and the grey line indicate the ridgeline obtained from the tapered 5 GHz EVN+e-MERLIN image.The red stars indicate the centroid positions of the Gaussian components E1 and C from the 5 GHz image.Right panel: the spectral index distribution along the ridgeline.A positive radius corresponds to positive right ascension coordinates, and vice versa.The Gaia position is set as the reference.The grey belt marks the uncertainty of the spectral indexes along the ridgeline.The red stars mark the locations of the Gaussian components E1 and C from the 5 GHz image.In both panels the blue asterisks indicate the Gaia position.

Figure 5 .
Figure 5.Comparison of positions between Gaia DR2 and VLBI component C in spectral index map.Here only the flat spectral region (α > −0.5) is shown.We take 3σ position error of component C here, where 1σ position error is estimated through the method described in Section 2.3.

Figure 6 .
Figure 6. 5 GHz flux density distribution along jet ridge line for I Zw 1.The data points with positive (blue) and negative (red) radii are the approaching and receding jets, respectively.Positional uncertainties were directly measured in fitting the ridge line in Section 2.4 and flux uncertainties are estimated by accounting for thermal noise errors and calibration uncertainties (see Section 2.5).
Note-Column 1: component name; Column 2: frequency; Column 3-4: right ascension and declination offset correspond to the Gaia DR2 position; Column 5: integrated flux density; Column 6: the angular size of components from a Gaussian model-fit; Column 7: resolution limit along the major and minor axis direction of the synthesized beam; Column 8: lower limit of the radio brightness temperature.

Figure 1 .
Figure 1.Model-fitting images of the phase calibrator J0056+1341 at 1.5 GHz (panel a and b) and 5 GHz (panel c).The images are produced using a two-dimensional Gaussian model fit with natural weights.The contours are plotted as 3σ × (−1, 1, 2, 4, 8, . . .), where σ is the root mean square (rms) noise.The white solid curves represent positive values and the red dashed curves represent negative values.The rms noise is 0.2 mJy/beam for both 1.5 and 5 GHz images.The model-fitting components are superimposed as yellow circles.The grey ellipses in the bottom left corner of each panel represent the full width at half-maximum (FWHM) of the restoring beam.The grey lines between panels c and b indicate the corresponding components without the core-shift effect, i.e. the optically thin components.

Figure 2 .
Figure 2. VLA A-array 8.4 GHz image of I Zw 1. Gaia DR2 position of I Zw 1 is set as the map center.The image peak is 0.798 mJy beam −1 .The contours are at 3σ × (−1, 1, 2, 2, 4, 8, . . . ) and 1σ = 0.038 mJy beam −1 , where positive contours are white and negative ones are red dashed.FWHM of restoring beam is 0.258 × 0.245 arcsec at 12.7 • and displayed as grey ellipses in the lower-left corner.The red arrow marks a possible southern bump.

Figure 3 .
Figure 3.The radio flux density of I Zw 1 over a collection area range from ∼0.04 to ∼50 arcsec.The integrated radio flux densities and uncertainties of I Zw 1 in L and C-band are shown (Table2).As I Zw 1 is not resolved in the given observations, the synthesized beams are taken to represent the collection area.The dashed lines and belts show power-law fittings and 95% confidence intervals, respectively, in the L (red) and C-band (green).
Grant No. 2022SKA0120102), and the science research grants from the China Manned Space Project with NO.CMSCSST-2021-A06.Scientific results from data presented in this publication are derived from the EVN project EY037 and the VLBA project BY145.The European VLBI Network (EVN) is a joint facility of independent European, African, Asian, and North American radio astronomy institutes.e-MERLIN is a National Facility operated by the University of Manchester at Jodrell Bank Observatory on behalf of STFC.The National Radio Astronomy Observatory is a facility of the National Science Foundation operated under cooperative agreement by Associated Universities, Inc.This work has made use of data from the European Space Agency (ESA) mission Gaia (https://www.cosmos.esa.int/gaia),processed by the Gaia Data Processing and Analysis Consortium (DPAC, https://www.cosmos.esa.int/web/gaia/dpac/consortium).Funding for the DPAC has been provided by national institutions, in particular the institutions participating in the Gaia Multilateral Agreement.

Table 1 .
Statistical parameters for different image deconvolution algorithms.

Table 2 .
Summary of historical observations and results for I Zw 1.

Table 4 .
Model-fitting results of the radio components detected in I Zw 1 with the VLBA 1.5 GHz and EVN+e-MERLIN 5 GHz observations.