Symbiosis of radio-emitting winds and a jet in a super-Eddington active galactic nucleus

Netherlands Super-critical accretion is the most powerful episode in nursing the black hole growth 1 and works in several types of objects 2–12 . Given that the inverse correlation between radio loud- 20 ness and Eddington ratio 12 , the super-Eddington active galactic nuclei (AGNs) hold the 21 extremely radio-quiet end of AGNs. Regarding the existence of jet in super-Eddington or 22 radio-quiet AGNs, it’s still unclear 13,14 . Years of studies indicate nearly all types of super- 23 Eddington accreting systems can launch a jet 3,5–11,15,16 , with one exception: no clear evi- dence to show jet in super-Eddington AGNs. Observations and theoretical works suggest that super-Eddington accretion can drive high-speed wind-like outﬂows , therefore through synchrotron (shocked winds) and bremsstrahlung mechanisms. How- ever, such a radio-emitting wind has not been observed in super-Eddington systems except for the Galactic micro-quasar SS 433 . In principle, high resolution very long baseline interferometry (VLBI) observation can directly map the inner structure of super-Eddington AGNs. Here, we report the discovery of the coupling of jet and radio-emitting winds in a nearby super-Eddington AGN, I Zw 1. Its parsec-scale jet exhibits a wiggling, we interpret this as a jet precession. All the features make I Zw 1 act as a scaled-up version of SS 433. The observations favour that jet can be launched in extremely radio-quiet AGNs and ubiquitous in super-Eddington accreting systems. The jet wiggling or precession can produce a large aperture-angle shock, which emphasises the jet’s contribution to gas feedback. As the jet precession was also discovered in other super-Eddington systems such as SS 433 3 and V404 5 , it is possible that there is a correlation with each other.

radio structure suggests Doppler boosting effect for the approaching jet. 66 The nucleus of I Zw 1 has an optical position (i.e. optical nucleus: R.A.= 00 h 53 m 34 s .933288± 67 0.000187, DEC.= +12 • 41 ′ 35 ′′ .93081 ± 0.00017, J2000) 28, 29 , closer to the component Wa ( Figure   68 1), which has the highest brightness temperature log T B = 12.39 K. An AGN core generally has 69 a flat radio spectrum due to the optically-thick synchrotron emission (α > −0.5 3 ). In the full 70 resolution image of C-band, the component Wa fits well with the optical nucleus position and has 71 a flat radio spectrum (Figure 2, α ∼ −0.4). Therefore, component Wa is most likely to be the 72 radio core. As we produced the ridgeline and radio spectral index map based on the full reso-73 lution L-band and tapered C-band images, only component W can be identified in the relatively 74 low-resolution images. Therefore, in the discussion of jet structures, we alternatively selected the 75 1 λ Edd ≡ L bol /L Edd ; L bol is the bolometric luminosity, L Edd = 3.2 × 10 4 (M BH /M ⊙ )L ⊙ is the Eddington luminosity, M BH is the black hole mass.
2 R = 1.3 × 10 5 (L 5 /L B ); L 5 is the radio luminosity at 5 GHz and L B is the optical luminosity of the nucleus at λ B = 4400Å; Radio-quiet AGNs are defined as R < 10.
3 S ν ∝ ν α ; S ν is the flux density and ν is the frequency 2 flattest radio spectral position (see Figure 2) along the ridgeline as the central reference point for 76 the jet trajectory, which is consistent with component W within 1-σ uncertainty. The location is 77 at ∆R.A.= 2.18 (mas, J2000), ∆DEC.= −2.18 (mas, J2000) with respect to the Gaia optical nu-78 cleus. Regarding the jet reference, the radio spectral index decreases rapidly in both extending 79 directions (right panel of Figure 2), further confirming the core environment and the jet originated 80 from there. The steepening of the spectra with increasing distance is due to the radiative losses of 81 synchrotron-emitting plasma 30 . 82 There is no theoretical consensus on whether there are jets in a super-Eddington system and 83 how the jets are produced. The observed inverse correlation between the radio loudness and the 84 Eddington ratio 12 suggests that the jet is progressively suppressed as the Eddington ratio increases, 85 mainly due to the Papaloizou-Pringle instability 31 . As the accretion rate approaches about one-86 third of the Eddington limit, the Papaloizou-Pringle instability destroys the jet-producing region 87 of the disc, suppressing the jet production. Relativistic winds arising from super-Eddington accre-88 tion disc are sometimes considered to be a secondary mechanism for further jet suppression 32, 33 .

89
However, recent studies have shown that a relatively rapidly developing magnetic field may induce 90 magneto-rotational instability 34, 35 that can stabilise the disc 36 , thus giving a chance to form the 91 jet. There is growing observational evidence for the existence of jets in super-critical accretion 92 systems, such as accreting X-ray pulsars 6 , black hole X-ray binaries 3,5,15,19 , and tidal disruption 93 events 7-11 . 94 In the 1.5 and 5 GHz images, we find an excess of transverse emission in the north-south 95 direction (denoted as N and S at 1.5 and 5 GHz). The transverse emission at 5 GHz is further 96 confirmed in the tapered image and well overlapped with the 1.5 GHz image (contours in the left 97 panel of Figure 2). We find the transverse emission seems to arise from the bright jet components 98 (W/Wa and E2), and their high brightness temperatures > 10 6 K make them unlikely to result 99 from star-forming activities. These features resemble the radio structure in SS 433 37 , in which the 100 transverse emission has a steep radio spectrum α = −1 19 , indicating that it originates from wind-101 like outflow (shocked winds) through the synchrotron process. The wind-like equatorial outflow 102 in SS 433 maintains at a projected distance of ∼ 20 mas 20, 37 from the core, corresponding to a size 103 of ∼ 10 8.5 R g , where R g is the Schwarzschild radius 38 . By comparison, the transverse emission in 104 I Zw 1 also has a steep radio spectrum (α ∼ −1, Figure 2) and the base of the transverse emission 105 maintains at a size of ∼ 10 mas, corresponding to ∼ 10 7.0 R g , which can be explained by a similar 106 synchrotron process of the relativistic electrons accelerated by the wind-driven shocks 39 .

107
The wind-like outflow in the super-Eddington regime has long been suggested by analytical 108 and numerical models 18 . It is found that the wind-like outflow is ubiquitous in black hole X-ray 109 binaries when it is in a high accretion state close to the Eddington limit 17 . There are also signs 110 of wind-like radio-emitting outflows in high accreting AGNs 40 . I Zw 1 do host strong multiscale

115
The radio emission from the wind-like outflow (shocked winds) in I Zw 1 may arise from the 116 jet collision with the massive outflow (perpendicular to an outer accretion disk) as was suggested in SS 433 44 or a parsec-scale dense molecular gas disc/torus (jet-induced wind-like outflow, when 118 jet-ISM collision, the jet plasma tends to flow laterally along the decreasing gradient), or both. The 119 jet-ISM collision interpretation is supported by the parameter η < 1 (Figure 6), which suggests 120 that a dense environment hinders the expansion of the jet plasma 45 . Further support is provided by 121 recent observations of I Zw 1 at other bands: the direction of the kilo-parsec-scale molecular gas 122 disc was found 46, 47 to be approximately along the parsec-scale jet direction here, and a radius of 123 ∼7 mas (∼8 pc) torus is also proposed 48 .

124
Furthermore, we know that jet direction is a good indicator of BH spin 49 or inner accre-125 tion disk orientation 50 , i.e. jet is either along the BH spin direction or inner accretion disk axis. The phase calibrator J0056+1341 shows a core and a jet extending ∼100 mas to the north 151 (see the 1.5 and 5 GHz images in Figure 5). We first performed self-calibration to the calibrator 152 and obtained its CLEAN model, which was then used as the input model to re-solve the phases in 153 AIPS. This operation can eliminate the phase reference errors due to the jet structure. We applied 154 the solutions obtained from the calibrator by re-running fringe fitting to the target. Next, the data 155 from the target source is exported to DIFMAP 52 for deconvolution. No self-calibration was used 156 for the target since it is too weak and resolved, with a signal-to-noise ratio of ∼ 30 and ∼ 14 at 1.5 157 and 5 GHz, respectively.

158
We performed CLEAN and model fitting algorithms to produce radio maps. The CLEAN 159 algorithm is suitable for compact and isolated emission structures, however, it has serious problems 160 in recovering a complex emission source when the uv sampling is poor, as is always common in 161 VLBI observations. In contrast, the model fitting procedure has the desirable property to take 162 into account the statistical details, and is suitable for reliable parametric representation of the 163 sparsely sampled visibility data. It is important to note that the solutions for model fitting are more 164 or less not unique when complex structures are to be handled. Here we use the delta function   Table 2), to ensure consistency in the data reduction, we performed a manual calibration for all available datasets using 195 the Common Astronomy Software Application (CASA v5.1.1) 53 . Our data reduction followed the induce an offset between two images, in the form of astrometric error. The alignment between the 220 images of the two observations can be done using an optically thin component, whose position is 221 less affected by the frequency-dependent opacity effect, as a reference.

222
For VLBI observations in both L and C bands, we used J0056+1341 as the phase-referencing 223 calibrator, which has a flat-spectrum radio core. The absolute astrometric position of J0056+1341 224 is derived from the VLBA X-band (7.6 GHz) observations. Since neither the C-and X-band data 225 from Astrogeo 7 are self-calibrated, the core-shift at C-band can be directly estimated as ∆R.A. metric excess noise error of 0.14 mas. The optical core I Zw 1 located in the radio emission region.

239
The Gaia position indicates the compact and brightest feature in the optical band, which 240 we use as the optical nucleus. In radio bands, we identify the radio core based on the compact  Figure 7).

253
Radio Spectrum To obtain the radio spectral index, we first checked the variability of I Zw 1.

262
For the high-resolution data observed with the VLBA and EVN+e-MERLIN, in order to 263 obtain its spectral index distribution, we obtained a spectral index map following the procedure 264 described in 30 . Here we used the uv-tapered image by 0.5 at 20 Mλ at 5 GHz and restored it to 265 match the 1.5 GHz map. The spectral index was calculated pixel by pixel between the 1.5 and 266 5 GHz total intensity maps. For a given frequency, pixels with an intensity less than 3σ rms were 267 removed. The spectral index map between 1.5 and 5 GHz is shown in the left panel of Figure 2.

268
Radio flux density and Brightness Temperature We estimated the uncertainties of integrated 269 flux density S i and peak flux density S p using the expression σ i = 2.5σ 2 rms + (0.01S i ) 2 55 and 270 σ p = σ 2 rms + (1.5σ rms ) 2 ≈ 1.8σ rms 56 , respectively, where σ rms is the RMS noise estimated in 271 8 https://gea.esac.esa.int/archive/ 7 a blank sky zone far away from the target source. The radio brightness temperature was estimated 272 by using the formula 57 : is only the upper limit, the radio brightness temperature should be considered as the lower limit.

284
We also estimated a total radio flux density using the tapered images at 1.5 and 5 GHz, which 285 is 2.636±0.098 and 0.614±0.047 mJy, respectively, the corresponding angular size is ∼50 and 286 ∼40 mas, respectively.

293
The total synchrontron emission from a single optically thin jet knot where the synchrotron-294 emitting plasma undergoes adiabatic expansion can be scaled 45 as where p = 1 − α and α is the spectral index of the emission (S ν ∝ ν α , ν is the frequency, R is 296 the radius of a jet knot, assuming symmetric expansion of the approaching and receding jets. By 297 taking the expansion of the plasma to be the form where r is the angular distance to the core and υ i is velocity on the sky plane, D is the distance from 299 the observer to the target, η indicates the linearity of jet expansion (R ∝ t η ), for linear expansion 300 η = 1 and for the deceleration of the expanding plasma, η < 1. We can get Then the ratio of radio flux densities as seen by the observer is where S app and S rec are the flux density for approaching and receding jet components, respectively, 303 β = υ/c is the jet speed in unit of speed of light c. t app and t rec are the times at which light leaves 304 the approaching and receding knots, respectively, and L app and L rec is their luminosities, k takes 305 value 2 in our case for a continuous jet. Here we know the velocities on the sky plane (or proper 306 motions) at approaching υ app and receding υ rec jet knots have the relation where S app and S rec are the flux densities of a corresponding pair of approaching and receding 309 knots, r(t app ) and r(t rec ) are their proper motion distances to the core. Here we fit our data with 310 the model above, η as a free parameter, and introduce a new free parameter f 0 as the flux density where f 0,app and f 0,rec are the flux density at the same distance in approaching and receding jet, 313 respectively. At the same jet direction we have The flux density scale f 0 of approaching and receding jet and the jet expansion parameter η can 315 be constrained based on the variation of the jet flux density with the distance from the core. In 316 Figure 6, we show the constraints on parameters, and we take one percent of the integrated flux 317 density as the error. Due to the limited data points and the effect of jet helical motion and colliding 318 with dense medium, for example, the radio flux density at radius ∼ 1 and ∼ 10 mas (which does 319 not take to constrain the parameters). In order to give a reasonable estimate of parameter f 0 and 320 η, we select bright components for fitting, which are marked with blue (approaching jet) and red 321 (receding jet) in Figure 6), we firstly fit the data points at the approaching jet with least-square,     Competing Interests The authors declare that they have no competing financial interests.  Table 1: Model-fitting results of the radio components detected in I Zw 1 with the VLBA Lband and EVN+e-MERLIN C-band observations. Column 1: component name; Column 2: frequency; Column 3: right ascension offset, the uncertainty is 0.47 and 0.12 mas for 1.548 and 4.926 GHz, respectively; Column 4: declination offset, the uncertainty is 1.14 and 0.31 mas for 1.548 and 4.926 GHz, respectively; Column 5: integrated flux density; Column 6: angular size of components (if given) from a Gaussian model-fit; Column 7: lower limit of the radio brightness temperature.  In the L-band image, the beam size is 11.5 × 4.67 mas and the rms noise is 0.025 mJy/beam. In the C-band image, the beam size is 3.22 × 1.14 mas and the rms noise is 0.007 mJy/beam. At both images, the red asterisks indicate the Gaia position, the uncertainty of the Gaia position is ∆α = 0.18 and ∆φ = 0.17 mas at R.A. and DEC., respectively, including an astrometric excess noise error of 0.14 mas. All the images are produced with a natural weighting and the map reference is at the Gaiaposition, the contours are at 3σ × (−1, 1, 1.41, 2, 2.83, ...). At the redshift of I Zw 1, 1 mas corresponding to 1.139 pc. The components are recognised based on the delta-model-fitted images except for the component X, which we cannot match with the model-fitted image. The region with radio flux density below 3σ was set as blank (white), i.e. the outer region of the red and white lines, respectively. Radio spectral indexes within both red and black lines are most reliable. The cyan asterisks indicate the jet ridgeline obtained from the tapered 4.92 GHz image; Right panel: the spectral index distribution along the ridgeline, the positive radius corresponding to the positive right ascension coordinates, and vice versa, and the Gaia position is set as the reference. The blue belt is the uncertainty of spectral index, in both panels the red asterisk indicates the Gaia position.

536
The jet wiggling The jets at both 1.5 and 5 GHz is strongly nonlinear, with a position angle cov-537 erage of 18 • at the approaching jet (the east branches). The non-linearity grows with the distance 538 from the core and the jet wiggling is more likely periodic than random (Figure 7), hints a helical 539 jet structure. We use the components from δ-model fitting of the tapered 5 GHz image to depict 540 the ridgeline of the jet limb (see Methods), and the ridgeline clearly shows a helical pattern. In the 541 5 GHz image (Figure 4), we also found the counter jet, which is more clear in the tapered image 542 (Figure 4), indicating the two-sided jet. Under the assumption of symmetrically approaching and 543 receding jets, we can estimate the inclination angle i to the line of sight and the jet velocity β 544 (β = υ/c, where υ is velocity, c is the speed of light, see Methods). Figure 6 shows the radio flux 545 density versus distance from the core. The radio emission at ∼ 1 and ∼ 10 mas region may be due in I Zw 1.

559
The helical jet structure in the super-Eddington accreting system I Zw 1 enhances the simi-560 larity between super-Eddington X-ray binaries such as SS 433 3 and V404 Cygni 5 , and with three 561 sources so far showing helical jets at super-Eddington accretions, we are almost coming to a con-562 sensus and a strong demand for a unified model. The most plausible mechanism in driving the 563 helical jet structures is accretion disk precession. If the axis of the accretion disc is tilted to the 564 black hole spin, the frame-dragging generated by the black hole rotating can lead to particle pre-565 cessing inside the disc, a.k.a. the Lense-Thirring effect 60 . The frame-dragging effect decreases 566 with distance and becomes negligible at large distances. Therefore an inclined viscous accretion 567 disk with respect to the spin of the black hole will produce warps in the disk which force the align- Large-scale radio emission Figure 10 shows I Zw 1 has a power-law radio spectrum over the 584 entire observed frequency range, indicating the dominance of synchrotron radiation, which does 585 not significantly affect by the collection areas, see Figure 11, where the intercept between the L- the L and C-band flux density in mJy, and r is the angular size in mas ( Figure 11).

596
Both star-forming activities and relativistic winds can produce large-scale radio emission 13 , 597 here the star-forming activities are preferentially referred to supernovae or supernova remnants due 598 to the power-law spectrum. Assuming the whole radio emission is from the star-forming activities, 599 we can estimate the star formation rate (SFR) with radio emission by using the SFR-radio relation   9×14.9 mas, the rms noise is 0.06 mJy/beam. The right panel shows the delta-model-fitted image with a beam size of 11.5 × 4.67 mas and the rms noise is 0.029 mJy/beam. In the delta-model-fitted image, the black crosses show the location of components, the red asterisks indicate the Gaia position, the uncertainty of the Gaia position is ∆α = 0.18 and ∆φ = 0.17 mas at R.A. and DEC., respectively, including an astrometric excess noise error of 0.14 mas. All the images are produced with a natural weighting and the map reference is at the Gaia position, the contours are at 3σ × (−1, 1, 1.41, 2, 2.83, ...). At the redshift of I Zw 1, 1 mas corresponding to 1.139 pc.   : Constraints on the flux density factor b and expansion factor η with the adiabatic jet expansion model for I Zw 1. The model is constrained on the radio flux density distribution along with the distance to the core, where the data points with positive and negative radius are the approaching and receding jet, respectively. The radio flux density is from clean components at 4.92 GHz, where only the bright components (marked with the blue and red circle) are used for least-square fitting, the components at radius ∼1 and ∼10 mas are not considered in the fitting. The blue and red solid line indicates best-fit results for approaching and receding jets, respectively, where the blue and red belt is the 1σ error.    Table 2, where the data with the same observing band (approximately equal central frequencies) and arrays/sub-arrays are concatenated to show the variability. Figure 10: 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 Table 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.88, the blue belt shows the 95% confidence intervals (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 (red) and 4.3∼5.3 (green), respectively, the green and red belts indicate their 95% confidence intervals, respectively. Figure 11: 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 bands are shown, which are taken from Table 2. As I Zw 1 is not resolved in the given observations, the synthesised beams are taken to represent the collection area. The dashed lines and belts show power-law fittings and 95% confidence intervals, respectively, at L (red) and C-band (green). The intervals between L and C-band flux density represent the radio spectral index.