Revisiting the Dragonfly Galaxy II. Young, radiatively e ffi cient radio-loud AGN drives massive molecular outflow in a starburst merger at z = 1 . 92

Radio-loud active galactic nuclei (RLAGNs) are a unique AGN population and were thought to be preferentially associated with supermassive black holes (SMBHs) at low accretion rates. They could impact the host galaxy evolution by expelling cold gas through the jet-mode feedback. In this work, we studied CO(6-5) line emission and continuum emission in a high-redshift radio galaxy


INTRODUCTION
RLAGNs are rare amongst all AGN populations (15 − 20%; Kellermann et al. 1989;Williams et al. 2018).At higher redshifts (1 < z < 2), there is an increasing space density of RLAGNs hosted by radiatively efficient (Eddington ratio λ Edd ≳ 3 × 10 −2 ; Merloni & Heinz 2008) SMBHs linked to high-excitation radio galaxies (see Hardcastle & Croston 2020 and references therein).These high-z galaxies hosting RLAGNs are often referred to as high-redshift radio galaxies (HzRGs) with rest-frame radio power L 500 MHz > 10 27.5 W Hz −1 (Miley & De Breuck 2008).The typical infrared (IR) luminosity of HzRGs is found to exceed 10 12 L ⊙ , in agreement with the classifications of ultra-luminous infrared galaxies (ULIRG; L IR ≥ 10 12 L ⊙ ), while some are found to be even brighter than 10 13 L ⊙ , entering hyper-LIRG (HyLIRG; L IR ≥ 10 13 L ⊙ ) regime (Drouart et al. 2014).U/HyLIRG populations are often massive, starburst systems and are thought to be triggered by galaxy mergers (Sanders & Mirabel 1996).The fact that HzRGs are also U/HyLIRGs suggests that powerful RLAGNs are often associated with major mergers, which is evident by accumulating observations (e.g., Ramos Almeida et al. 2012;Chiaberge et al. 2015;Pierce et al. 2022Pierce et al. , 2023).An investigation into the cold molecular gas -the raw fuel that feeds both the growth of SMBHs and star formation -provides insights to enhance our understanding of galaxy mergers, starburst galaxies (SBGs), and AGN/SMBHs as the basis for co-evolution in the high-z universe.
Recent studies of HzRGs find that the most powerful RLAGNs are likely to be associated with AGN with bolometric luminosity (L AGN,bol ) comparable to high-z QSOs as the SMBHs are in a fast growth phase (e.g., Nesvadba et al. 2017;Ichikawa et al. 2021).For QSOs, the accreting gas around the SMBHs usually forms an optically thick, geometrically thin accretion disc (Shakura & Sunyaev 1973;Narayan et al. 1998).The question then arises, that, a standard thin disc (0.01 ≲ λ Edd ≲ 1) is widely considered jet phobic, possibly because the large-scale magnetic field that powers the radio jets is diffused out of the accretion disc faster than being dragged inward (Lubow et al. 1994;Guilet & Ogilvie 2012) (however, recent simulations show that thin discs can sustain large-scale poloidal magnetic fields around the BH; see, e.g., Rothstein & Lovelace 2008;Liska et al. 2019).One possible explanation for the origin of the powerful jets is that these most powerful HzRGs (jet kinetic power L jet ≳ 10 47 erg s −1 and L AGN,bol ≳ 10 46 erg s −1 ) may in general host over-massive SMBHs (10 −2 ≲ M BH /M ⋆ ≲ 2 × 10 −1 ; Ichikawa et al. 2021) that lie far above the BH mass (M BH ) -stellar mass (M ⋆ ) scaling relations.The M BH estimated based on the scaling relations is underestimated, thus the Eddington ratio corresponding to a standard thin disc is an overestimation.These SMBHs are actually surrounded by optically thin, geometrically thick sub-Eddington accretion discs (Ichikawa et al. 2021).However, using rest-frame optical emission lines (Hα, Hβ, Mg II, and C IV), Poitevineau et al. (2023) estimated the BH masses of a sample of RLAGNs at 0.3 < z < 4 with L jet < 10 47 erg s −1 .They found a good agreement of the M BH -M ⋆ scaling relations of RLAGNs and other AGN populations in samples across different redshifts.Therefore, these HzRGs seem to follow the scaling relations statistically.This raises another possibility, that is, many HzRGs host SMBHs accreting above the Eddington limit, and thus the accretion discs are both optically and geometrically thick, capable of powering the jets (Tchekhovskoy 2015).
AGN feedback has been an important topic to study AGN-host galaxy co-evolution because it may lead to an outflow (compression) of the gas to suppress (enhance) the host galaxy star formation (e.g., Cano-Díaz et al. 2012;Cicone et al. 2014;Shin et al. 2019;Duncan et al. 2023).Outflows in ionized and/or molecular forms are often observed in galaxies hosting AGNs with L AGN,bol ≳ 10 44 erg s −1 from local to high redshifts (e.g., Fiore et al. 2017;Fluetsch et al. 2019;Bischetti et al. 2019), as well as those with lower AGN luminosities (Jarvis et al. 2021).These outflows are often ascribed to the radiative-mode AGN feedback no matter it is through the radiation pressure on dust or through shocks generated from inner AGN winds (Costa et al. 2014;Ishibashi et al. 2018).This radiative-mode AGN feedback has the potential to remove a significant fraction of molecular gas from the AGN host galaxy (Zubovas & King 2012).Without molecular gas as the fuel to support ongoing star-forming activities, the SFR can decrease and these galaxies will finally become quiescent.In addition to the radiative-mode, the jet (kinetic)mode AGN feedback, through the couplings between jets and the host galaxy interstellar medium (ISM), has also been found to be efficient in accelerating gas with colossal energy injections to power the outflow (Wagner et al. 2012;Mukherjee et al. 2018;Meenakshi et al. 2022b).Both simulations and observations find that the jetmode feedback can lead to the quenching of host galaxies (Nesvadba et al. 2006;Huško et al. 2024).Because of the simultaneously high L jet and L AGN,bol in powerful HzRGs, jet-and radiative-mode feedback co-exist in the host galaxies (Huško et al. 2024).This makes the outflow mechanism mysterious: which is the main driver of the outflow?Additionally, do outflows in these HzRGs show distinctly different molecular outflow properties when compared to QSOs and low-z AGNs?These questions remain unclear and require accumulating observational studies to investigate.
RLAGNs at an early evolution phase, at a time of radio jets remaining powerful and efficiently injecting their energy into the ISM (Meenakshi et al. 2022b), are good targets to investigate the jet-mode feedback.Moreover, by investigating the morphology of RLAGN host galaxies at z > 1 using HST WFPC3 images, Chiaberge et al. (2015) found that 92% of these RL galaxies show recent or ongoing merging events.This finding is suggestive of possible mergertriggered AGN activities.Therefore, young RLAGNs at high-z serve as ideal proxies to scrutinize the co-evolution of AGN and its host galaxy through processes like galaxy interaction/merging and AGN feedback since the stochastic gas inflows accompanied with these events may both fuel star formation and trigger AGN activities (Hardcastle & Croston 2020;Stemo et al. 2021).
In this work, we study a HzRG -MRC 0152-209 (L 147 MHz = 3.2(±1.0)× 10 28 W Hz −1 ), named Dragonfly galaxy, which has starbursts (SFR ∼ 3000 M ⊙ yr −1 ; Drouart et al. 2014).It is a HyLIRG (L IR ∼ 2 × 10 13 L ⊙ ) and a major merger comprising three components: the North-West (NW) galaxy, the South-East (SE) galaxy, and a possible companion galaxy.Its double radio hotspots were identified by the Very Large Array (VLA) at 4.7 and 8.2 GHz (Pentericci et al. 2000).Zhong et al. (2023) (referred to as Paper I) further investigated the radio hotspots combining high-resolution and highfrequency VLA and Atacama Large Millimeter/submillimeter Array (ALMA) observations.They found that the Dragonfly galaxy can be classified as a Compact-Steep-Spectrum source.The radio hotspots have an age of ∼ 2 × 10 5 yr, in line with the typical order of magnitude of young RLAGNs.We further present high-angular resolution ALMA observations (0.08 ′′ for Cycle 4 and 0.02 ′′ for Cycle 6) of CO(6-5) line emission from the Dragonfly galaxy to investigate the molecular gas within sub-kpc regions.Our new study provides a unique view of the molecular gas distribution among this merging system, revealing details of the AGN-driven outflow.

OBSERVATIONS
2.1 ALMA Band 6 ALMA Cycle 4 observations (Project ID: 2016.1.01417.S, PI: Bjorn Emonts) were conducted on 9 and 17 August 2017 for 1.2 hours onsource time with 45 antennas and baselines of 12-m array up to ∼ 3.6 km.ALMA Cycle 6 observations (Project ID: 2018.1.00293.S, PI: Bjorn Emonts) were conducted on 23 June 2019 for 2.2 hours onsource time with 48 antennas and baselines of the 12-m array up to 11.5 km.For both observations, there are four spectral windows configured to cover two 3.75 GHz bands, one of which includes 235.83−239.58GHz to observe the redshifted CO(6-5) line emission (ν rest = 691.47GHz) and another includes 251.20−254.95GHz such that only continuum is observed.The redshifted frequency is estimated based on z = 1.92 determined by the observation of CO(1-0) line emission (Emonts et al. 2011).The data calibrations were performed via the ALMA pipeline built within CASA (Common Astronomy Software Applications; McMullin et al. 2007;CASA Team et al. 2022) version 4.7.2 for Cycle 4 and 5.4.0 for Cycle 6, respectively, by running the calibration scripts supplied with the data by the North American ALMA Science Center (NAASC).
Prior to imaging the CO(6-5) line emission, we subtracted the continuum emission in the uv-plane.To do so, we first flagged the channels that include the real line emission, as well as those that include pseudo-line emissions because of the strong atmospheric absorption between 237.15 GHz and 239.1 GHz (private communication with Hiroshi Nagai through ALMA helpdesk, Jul 27, 2022).We then estimated the continuum emission by a linear function to the line-free channels and subtracted it in the uv-plane by using the task uvcontsub.We adopted the same methodology for Cycle 4 and 6 observations to image the line and continuum emissions in this work.First, we created a dirty image without any clean to calculate the root-meansquare noise (σ) under the 'briggs' weighting with a robustness parameter +0.5.Then, we cleaned the image non-interactively by setting 3σ as the stop threshold of the cleaning.We used the 'högbom' deconvolution algorithm to produce the restored image and applied a primary beam correction on the restored image.For the CO(6-5) line emission, the channel width was set to be 20 km s −1 for both Cycle 4 and 6 observations, and the reference frequency was set as 236.69GHz, corresponding to the line at z = 1.9214.We further imaged the CO(6-5) line emission using concatenated Cycle 4 and 6 (C46 hereafter) data.A robustness parameter of +2.0 was chosen for the 'briggs' weighting and a spectral resolution of 15 km s −1 was adopted.We chose this natural weighting to balance the beam size and noise level.The basic information and properties of the cleaned images of line and continuum emissions are summarized in Table .1.
The properties of the line profiles and their associated physical properties based on the Cycle 4 and 6 combined data are not presented in this paper due to an abnormal elevation in flux densities.This was considered a software issue of CASA, which has been reported to the CASA development team (private communication with Hiroshi Nagai through ALMA helpdesk, Mar 25, 2022).This abnormal elevation only influences the combined dataset and has no impact on Cycle 4 and 6 observations individually.

VLA Band Q
The VLA observations reported in this work were conducted in BnA-configuration centered at 44 GHz with an effective bandwidth of 7.5 GHz.The total on-source time is 42 min.
The observations were calibrated by requesting pipeline calibrations through the NRAO Science Helpdesk.We imaged the radio continuum emission using the 'högbom' deconvolution algorithm.A 'briggs' weighting scheme with a robustness parameter +0.5 was chosen to multiply the visibility value during the gridding of the dirty image.To create a clean image, we put a mask on the strongest signal and manually iterated until the strongest signal reached the noise level in order to circumvent artefacts as a result of overcleaning Due to the low signal-to-noise ratio of each spectral window, no self-calibration in any dataset can be applied, leaving lowlevel sidelobe contamination on the clean images.

Morphology
In panels (a) and (c) of Fig. 1, we show the images of the dust continuum, and in panels (b) and (d), we show the integrated intensity (mom0) map of CO(6-5) line emission integrated over −500 to +250 km s −1 , for both Cycle 4 and 6 observations.All line and continuum images and contours are overlaid directly based on the world coordinate system (WCS) after the clean procedure without any manual manipulation to associate the features observed in different bands.A systematic offset in the WCS of VLA may exist because of the astrometry worsened by the very extended baseline configurations and atmospheric conditions2 .Under a typical condition3 , the VLA and ALMA observations are consistent in positions within 2σ (see also Paper I, Zhong et al. 2023) NW and SE galaxies show individual, putative disc structures in both line and continuum emissions of Cycle 4 observations.The molecular tidal bridge identified in ALMA Cycle 2 observations between SE and NW galaxies is not observed in both Cycle 4 and 6 observations because of its diffuseness and low molecular gas mass (see Emonts et al. 2015b;Lebowitz et al. 2023 for details).In addition to these two galaxies, the Cycle 4 dust continuum has identified a companion component (see panel (a) in Fig. 1), which is undetected in both Cycle 4 and 6 mom0 maps.In Cycle 2 observations, this companion was identified in CO(6-5) line emission while missing in the dust continuum, and argued to be a small companion galaxy (Emonts et al. 2015b).
The highest resolution Cycle 6 observations have revealed the detailed morphology of the Dragonfly galaxy.The CO(6-5) line emission of the NW galaxy is still disc-like but is highly compact with a 2D Gaussian component of 0.091 ′′ (±0.01 ′′ ) × 0.083 ′′ (±0.01 ′′ ), corresponding to a radius of R NW ∼ 0.8 kpc, after being deconvolved from the beam.This size is more compact than R NW ≳ 1 kpc based on Cycle 2 and 4 observations.To investigate whether such compactness is attributed to the non-detection of a significant fraction of extended molecular gas in Cycle 6 observations, we compare the integrated intensity I CO(6−5) calculated in three cycles (see §3.2 and the rightmost panel in Fig. 1).In Cycle 6, the calculated I CO(6-5) = 1.79 ± 0.13 Jy • km s −1 is consistent with I CO(6-5) = 2.0 ± 0.2 Jy • km s −1 calculated by Emonts et al. (2015b) in Cycle 2 and I CO(6-5) = 1.98 ± 0.13 Jy • km s −1 calculated in Cycle 4 within 1σ uncertainty.This suggests that the bulk of molecular gas is concentrated within a sub-kpc scale disc, though there can be ∼ 10 per cent of the diffuse molecular gas missed by the high-resolution beam.In the dust continuum, the NW galaxy has its flux density concentrated within a radius of ∼ 0.07 ′′ (≈ 0.6 kpc), agreeing well with the  [3,5,7,9,11,13]×σ.Throughout all images, the world coordinates are identical.The black contours indicate the VLA 44 GHz observation in BnA-configuration (σ = 18 µJy/beam) with contour levels of [5,7,20,40,60,80] × σ with the detection coinciding with the NW galaxy indicating the radio core and the detection adjacent to the SE galaxy indicating the SE hotspot.The yellow dashed line indicates the jet axis that links the double radio hotspot (see also Fig. 1 in Paper I; Zhong et al. 2023).The beam size of the corresponding image is marked in the bottom left corner.In panels (e) -(h), we show line profiles extracted from a circular aperture that encloses contours defined by 3σ.molecular gas distribution, but shows diffuse emissions extended to the east.As discussed in §7.1.2,this extended feature may originate from the interactions between the two galaxies.
The CO(6-5) line emission distribution of the SE galaxy in Cycle 6 significantly differs from that observed in Cycle 2 and 4 observations, and the zoom-in is shown in Fig. 2. The main feature of the SE galaxy is a linear structure (labeled by #Main) elongated along the southeast to northwest direction.To the north of the centroid of the SE galaxy, there are several clumps of CO(6-5) line emission (white contours in region #Tidal), which are in agreement with the distribution of the dust continuum.These clumps viewed at ∼ 100 pc-scale form a filamentary structure on a ∼ 600 pc-scale and may come into existence through the interactions between two galaxies (see §7.1.2for a detailed discussion).On the opposite side of these clumps, there is no line emission detected.An equivalent lack of emission is seen in the dust continuum as well, save for a sub-component (see next paragraph).Additionally, in HST NICMOS F160W imaging (Fig. C1), the SE galaxy is associated with a tidal tail extended over 10 kpc.Such a long tidal tail is commonly observed in wet mergers whose bulge separation is smaller than 5 kpc (Ren et al. 2020) due to the tidal stripping (Mo et al. 2010), that is, the tidal force strips a significant fraction of the molecular gas away from the galaxy.
The fact that the double radio hotspots identified in VLA 4.7, 8.2, and 44 GHz observations are symmetric relative to the NW galaxy suggests that the RLAGN resides in the NW galaxy (see Fig. 1 in Paper I; Zhong et al. 2023).This is supported by the identification of the radio core that overlaps with the NW galaxy (see Fig. 1), though it is offset from the centroid of ALMA Cycle 6 dust continuum by 0.04 ′′ .The 237 GHz continuum at 0.024 ′′ -resolution further reveals a sub-component adjacent to and offset from the bulk of SE galaxy dust continuum by ∼ 0.15 ′′ , corresponding to ∼ 1.3 kpc (Fig. 2).This sub-component has its location coincided with the SE hotspot with an offset of merely 0.008 ′′ .It is then argued to be the radio hotspot originating from the interaction between the medium and radio jet launched by the AGN and its flux density is dominated by the synchrotron radiation (Paper I, Zhong et al. 2023).If the medium at play is the ISM of the SE galaxy, the cool gas can be either blown away by the kinetic energy or heated up to a higher temperature by the thermal energy of the radio jet.

Line Profiles
In the third column of Fig. 1, we show the spectra of NW and SE galaxies observed in Cycle 4 and 6.For both NW and SE galaxies, the line profiles are extracted using a circular aperture that encloses the spatial area defined by 3σ-level contours.The physical properties measured and derived from the line profiles are listed in Table .2.
The NW galaxy shows a clear line splitting into two velocity components and a broadening of the main component, and thus the line profile can be fitted by two Gaussian components with velocity center (v cen ), amplitude, and full width at half maximum (FWHM) as fitting parameters.In Cycle 4, the main component (NW -2 in Table.2) has v cen = −17 ± 9 km s −1 and FWHM = 310 ± 22 km s −1 , consistent with the measurements of v cen = −30 ± 10 km s −1 and FWHM = 360 ± 20 km s −1 in low-resolution Cycle 2 observations (Emonts et al. 2015b).In Cycle 6, we first applied Hanning smooth to the spectrum and then adopted a 2-component fitting by fixing v cen of the main component (NW -2) to -20 km s −1 found in both Cycle 2 and 4 observations within 1σ uncertainty.Otherwise, given the low signal-to-noise ratio, the fitting algorithm returned a third component with uncertainties exceeding 100%.The resultant Cycle 6 NW -2 has an FWHM of 414±32 km s −1 , larger than that of Cycle 4 but consistent with Cycle 2 within 1σ uncertainty.The blueshifted component (NW -1) is offset from the main component by at least 300 km s −1 in all observations.It has a broader FWHM in Cycle 6 than in Cycle 4 but almost the same amplitude, leading to a slightly higher I CO(6−5) by 0.08 Jy • km s −1 .As we will discuss in §7.3.4,this blueshifted component (NW -1) may represent the outflow driven by the expanding bubble attributed to the radio jet in a 3D shell-like geometry (e.g., García-Burillo et al. 2019).
The SE galaxy has its line profile decomposable into three Gaussian components in Cycle 4, whereas large uncertainties are left, and the sum of the integrated intensities I CO(6−5),SE = 0.86 ± 0.14 Jy/beam • km s −1 lies far below the value of 1.4 ± 0.2 Jy/beam • km s −1 observed in Cycle 2 data by (Emonts et al. 2015b).In Cycle 6, I CO(6−5),SE is 1.3 ± 0.37 Jy/beam • km s −1 , showing no significant deviations from that in the literature, though the SE galaxy shows a much more complex molecular gas distribution in Cycle 6.Additionally, in Cycle 2 and 4 observations, the brightest component has a center of ∼ 60 km s −1 (SE -2 of Cycle 4 in Table .2).However, the blueshifted component (SE -1 of Cycle 6 in Table .2) has its I CO(6−5) dominant over the other two components in Cycle 6 observations.To understand the origin of such a change in the brightest component, we zoom in on the SE galaxy and divide its mom0 map into two regions: #Main and #Tidal, as shown in the left column of Fig. 2. #Main indicates the main structure in the SE galaxy and the region #Tidal is named after its possible tidal origin to be discussed in §7.1.The spectra extracted from these regions are shown in the right column of Fig. 2, and the measured properties are listed in Table .2. The #Main region shows a bimodal distribution of the line profile characteristic of a rotating structure.The #Tidal region contains ∼ 30 per cent of the molecular gas in the SE galaxy and is globally blueshifted, which significantly contributes to the blueshifted SE -1 in Cycle 6.

Molecular gas mass
We calculated the integrated source brightness temperature from CO(6-5) line intensity using the following equation (Carilli & Walter 2013): where S CO(6−5) ∆ν is the integrated intensity of the CO(6-5) line in Jy • km s −1 , D L is the luminosity distance in Mpc, and ν obs is the observed frequency of the CO(6-5) line in GHz.To estimate the molecular gas mass, L ′ CO(6-5) has to be converted to L ′ CO(1−0) via the line intensity ratio r 65/10 = S CO(6−5) ∆ν CO(6−5) /S CO(1−0) ∆ν CO(1−0) .And the molecular gas can then be calculated using (Carilli & Walter 2013): where α CO is the CO-to-H 2 conversion factor.An α CO ∼ 0.8 M ⊙ (K km s −1 pc 2 ) −1 (Downes & Solomon 1998) found in the starburst nuclei of ULIRGs on scales < 1 kpc is adopted as a conservative estimation.
The CO(1-0) line emission was observed by the Australia Telescope Compact Array (ATCA) with a beam size of 4.0 ′′ × 1.3 ′′ , incapable of resolving two galaxies (Emonts et al. 2015a).By tapering CO(6-5) observation from Cycle 2 to the same resolution as the CO(1-0) data, Emonts et al. (2015b) compared CO(6-5) and CO(1-0) spectra extracted from the same region covering the entire Dragonfly galaxy and found that the bulk molecular gas component has a line intensity ratio r 65/10 ∼ 13.Scaling the CO(1-0) spectrum by 13 and subtracting it from the CO(6-5) spectrum, Emonts et al. (2015b) found large CO(6-5) residuals corresponding to the high-excitation gas at the blueshifted side with v ≤ −200 km s −1 .Including 2σ measurement uncertainties as an upper limit for I CO(1−0) , Emonts et al. (2015b) loosely constrained 17 ≤ r 65/10(blue) ≤ 36 for the blueshifted, high-excitation components.We adopt this high ratio for molecular gas mass estimation of the blueshifted components and the corresponding M H 2 are listed in Table .2.

Optical-to-FIR SED
The electromagnetic radiation emitted by galaxies at multiwavelength provides proxies to investigate the formation and evolution of the galaxy.Fitting the observed spectral energy distribution (SED) with a combination of various templates allows us to disentangle the galaxy emission complexity and compute both host galaxy and AGN physical properties.We performed SED fittings using CIGALE (Code Investigating GALaxy Emission; Burgarella et al. 2005;Noll et al. 2009;Boquien et al. 2019).Based on the new photometric data from optical-to-near IR (NIR) collected from the Dark Energy Camera (DECam) constructed for the Dark Energy Survey (DES; Abbott et al. 2018Abbott et al. , 2021) ) in addition to the existing IR photometry with Spitzer and Herschel, we re-examine the physical properties including SFR and M ⋆ and compare them with those reported in the literature (Drouart et al. 2014;De Breuck et al. 2010;Falkendal et al. 2019).All photometric data used for the SED fitting are summarized in Table .A1.We note that, since there are no spatially resolved optical-to-far IR (FIR) data available, the SED fitting thus can only compute the physical properties of the entire system, including NW and SE galaxies.
Recent studies of AGN host galaxies and LIRG populations through SED fitting favor a star formation history (SFH) with a recent burst (Toba et al. 2021;Dey et al. 2022;Georgantopoulos et al. 2023).This is naturally expected for galaxy mergers, especially latestage mergers that have small bulge separations.In these systems, the gas inflow -attributed to angular momentum removal by the gravitational torque -towards the galaxy center may enhance the SFR within this central region (Moreno et al. 2015(Moreno et al. , 2021;;McElroy et al. 2022;Shah et al. 2022).Therefore, we assumed a delayed SFH model with an additional burst for the host galaxy (Ciesla et al. 2015).In addition to the e-folding time (τ main ) and age (t main ) of the main stellar population, τ burst , t burst , and a mass fraction of the late burst population ( f burst ) can be parameterized.We adopted Bruzual & Charlot (2003) stellar population synthesis library to model the stellar emission with a solar metallicity 0.02, Chabrier initial mass function (Chabrier 2003), and a standard model for nebular emission (Inoue 2011).Although AGNs would enrich their environments, resulting in higher chemical abundance, we consider such an enrichment negligible in this HyLIRG (Zubovas et al. 2013;Taylor & Kobayashi 2015).The dust emission was modeled adopting a dl2014 template (Draine et al. 2007(Draine et al. , 2014)).This template includes a parameter q PAH which optimizes the mass fraction of the polycyclic aromatic hydrocarbon (PAH) emission that dominates the MIR emission, a radiation field U min that models the diffuse dust emission, and a power-law distribution of the dust corresponding to the star-forming regions with index α.The α is defined in a way such that dM d (U)/dU ∝ U −α , where M d is the dust mass and U is the radiation field intensity.The initial estimations of these parameters are referred from recent studies of AGN host galaxies (Yamada et al. 2023;Buat et al. 2021).
UV and optical emissions from the AGN will be absorbed and then re-emitted at longer wavelengths owing to the existence of the obscuring structure, i.e., the dusty torus.To simulate this reddening effect, we adopted SKIRTOR, a clumpy two-phase torus model derived from a modern radiative-transfer method, assuming that the dusty torus is made of dusty clumps rather than a smooth structure (Stalevski et al. 2012(Stalevski et al. , 2016)).This model depends on several parameters including the average edge-on optical depth at 9.7 µm, the half-opening angle of the dust-free cone, inclination (i, 0 • for face-on and 90 • for edge-on), and the AGN fraction defined as the fraction of AGN IR luminosity to total (AGN+host) IR luminosity.
The models described above and the values of corresponding parameters are summarized in Table .3.

Radio SED
In the radio regime, CIGALE models the flux densities attributed to different mechanisms, including synchrotron radiation from AGN and star formation and nebular continuum emission that involves mainly free-free emission (Boquien et al. 2019).There are four parameters controlling the model of the synchrotron radiation: the FIRto-radio correlation parameter q IR , radio loudness R AGN defined as the ratio of L 5GHz /L 2500Å where L 2500Å reflects the AGN disc luminosity, and the spectral index of star formation-(α SF ) and AGNrelated synchrotron radiation α AGN .
In normal star-forming galaxies (SFGs), the synchrotron radiation is primarily produced by the electrons accelerated in supernova remnants.In the Milky Way and many other SFGs, the overall slope of the observed radio SED has been known to have a mean value of α SF ≈ −0.75 at ν ≈ 1 GHz in a negative convention (S ∝ ν α ) (Condon & Ransom 2016; Klein et al. 2018).Hence, we fix α SF to -0.75 for the synchrotron radiation relevant to the star formation.
In the Dragonfly galaxy, the bulk of the synchrotron radiation originates from the SE hotspot and its associated diffuse radio lobe and can be significantly enhanced (see §6.2).Therefore, the observed radio emission no longer traces the intrinsic AGN activities but the physical conditions in situ.Thanks to VLA 44 GHz observations, the radio core associated with the RLAGN in the NW galaxy has been identified (Fig. 1).Such a core component generally has a spectral slope flatter than those of radio hotspots and lobes because of synchrotron self-absorption and/or free-free absorption arising from the dense environment.We adopt a spectral index of α AGN = 0, which is a commonly observed value in flat-spectrum radio quasars that are core-dominated populations (Chen et al. 2009).R AGN and q IR are fixed at 10 and 2.3, respectively.
So far, the SED fittings for powerful radio galaxies primarily make use of S 1.4GHz without spatially decomposing the emission source and considering enhancements on the flux densities.The applications of q IR and R AGN on the radio fitting are highly motivated in this manner.Therefore, a consideration of only the radio core, excluding radio hotspots and lobes, for the Dragonfly galaxy may lead to a systematic bias.We then perform a fitting using 1.4 GHz data instead of the 44 GHz data assuming a spectral index of α AGN = 1.0 as a comparison.In this fitting, we fixed R AGN = 1000 and q IR = 0.
The spatially resolved 44 GHz flux density of the radio core that overlaps with the NW galaxy is used for the fitting.
b The unresolved total 1.4 GHz flux density, including the radio core, hotspots, and lobes is used for the fitting.

Fitting Results
We have performed three fits to the SED: firstly based on S 1.4GHz of the total radio emission, then based on S 44GHz,core of the radio core, and finally without radio data.We summarize the best-fitting results in Table .4 and plot the modeled SEDs in Fig. 3.The stellar mass of NW plus SE galaxies is M ⋆,NW+SE ∼ (9 − 11) × 10 10 M ⊙ with an instantaneous SFR of ∼ 2000 − 3000 M ⊙ yr −1 .For all fittings, the instantaneous SFR is close to the SFR averaged over 10 Myr but beyond that averaged over 100 Myr, suggesting recent starbursts.We performed a test without an additional starburst by fixing the starburst mass fraction to 0. Different from the fiducial fitting with the recent burst, the test fitting includes only a single population that forms within 100 Myr, still suggesting recent star-forming activities.

Comparisons with literature
Using Spitzer and Hershcel data (3-500 µm), Drouart et al. 2014 decomposed the IR luminosity into AGN and starburst components by simultaneously calculating AGN and host galaxies contributions to the SED.They found L IR SB = 1.72 × 10 13 L ⊙ attributed to star-forming activities and estimated SFR ∼ 3000 M ⊙ yr −1 based on the SFR and IR luminosity relation ( Kennicutt 1998).Using a non-starburst SED template and combining low-frequency VLA observations at 1.4, 4.7, and 8.2 GHz, Falkendal et al. (2019) found SFR ∼ 1900 M ⊙ yr −1 .Our SED fitting analysis involving an SFH with an additional starburst finds SFR ∼ 2100 − 2700 M ⊙ yr −1 that agrees with Drouart et al. (2014) and slightly higher than Falkendal et al. (2019).The stellar mass of our fitting is merely one-fifth the value of M ⋆ ∼ 5×10 11 M ⊙ reported in De Breuck et al. (2010).Such a difference could be explained by the choice of elliptical galaxy templates instead of a star-forming one in the SED fitting performed by De Breuck et al. (2010), the lack of IR data at λ ≥ 70 µm, and the lack of optical photometry.An elliptical galaxy SED dominated by old stellar populations is a reasonable choice for low-redshift populations since powerful radio galaxies in the nearby Universe are found to preferentially reside in massive, early-type galaxies (ETGs), typical of the final evolution stage of an AGN (Hopkins et al. 2008;Shankar et al. 2020).However, RLAGNs at high redshifts are often hosted by mergers that are often accompanied by intense star-forming activities, thus the choice of SED templates of ETGs is not always proper for HzRGs.
The AGN has an IR contribution of ∼ 0.1 to the total IR luminosity given by our fitting including 44 GHz data, significantly smaller than the value of ∼ 0.4 given by Falkendal et al. (2019).Our fitting including the observed 1.4 GHz flux density, as well as the one with-out radio data, returns a similar AGN fraction of 0.3 compared with ∼ 0.2 given by Drouart et al. (2014) and ∼ 0.4 given by Falkendal et al. (2019).These results are naturally expected since when the radio data is included in our fitting, the AGN disc luminosity depends on the normalization of the luminosity at 2500 Å using R AGN .The 1.4 GHz data can obviously lead to a higher luminosity at 2500 Å than 44 GHz data does because of a much larger radio loudness.The AGN component dominates over the dust component up to ∼ 80 µm in the 1.4 GHz case and merely up to ∼ 20 µm in the 44 GHz case.Accordingly, the AGN contributes more to the total IR luminosity, resulting in a larger f AGN .

Caveats to fitting results
Since there is no spatially resolved optical-to-FIR data available for the SED fitting, the reported physical properties regarding the AGN may not be the proxy to the RLAGN that resides in the NW galaxy solely.If the SE galaxy also has an AGN, though we find no signature yet, the AGN component presented here is the sum of the NW AGN and the potential SE AGN.
The fitting of synchrotron emission does not include detailed physical processes -how jets interact with the medium in situ -but depends on the empirical q IR and R AGN for radio galaxies with lobes and/or hotspots.And the finally derived L AGN,bol will be regulated by the R AGN when radio data is included.Therefore, a fitting based on the radio core but neglecting radiation from radio hotspots/lobes requires values of q IR and R AGN for weaker AGNs, putting a lower limit on the L AGN,bol in our target.On the other hand, the imbalanced flux densities between SE and NW hotspots at 1.4 GHz are suggestive of possibly enhanced synchrotron radiation due to the environmental difference or the Doppler boosting effect (see Paper I, Zhong et al. 2023 and §6.2).Then, the fitting requires values of q IR and R AGN for more powerful AGNs, which may result in overestimated L AGN,bol .After removing the radio data, we find f AGN and L AGN,bol consistent with the results based on 1.4 GHz.This may suggest that the application of q IR and R AGN is reliable even when dealing with powerful radio galaxies with biased flux densities.However, we still call for caution when trying to retrieve L AGN,bol for these HzRGs.
In conclusion, the Dragonfly galaxy is a massive, starburst one with M ⋆ ∼ (0.9 − 1.1) × 10 11 M ⊙ and an SFR of ∼ 2000 − 3000 M ⊙ yr −1 .The RLAGN reaches a QSO-level luminosity of L AGN,bol ∼ (0.9 − 3) × 10 46 erg s −1 .We can take the fitting using 44 GHz as the lower limit for AGN properties, which sufficiently suggests that the RLAGN is radiatively efficient (λ Edd > 3 × 10 −2 ; see §6.1 for detailed discussions).In §7.2 and §7.3, we will further , where the degree of freedom is not involved.We re-calculated χ 2 red considering the number of fitted free parameters.
discuss the AGN activities by investigating the Eddington ratio and outflow properties.

Modeling with 3D Barolo
We investigated the molecular gas kinematics by fitting a kinematic tilted-ring model to the CO(6-5) line emission data cube using 3D Barolo (Di Teodoro & Fraternali 2015).3D Barolo assumes that the rotating gaseous disc of a galaxy is made up of a series of concentric rings, where each ring is described by geometrical and kinematical parameters, including spatial center (x 0 , y 0 ), systemic velocity (V sys ), position angle (PA), inclination angle (i), scale height of the disc (z 0 ), rotation velocity (V rot ), and velocity dispersion (σ V ).
The modelings are based on the C46 dataset.We adopted 8 rings for the NW galaxy and 9 rings for the SE galaxy.We chose a radial separation of 1/5 along the minor axis of the beam (Ramos Almeida et al. 2022).We tested modelings with larger separations up to beam minor /2 and found similar results.The outermost ring of each galaxy was limited to be smaller than the aperture size that includes 3σ contours in the mom0 map to exclude possible sub-structures such as tidal tails/spiral arms.None the less, as we will see in §5.3 and Fig. 4, these features still contaminate the velocity field on the redshifted side of the NW galaxy.We set the scale height as 0 adopting a thin disc approximation where the effect of any vertical structure is neglected (Ramos Almeida et al. 2022;Lelli et al. 2022).We fixed the kinematic centers of both galaxies to the coordinates of 1.2 mm dust continuum peak positions.These coordinates were uniform such that the rings remained concentric.The rotational velocity, velocity dispersion, systemic velocity, position angle, and inclination angle were left as free parameters.We adopted a twostage fitting strategy such that the algorithm will perform a first fitting stage to estimate all free parameters.Then, after regularizing the free parameters, the algorithm proceeds to a second stage where the radial profiles of inclination and position angles follow a Bezier function and kinematical parameters (V rot and σ V ) are allowed for free parameters (Di Teodoro & Fraternali 2015).
In the upper row of Fig. 4, we show the velocity field (mom1) and velocity dispersion (mom2) maps directly generated from the observed CO(6-5) data cubes by excluding pixels under 3σ.The noise level is calculated by averaging the standard deviation of each channel of the data cube.The position-velocity diagrams (PVDs) are shown in the lower row of Fig. 4 for NW and SE galaxies, respectively.The PVDs are extracted along the kinematic major and minor axes indicated by arrows in the mom1 and mom2 maps with a slit width the same size as the FWHM of the synthesized beam along the major axis, which are automatically produced by running 3D Barolo.The width of the region used for the extraction is the same as the region shown in Fig. D1 and D2.

SE Galaxy
The PVDs of the SE galaxy clearly show ordered rotation along the major axis.A clear velocity gradient, as well as the zero velocity region, can be seen in the main structure in the mom1 map generated from Cycle 6 observations (see Fig. B1), which reveals a PA of 132 • .This value is consistent with PA = 130 • retrieved by 3D Barolo, suggesting that the rotating structure in the PVD corresponds to the structure in #Main.This is also supported by that the main structure has a bimodal line profile characteristic of rotation (Fig. 2).In addition to these rotational molecular gas, the PVD along the minor axis shows non-circular motions with offset > 0.05 ′′ .The component with V LOS ∼ −100 km s −1 corresponds to the #Tidal but is partially smeared because of the larger beam size of C46.Because of the tidal effect, which may give rise to the tidal bridge between the two galaxies (see §7.1), the filamentary structure in #Tidal cannot be simply treated as a part of rotation subject to the gravitational force.
The modeling returns an inclination angle of 28 • .This value varies by only 1 • from the innermost to outermost rings despite the gas being distorted by the tidal force.With the inclination angle, we can calculate the ratio of the rotational velocity (V rot = V LOS / sin i) to velocity dispersion for a series of concentric rings returned by 3D Barolo.These rings have a minimum ratio of V rot /σ V ≥ 5, indicative of a rotation-dominated disc (Simons et al. 2017).

NW Galaxy
The high-resolution observations reveal non-circular and tangential motions that are diluted because of the beam-smearing effect in Cycle 2. These non-circular motions include the massive (M H 2 ∼ 2 × 10 9 M ⊙ ), lopsided molecular outflow (see §7.3), which corresponds to NW − 1 in Fig. 1.This outflow component is clearly reflected in the NW galaxy PVDs with V LOS ≤ −300 km s −1 (Fig. 4).The diffuse molecular gas in the northeast part (see C46 moment maps in Fig. 4 and Off-2 in Fig. 2 of Lebowitz et al. 2023) is another origin of non-circular motions and is labeled in the PVD along the minor axis.This additional CO(6-5) line emission has −100 ≲ V LOS ≲ 0 km s −1 with an offset < −0.05 ′′ as clearly shown in the C46 PVDs along the major axis.In the current stage, lack of observations primarily tracing cold molecular gas, we cannot confirm whether this feature originates from the spiral arm, thus belonging to part of the rotation, or stems from tidal effects.
In the C46 mom1 map, as well as indicated by the green curve in the observational data constructed by 3D Barolo ("VELOCITY" in Fig. D1), the projected zero-velocity region is curved.This curvature feature has been reproduced in the simulated rotating disc plus jetdriven outflows, showing that the kinematic centers could be curved rather than being a straight line (Meenakshi et al. 2022b).Consequently, the kinematic major axis and inclination angle of the CO disc cannot be unambiguously defined, as we will discuss below.
The non-circular motions make it more difficult to model the rotational velocity field.Depending on the initial guesses of position and inclination angles, the modeled PA ranges within [69 • , 82 • ], and the inclination angle ranges within [27 • , 55 • ].By exploring the PA−inclination parameter space using the SPACEPAR task of 3D Barolo, we found that the PA and inclination angles have degenerated.The main kinematic parameter influenced by this loosely constrained parameter space is the rotation velocity V rot as it is calculated by V rot = V LOS / sin i, where V LOS is the observed velocity.We then directly estimated the inclination angle using (Courteau et al. 2014) where a is the semi-major axis, b is the semi-minor axis, and q 0 is the axial ratio of a galaxy viewed edge-on.Adopting q 0 = 0.13 for late-type galaxies (Hall et al. 2012), the inclination is 24 • based on Cycle 6 observations where the NW galaxy is spatially resolved.
Hence, although we cannot determine the precise geometry of the NW galaxy, an i ∼ 27 • could be used to set an upper limit for V rot , which is ∼ 520 km s −1 .
The modeled rings have a minimum V rot /σ V of ∼ 3 and a maximum of ∼ 50, suggesting that the NW galaxy is rotationdominated as well.Although this ratio is highly varied because of the large uncertainties arising from the non-circular motions, in lowresolution Cycle 2 observations where the most non-circular motions are diluted, the velocity field map of the NW galaxy shows distinct blueshifted and redshifted regions, indicating a rotating disc (Emonts et al. 2015b;Lebowitz et al. 2023).

Dynamical Mass
As argued by Emonts et al. (2015b), the two galaxies should have low inclinations (10 − 20 degrees) such that the dynamical mass is comparable with the high stellar mass M ⋆ ∼ 5.8×10 11 M ⊙ computed by De Breuck et al. (2010).However, based on our SED fitting results, the stellar mass is M ⋆ ∼ (0.9 − 1.1) × 10 11 M ⊙ , which differs from the literature value by a factor of 5 but should be a more robust estimate (see §4.3).In this case, neither the ALMA Cycle 2 observation has significantly missed the extent of rotating discs nor an extremely low inclination is required for the dynamical mass to match the stellar mass.The inclination angle estimated from modeling gas kinematics using 3D Barolo is 28 • for the SE galaxy and 27 • , which is likely to be the lower limit (see §5.3), for the NW galaxy.We then estimate the dynamical mass of each galaxy by (Courteau et al. 2014) where R is the radius along the major axis, V obs is the observed rotation velocity, and i is the inclination angle.R is determined by fitting a 2D Gaussian component to the mom0 map of Cycle 6 observations and being deconvolved from the beam size.For the NW galaxy with R = 0.8 kpc and i = 27 • ± 2 • , the observed rotation velocity is V obs = 236 +23 −32 km s −1 , which gives a dynamical mass of M dyn,NW = 5.0 +1.1 −1.2 × 10 10 M ⊙ .And for the SE galaxy with R = 0.8 kpc, V obs = 155 +8 −10 km s −1 , and i = 28 • ± 2 • , the dynamical mass is M dyn,SE = 2.0 +0.3 −0.3 × 10 10 M ⊙ .The dynamical mass can be compared with other mass estimates (e.g., Tripodi et al. 2022).Without considering the BH mass, the NW galaxy has a stellar mass of M ⋆,NW = M dyn,NW − M H 2 ,NW ∼ 3×10 10 M ⊙ within a 0.8 kpc radius.This rough estimate likely rejects a large inclination angle (i ≳ 40 • ).Otherwise, the dynamical mass will be even smaller than the molecular gas mass.For the SE galaxy, M ⋆,NW ∼ M dyn,SE − M H 2 ,SE ∼ 1 × 10 10 M ⊙ within a 0.8 kpc radius.Within a 0.8 kpc radius, M ⋆,NW+SE is expected to be ∼ 4 × 10 10 M ⊙ , significantly smaller than the value 9.2 × 10 10 M ⊙ returned by the SED fitting.This suggests the existence of extended stellar components out of the CO discs, which can be seen from Fig. C1 that the stellar continuum extends well beyond the CO(6-5).A size comparison between stellar and molecular gas components in local LIRGs has been done by Bellocchi et al. (2022), finding out that the stellar size traced by ionized gas is typically three times that of the molecular gas.

BH mass
One crucial quantity in understanding AGN evolution is the BH mass.To provide an estimate of M BH related to the RLAGN hosted by the NW galaxy, we start with the scaling relation between central SMBH masses and the host galaxy stellar masses (Kormendy & Ho 2013).The SED fitting returns M ⋆ = (9.2± 0.4) × 10 10 M ⊙ including NW, SE, and the possible companion galaxy.Since there are no spatially resolved data allowing for an estimate of the individual galaxy, we then set an upper limit of M ⋆,NW ∼ 8 × 10 10 M ⊙ for the NW galaxy such that it does not exceed the total stellar mass in this system after considering M ⋆,SE ∼ 1 × 10 10 M ⊙ for the SE galaxy within the CO disc.The lower limit is 3 × 10 10 M ⊙ , which matches M ⋆,NW ∼ M dyn,NW − M H 2 ,NW estimated from the dynamical mass of ∼ 5 × 10 10 M ⊙ (see 5.4 for details).We note that these assumptions do not lead to significant quantitative changes in the estimates of the Eddington ratio to be discussed in §7.2.1 that are directly related to the BH mass and all successive discussions because the scaling relations are intrinsically scattered.

Jet kinetic power
The total kinetic power of AGN jets can be estimated from the flux densities of the radio hotspots.We have confirmed the existence of the SE and NW hotspots in Paper I, and the NW hotspot is situated towards the northwest of the radio core, along the jet axis as shown in Fig. 1 (see also Fig. 1 in Zhong et al. 2023).The observed flux density of the SE hotspot is larger than that of the NW one by at least a factor of 6, which is argued to be a result of Doppler boosting and/or a jet-ISM interaction in the SE galaxy (Zhong et al. 2023).Therefore, a restoration of the intrinsic flux density is required to estimate L jet .In a Doppler boosted case, the intrinsic (S int ) and observed (S obs ) are related by δ 2−α < S obs S int < δ 3−α , where α is the spectral index defined as S ν ∝ ν α , δ ≡ [γ(1 − β cos θ)] −1 is the Doppler factor, and γ ≡ (1 − β 2 ) −1/2 is the Lorentz factor.Following the scenario proposed in Paper I, that is, the SE hotspot has its flux density enhanced whereas the NW one dimmed, the upper limit of intrinsic flux density is then the case the NW hotspot is dimmed by a factor of four, i.e., δ = 4. On the other hand, the imbalanced flux densities may be only attributed to that NW and SE jets interact with the medium of different electron densities.In this case, we can set a lower limit for the intrinsic flux density, that is δ = 1.Assuming α = −1.2 between 1.4 and 4.7 GHz, the intrinsic flux density for a single jet at 1.4 GHz is then S int,1.4= 123 mJy for δ = 4 and 31 mJy for δ = 1, respectively.Adopting a relation between L jet and the radio power at 1.4 GHz (P 1.4GHz = 4πD 2 L (1 + z) −(α+1) S int,1.4ν 1.4GHz ; Cavagnolo et al. 2010): L jet ≈ 5.8 × 10 43 (P 1.4GHz /10 40 ) 0.7 erg s −1 , (5) the total jet kinetic power then ranges from 2 × 10 46 erg s −1 for the Doppler boosted case and 4 × 10 46 erg s −1 for the jet-medium interaction case, respectively.

DISCUSSION
7.1 The Dragonfly Galaxy as a Merger

A Late-stage Merger
Low-resolution ALMA Cycle 2 observations have revealed a bridgelike structure that connects NW and SE galaxies and was interpreted as tidal debris that arises from the gravitational interaction between the two rotational gaseous discs (Emonts et al. 2015b).This speculated structure has been further investigated by Lebowitz et al. (2023), arguing that this structure is the tidal bridge with M H 2 = (3 ± 1) × 10 9 M ⊙ , which is often observed in major mergers in both simulations and observations (Scoville et al. 2017;Sparre et al. 2022).Combining our kinematic modelings (see §5), these results support the consensus that the Dragonfly galaxy is a major merger constituting two likely rotating discs (Emonts et al. 2015b).A merger will undergo a pair phase between the first and the second pericentric passage and a merging phase after the second pericentric passage and before a full coalescence (McElroy et al. 2022).During the passage, gravitational torques will exert on the gas, (Mihos & Hernquist 1996;Barnes & Hernquist 1996), leading the gas material to lose angular momentum and inflow towards the centers of the merging galaxies, fueling star formation activities.The tidal bridge seen in Cycle 2 and the tails revealed by HST/NICMOS F160W image (Pentericci et al. 2001;Emonts et al. 2015b) serve as signatures by which the Dragonfly galaxy at least can be classified as a major merger at stage 2 characterized by obvious tidal bridges and tails (Larson et al. 2016).They are evident that the Dragonfly galaxy has already undergone the first pericentric passage.Additionally, recent studies (e.g, Stemo et al. 2021) of dual and offset AGNs linked to major mergers have found significantly enhanced AGN activities at bulge separations ranging from 14 − 11 kpc, attributed to the first pericentric passage, and 4 − 2 kpc, attributed the second passage.Considering also the small nuclear separation of ∼ 4 kpc and high IR luminosity (L IR ∼ 2 × 10 13 L ⊙ ), the Dragonfly galaxy could be a late-stage merger on its way to coalescence.
The starburst activities are thus an outcome of the second pericentric passage.This is supported by our results of SED fitting adopting an SFH with an additional recent burst, which shows that there is no significant difference between the instantaneous SFR and SFR 10Myr averaged over 10 Myr, but a comparably low SFR 100Myr averaged over 100 Myr for all three fittings (see Table. 4).In this scenario, the intense gas inflow accompanied by this passage results in a starburst that contributes to at least thirty per cent of the total stellar mass, which is reflected in the recent burst fraction ≳ 0.3 in Table 4.

Beads-on-a-string?
In late-stage mergers that truly have a final coalescence, overlapped gaseous and stellar discs are commonly observed (Di Matteo et al. 2007).Such overlapping of the stellar components between SE and NW galaxies traced by rest-frame UV and optical photometry can be found in HST/WFPC2 F814W and NICMOS F160W images (Pentericci et al. 2001;Emonts et al. 2015b).
From the aspect of gaseous discs, both SE and NW galaxies have their molecular contents concentrated within a sub-kpc scale region, showing no recognizable overlap other than the tidal bridge (Lebowitz et al. 2023).In Cycle 6 observations, the SE galaxy has a significant fraction (≳ 30 per cent) of the molecular gas originating from the filamentary structure in #Tidal, which is constituted of several clumps, as shown in Fig. 2. In the image of the dust thermal continuum, the NW galaxy has a diffuse and extended structure elongated towards the east.The potential tidal bridge, argued by Emonts et al. (2015b) and Lebowitz et al. (2023), links these two structures in SE and NW galaxies.Based on these features, we further discuss the Dragonfly galaxy as a high-redshift analog to a local late-stage merger, the famous Antennae galaxy (NGC 4038/4039) that has a separation of 6.6 kpc and shows a beads-on-a-string morphology of the molecular gas distribution (Elmegreen & Elmegreen 1983;Whitmore et al. 2014).
In a beads-on-a-string scenario, the filamentary structure forms from the tidal effects ascribed to the recent pericentric passage, as suggested for the Antennae galaxy (Espada et al. 2012).Each clump in #Tidal can be treated as a molecular filament, which is the string, that extends over ∼ 400 pc in the projected plane.Each filament may embed two or more supergiant molecular clouds (SGMCs) at a 100 pc scale, namely the beads, that host star clusters.Such beads and strings are reproduced in the high-resolution hydrodynamic simulations of a major merger of disc galaxies with a total mass of 10 8 M ⊙ (Teyssier et al. 2010).They are argued to be the result of the gas turbulent motions and will gain more masses through interactions between the merging pairs.Shocks originating from galaxy-galaxy collisions/interactions can result in higher turbulence such that the observed line width is determined by the turbulent motions and kinetic temperature (T kin ) of gas.In this case, we can estimate the Mach number -a measure of the velocity dispersion of the molecular cloud -following Leroy et al. (2016): where T 25K is the kinetic temperature of molecular gas divided by 25 K and σ is equal to FWHM/2.35.Assuming a typical T kin = 45 K for the warm molecular gas in high-redshift star-forming galaxies (Birkin et al. 2021), the Mach number reaches ∼ 400 with FWHM ∼ 280 km s −1 for the line profile extracted from #Tidal.This high Mach number is in line with the scenario in which the collisions between cold clouds could be supersonic (Struck 1999).Such a large value is more likely to be a result of the unresolved large-scale structure.A more reasonable scenario is that, under a 'beads-on-a-string' speculation, the observed FWHM ∼ 280 km s −1 can be a composition of several SGMCs with narrower FWHMs and different velocity centers (e.g., Whitmore et al. 2014).Assuming an FWHM ∼ 60 km s −1 which is a typical value for the SGMCs in the overlap regions of the Antennae galaxy and T kin = 45 K, the resultant Mach number is ∼ 80, in agreement with the values found in molecular clouds at 60 pc scale in mergers and SBGs (Leroy et al. 2016).This scenario requires future observations with similar resolutions and higher sensitivities to confirm.

A growing SMBH
Apart from starburst activities, the gas inflow accompanied by merger events can also fuel the growth of SMBHs (e.g., Lin et al. 2023).Therefore, we wonder whether this RLAGN can be associated with an active BH in a rapid growth phase.We, therefore, need to derive the corresponding Eddington ratio.For the reason discussed in §4, we use the SED fitting results adopting 44 GHz radio data as a conservative estimate.The RLAGN host galaxy refers to the NW galaxy throughout the discussions in this section.
AGNs can be classified as radiatively efficient (high accretion rate) and inefficient (often low accretion rate) populations based on whether its Eddington ratio is above λ Edd = 3 × 10 −2 or not, respectively (Merloni & Heinz 2008).The Eddington ratio is defined by λ Edd = L AGN,bol /L Edd , where L Edd = 1.26 × 10 38 M BH M⊙ .Since there exist no obvious observational features of AGN activities in the SE galaxy, such as high-excited molecular gas, molecular outflows, and PSF-dominated morphology in HST images (e.g., Zhong et al. 2022), the RLAGN that resides in the NW galaxy is treated as the primary contributor to the computed L AGN,bol .The Eddington ratios are then λ Edd ∼ 0.07−4 adopting M BH ∼ (0.2−10)×10 8 M ⊙ based on the BH mass estimates presented in §6.1 These values suggest that the RLAGN is likely to be radiatively efficient and the central SMBH is in a fast growth phase and may even enter the super-Eddington regime.
We also note that this λ Edd is the most conservative estimate and λ Edd could be larger by a factor of ∼ 3 if we adopt the fitting including 1.4 GHz data or without radio data.It should also be noted that the scaling relation is broadly scattered for all populations of AGNs, precipitating inevitable uncertainties in λ Edd .If our target RLAGN is a radiatively inefficient one (λ Edd < 0.03), the SMBH should at least have M BH (= L AGN,bol 1.26×(λ Edd =0.03)×10 38 ) ≥ 7 × 10 9 M ⊙ even adopting L AGN,bol = 0.9 × 10 46 erg s −1 as the lower limit.However, Poitevineau et al. (2023) has studied the RLAGN host galaxy properties at 0.3 < z < 4 and found statistical consistency in the scaling relations between RLAGNs and other AGNs populations.This colossal mass will make the NW galaxy an outlier extremely deviating from the scaling relations, and only a few such objects are discovered.We, therefore, forward our discussions without considering this uttermost possibility and we rely on the radiatively efficient scenario for the RLAGN in the NW galaxy.

Dragonfly as an Extremely Radio-loud Galaxy
Defining the ratio of L jet to L AGN,bol as the intrinsic radio loudness R int , we find that log R int is likely to lie above 0, following L jet calculated in §6.2 and L AGN,bol ∼ 0.9 × 10 46 erg s −1 listed in Table .4. We further calculate the specific black hole accretion rate, which is defined as sBHAR The obtained values of R int and sBHAR are similar to the average values (⟨log R int ⟩ = 0.0 and ⟨log (sBHAR/ erg s −1 M ⊙ −1 ⟩ = 35.3) of extremely radioloud galaxies (ERGs) defined by log( f 1.4GHz,rest / f g band,rest ) > 4 with log L 1.4GHz > 10 24 W Hz −1 (Ichikawa et al. 2021).Therefore, the RLAGN in the NW galaxy is classified as an ERG.
At low redshifts, it has been well established that the powerful radio galaxies with jets are representative of the final evolution stage of AGNs (Hopkins et al. 2008).These AGNs host SMBHs surrounded by accretion discs with weak gas inflow, so-called advection-dominated accretion flow (ADAF), which leads to low bolometric luminosity but hard spectrum.This accretion state is called as the low/hard state and the accretion disc is in the sub-Eddington phase.Based on the findings presented here, however, it appears that the Dragonfly galaxy, including ERGs, represents a distinct population, which involves AGNs with rapidly growing SMBHs at high redshifts.

Origin of jets from a radiatively efficient SMBH
The estimates for λ Edd in 7.2.1 merely represent a statistical expectation since it is directly related to M BH .However, the intrinsic dispersion around the M BH − M ⋆ scaling relation may still push the SMBH residing in the NW galaxy (λ Edd ∼ 0.07 − 4) into the super-Eddington accretion regime (λ Edd ≳ 1).In this case, the jets in the Dragonfly galaxy, as well as ERGs, might be launched from a super-Eddington accretion disc (or slim disc; Abramowicz et al. 1988) being both optically and geometrically thick.The disc thickness prevents the diffusion of magnetic fields, and these jets are powered by large-scale magnetic fields in the innermost region and rotating SMBHs and can have an Eddington-order luminosity.Simulations of jet-disc connections show that jets can be persistently driven by the Blandford-Znajek mechanism (BZ mechanism; Blandford & Znajek 1977) when the accretion flow is super-Eddington (McKinney et al. 2014;Qian et al. 2018).
The super-Eddington accretion is short-lived and episodic with a timescale ranging 10 4 − 10 7 yr, depending on the replenishment of cold gas from host galaxies (Volonteri et al. 2015).This is because the gas accretion timescale is inversely proportional to the mass accretion rate squared (Kato et al. 2008).The accretion rate of a super-Eddington thick disc launching powerful jets will decrease as time evolves and finally enter the standard 'thin' disc regime.Shown by the simulation performed by Ricarte et al. (2023), following this state transition, the jet power shows a rapid decrease by more than three orders of magnitude.The jets are then incapable of maintaining the steady state of the synchrotron radiation spectrum.As a result, the jets appear to be switched off and no longer interact with the medium, which is speculated in Paper I (Zhong et al. 2023).
From this aspect, RLAGNs showing high accretion rates at high redshifts such as the Dragonfly galaxy may represent an indispensable stage of the galaxy-BH co-evolution, through which BHs rapidly gain their masses.At the same time, the intense star-forming activities lead to a rapid depletion of the molecular gas in the merger as the merging pairs gradually coalesce.In the end, the accretion rate drops down and the accretion disc is finally settled in the low/hard state.Correspondingly, the slowly accreting SMBH then resides in a massive, low-SFR ETG with giant radio lobes, which is the typical picture we observe in the nearby Universe.
Another possibility is established upon RLAGNs being analogs to black hole X-ray binaries (BHBs) and microquasars that launch transient jets.It is widely argued that jets are likely to be launched, which are called continuous jets when the accretion disc is geometrically thick (e.g., Narayan et al. 1998;Meier 2001;S ądowski et al. 2014;S ądowski & Narayan 2015;Tchekhovskoy 2015;Blandford et al. 2019), while geometrically thin discs have low continuous jetlaunching probabilities since the magnetic fields will diffuse outward at a rate faster than being dragged inward, incapable of powering the jets (Lubow et al. 1994;Guilet & Ogilvie 2012).However, by modeling jet-disc coupling in BHB systems, Fender et al. (2004) found that transient jets are launched as the accretion disc transitions from low/hard to high/soft state.In the low/hard state, the system launches weak continuous jets.Then the disc luminosity increases accompanied by the emergence of powerful transient jets, and the disc is in an intermediate state.When the accretion disc is finally settled in the high/soft state, there will be no jets (Fender et al. 2004).The power of transient jets is argued to correlate with the BH spin, but how come they are correlated remains a mystery.Observational evidence of such a transitioning phase is limited to BHBs or microquasars for which the state transition timescale is short enough.State transitions in AGN or quasars, if scaled by the orders of magnitude of central BH masses against BHBs, have a period of 10 4 − 10 7 yr (Tchekhovskoy 2015).As studied in Paper I, the RLAGN may have transient or intermittent activities (Zhong et al. 2023).The synchrotron age, i.e., the lifetime of radio jets, has an order of ∼ 10 5 yr.Therefore, the galaxy merging results in the fueling of the BH and we are likely to be observing an RLAGN whose accretion has transitioned from a low/hard state to a high/soft one.

Outflow geometry
To better investigate the high-velocity component (NW−1 in Fig. 1), we have generated the mom0 map of the NW galaxy integrated over −500 to −280 km s −1 , which is shown in Fig. 5.This makes it clear that the high-velocity component is likely to extend along a direction perpendicular to the jet axis.This component also shows lopsidedness -only towards the southwest direction.A similar lopsided outflow cone has been observed in the radio galaxy B2 0258+35 and reproduced in simulations (Murthy et al. 2022).In B2 0258+35, the non-detection of the outflow on the opposite side of the monopolar outflow may be attributed to the obscuration by the CO disc or a much larger clumpiness of the ISM.This is likely to be the case for the NW galaxy, as evidenced by the diffuse dust continuum -which could be a result of galaxy interactions (see §7.1.2) -that extends towards the east direction (panel (c) in Fig. 1).
What is even more interesting revealed by Cycle 6 observations is that, along the jet axis and towards the northwest, there exists weak but extended CO(6-5) line emission.Lebowitz et al. (2023) also found a broad, high-velocity component at this position based on combined Cycle 2 and 4 observations tapered to a low resolution of 0.19 ′′ .Additionally, at this position, Lebowitz et al. (2023) detected a faint blob in the 237 GHz dust continuum (see their Fig.8), which appears to be an extended, diffuse emission just above the centroid of the NW galaxy in our Cycle 4 dust continuum image (see panel (a) in Fig. 1).If these observational features in Cycle 4 and 6 observations are real, then they may be indicative of an interaction between the NW jet and CO disc in a few tens pc-scale, which is proposed to explain the high-velocity component coincided with the faint blob by Lebowitz et al. (2023).Because of this jet-ISM interaction, the propagation direction of the NW jet is then distorted, resulting in a further misalignment in the inclinations of approaching and receding jets relative to the LOS, as discussed in the Doppler boosting effect in Paper I (Zhong et al. 2023).Numerical simulations show that an inclined jet could strongly couple with the ISM, leading to sub-relativistic and wide-angled outflows (Mukherjee et al. 2018;Talbot et al. 2022).Such a jet-ISM interaction in the innermost region is proposed for the geometry of the narrow-line region of NGC 1068, where splitted CO line emission is observed from the torus to the edge of the circumnuclear disk (García-Burillo et al. 2019).Due to the low signal-to-noise ratio (SNR = 4) of the detection in the current Cycle 6 dataset, this circumnuclear disc-scale jet-ISM interaction is more like a hypothesis and requires future observations with higher sensitivities to confirm.

Outflow energetics
We further investigate the outflow properties and energetics by adopting the time-averaged, thin-shell approximation (Rupke et al. 2005) as a conservative estimate instead of a scenario in which the outflow occupies a spherical volume with a uniform density and has a mass outflow rate larger than the thin-shell scenario by a factor of 3 (Maiolino et al. 2012;González-Alfonso et al. 2017).The mass outflow rate is given by where v out ∼ 450 km s −1 is the maximum outflow velocity defined as ∆v broad + FWHM/2 and ∆v broad is the line center shift between the narrow and broad components.We determine the outflow radius by fitting a Gaussian component to the mom0 map (Fig. 5) integrated from −500 to −280 km s −1 , by which R out ∼ 0.5 kpc based on Cycle 6 observations.The corresponding lower and upper limit, given the molecular gas mass M H 2 ∼ (1.3 − 2.8) × 10 9 M ⊙ (NW − 1 in Table 2), for the mass outflow rate is 1200 +300 −300 ≲ Ṁout ≲ 2600 +600 −600 M ⊙ yr −1 .The associated outflow kinetic power is Ėout = 1 2 Ṁout × v 2 out ∼ (8 − 16) × 10 43 erg s −1 and momentum boost factor ranges within 11 − 24 for L AGN,bol ∼ 0.9 × 10 46 erg s −1 and 3 − 7 for L AGN,bol ∼ 3 × 10 46 erg s −1 , where ṖAGN is the AGN radiation momentum rate.
We compare these outflow properties with the values in literature in Fig. 6, including a sample of low-redshift AGNs (Fluetsch et al. 2019), four radio-loud quasar host galaxies at z = 1.4 − 2.3 (Vayner et al. 2021b), a Type 2 QSO at z = 0.085 (Audibert et al. 2023), and AGNs at z ∼ 0.1 − 3 (Fiore et al. 2017).The Dragonfly galaxy shows a high molecular outflow rate comparable with the values found in L AGN,bol -matched AGN host galaxies for both molecular and ionized outflows.This high rate is a combined result of the compactness (∼ 0.5 kpc) and massiveness (∼ 2×10 9 M ⊙ ).However, this mass represents merely ≲ 10% of the total molecular gas in the NW galaxy.Additionally, the starburst activity of the Dragonfly galaxy consumes the molecular gas reservoir at a rate (SFR ∼ 2000 − 3000 M ⊙ yr −1 ) comparable with Ṁout .These results suggest that, although the AGN leads to instantaneous negative feedback in the Dragonfly galaxy, it is the long-term and combining effects of AGN and star-forming activities that are responsible for the final quenching, in line with the ideas proposed in recent simulations (e.g., Bollati et al. 2023).

Radiative-mode
Based on the estimates of outflow energetics, we first discuss the situation that the outflow stems from the radiative-mode feedback.In this scenario, the black hole wind drives the expansion of a hot shocked wind bubble that interacts with the host galaxy ISM surrounding the central BH (King 2003;King & Pounds 2015).If this energy bubble suffers efficient radiative cooling and only transmits the ram pressure to the ISM, the outflow is momentum-driven.And if the bubble expands adiabatically, the outflow is energy-driven.
These results yield previous studies of multiphase outflows in L AGN,bol -matched AGNs (e.g., Fiore et al. 2017;Fluetsch et al. 2019) and QSOs at different redshifts (e.g., Bischetti et al. 2019;Vayner et al. 2021a), showing that AGN radiative activities are sufficient to explain the outflows observed in systems.They are also in line with the molecular outflows found in L jet -and ∼L AGN,bol -matched radio-loud quasars at z ∼ 2 (Vayner et al. 2021b) that reveal energyconserving outflows.

Jet-mode
However, as shown in §6.2, L jet is comparable and can even exceed L AGN,bol , thus we have to consider a second situation that the outflow stems from the jet-mode feedback.The radio jets can strongly couple with its host galaxy ISM if the jet has a small angle relative to the disc plane of the host galaxy (Wagner & Bicknell 2011;Wagner et al. 2012;Mukherjee et al. 2016).The dense molecular clouds can be dispersed away from the jet axis through the ram pressure, forming an expanding cocoon of the molecular gas along the jet path.The jet flow will channel in various directions because of the porosity of the ISM, resulting in a larger-scale impact on the host galaxy.However, the outflow component in the NW galaxy extends no more than 0.5 kpc in the projection.Along with the enhanced velocity dispersion perpendicular to the jet axis, this compactness, as indicated by the simulation performed by Meenakshi et al. (2022b), may represent the early phase of the jet-driven outflow.
In many previous studies of outflows in nearby jetted systems with radio lobes, the structures of multiphase outflows traced by integrated intensity maps are often found to be elongated along the jet axis (Oosterloo et al. 2017;Girdhar et al. 2024).On the other hand, the [O III]λ5007 line velocity widths4 , which is indicative of enhanced velocity dispersion, are found to extend along the direction perpendicular to the jet axis (see Venturi et al. 2021;Meenakshi et al. 2022b and references therein).In the case of our target, the outflowing molecular gas is perpendicular to the jet axis in both integrated intensity and velocity dispersion maps.A similar observational example is the local type 2 QSO SDSS J1430 (L AGN,bol ∼ (0.7−6)×10 45 erg s −1 ; Venturi et al. 2023) whose radio jets have kinetic power of L jet ∼ (1 − 3) × 10 43 erg s −1 .The molecular gas outflow is argued to emanate from the lateral outflow/turbulence constituted of warm dense gas excited by the cocoon of shocked gas (Audibert et al. 2023).
The jet-medium coupling efficiency is only η = Ėout /L jet ∼ 0.03 for the Dragonfly galaxy.Considering the inevitable large uncertainties in the empirical L jet − P 1.4GHz relation, this value represents an upper limit for η.The 1.4 GHz flux density is corrected for the particle acceleration and a direct estimate of the jet power could have an order of magnitude increase (≳ 2 × 10 47 erg s −1 ), resulting in an extremely low coupling efficiency of η ∼ 0.001.These results suggest that, if the molecular outflow is the dominant phase, the radio jets are weakly coupled with the ISM.None the less, adopting η = 0.01 and L jet = 2 × 10 46 erg s −1 , the jets provide an total energy injection of ∼ 2 × 10 57 erg during its active timescale of ∼ 2 × 10 5 yr, comparable with the outflow kinetic energy of E out = 1 2 M out × v 2 out ∼ 3 × 10 57 erg.Albeit this is no more than an order of magnitude estimate, it further shows that the jets are capable of expelling the cold gas during its short active phase.Based on the recurrent jet scenario proposed in Zhong et al. (2023), the jets may finally remove all molecular gas from its host galaxy through long-term, episodic bursts.
On the other hand, both the outflow mass and kinetic power could be governed by the ionized rather than the molecular phase, which has been observed in RL quasars (Vayner et al. 2021b).For the target with log(L AGN,bol ) ∼ 47 erg s −1 shown in Fig. 6, the ionized phase has a mass outflow rate larger than that of the molecular phase by almost an order of magnitude and requires more energy to be powered.Similarly, in HzRGs with both L jet and L AGN,bol ≳ 1 × 10 47 erg s −1 at z ∼ 2 (Nesvadba et al. 2017), Ėout of the ionized outflow can be two orders of magnitude larger than those derived for the Dragonfly galaxy.In addition, as studied by Venturi et al. (2023), although SDSS J1430 is weak in molecular outflow, its ionized outflow may require a hundred per cent jet-ISM coupling efficiency in a jet-mode scenario, while the ionized outflow kinetic power and momentum boost factor agree well with the radiative-mode scenario.Actually, the simulation shows that when L jet ∼ L AGN,bol , jets can more efficiently shock-ionize gas than AGN radiation even jets are aligned 70 • relative to the galactic plane (Meenakshi et al. 2022a).Therefore, for the Dragonfly where L jet ≳ L AGN,bol , observations of ionized gas are then essential to investigate whether and how -in outflow mass and/or kinetic power -ionized outflow dominates over the molecular phase.
All in all, the jet-mode feedback is capable of driving massive molecular outflow within a short timescale without additional energy supply from AGN radiation.On the other hand, the radiative activity of this RLAGN is adequately responsible for the outflow without the assistance of powerful jets.This makes HzRGs such as the Dragonfly galaxy a unique AGN population amongst which jet-(characterized by the existence of jets) and radiative-mode (characterized by high L AGN,bol ) feedback should co-exist while their relative importance in shaping host galaxies remain unsettled.Additionally, unlike the RL quasars studied by Vayner et al. 2021b showing consistent distributions between jet propagation axes and outflows, the fact that the outflow in the NW galaxy is perpendicular to the jet axis differs from the common picture.Future statistical studies of multiphase outflows in radiatively efficient RLAGNs are required to explore the connections between L jet , L AGN,bol , jet ages, and outflow kinematics.

CONCLUSION
In this work, we have studied the molecular gas traced by CO(6-5) line emission in the Dragonfly galaxy, a hyper-luminous infrared, starburst major merger at z = 1.92 using ALMA and VLA observations.We have studied the gas kinematics using 3D Barolo and performed SED fitting covering optical-to-radio using CIGALE.Our major findings are as follows: (i) The SED fitting using optical-to-44 GHz photometric data gives an SFR ∼ 3100 M ⊙ yr −1 averaged over 10 Myr and SFR ∼ 600 M ⊙ yr −1 averaged over 100 Myr, implying a recent starburst that contributes to thirty per cent of the total stellar mass.The stellar mass of the Dragonfly galaxy computed by the new fitting is ∼ 9×10 10 M ⊙ , which is one-fifth the value computed using a SED template for early-type galaxies in previous studies (De Breuck et al. 2010).The fitting including the radio core flux density returns an AGN fraction of 0.1.If using the flux density at 1.4 GHz originating from radio hotspots and lobes in the SED fitting, the AGN fraction becomes 0.3.The corresponding AGN bolometric luminosity is L AGN,bol ∼ 2.9 × 10 46 erg s −1 , which is larger than the cases without 1.4 GHz data where L AGN,bol ∼ 0.9 × 10 46 erg s −1 but consistent with the fitting without radio data.(ii) The bulk of the molecular gas in both SE and NW galaxies is concentrated within a radius of ∼ 0.8 kpc.The NW galaxy shows a simple, circular distribution of the bulk of the molecular gas while some extended, diffuse molecular gas exists due to tidal effects.For the SE galaxy, the molecular gas is primarily constituted of two structures, named #Main and #Tidal (Fig. 2).One shows a double-peaked line profile characteristic of a rotating structure.Another one is likely to be associated with tidal effects and has its line center significantly blueshifted by ∼ 200 km s −1 In a 'beads-on-a-string' speculation, this structure may embed several supergiant molecular clouds with a Mach number of ∼ 80, in line with the values derived for cold clouds at 60 pc-scale in mergers and starburst galaxies.(iii) The gas kinematic modelings of the SE galaxy find a position angle of 130 • , consistent with the orientation of extension of the main structure named #Main.The SE galaxy has a rotation velocity to velocity dispersion ratio of at least V rot /σ V ≥ 5, indicative of being rotation-dominated.The molecular gas in the NW galaxy is highly disturbed because of the AGN activities and possible tidal effects, leading to a loosely constrained inclination angle of ∼ 28 • , in line with the value of ∼ 24 • directly derived from the integrated moment map based on the aspect ratio.The NW galaxy has V rot /σ V ≳ 3 and is therefore likely to be rotating as argued by Emonts et al. (2015b) and Lebowitz et al. (2023).(iv) Using the M BH − M ⋆ scaling relations for X-ray-selected AGNs at high-redshifts, we find M BH ∼ (0.2 − 10) × 10 8 M ⊙ , corresponding to an Eddington ratio of λ Edd ∼ 0.07 − 4. If the SMBH is accreting at the super-Eddington rate, the powerful jets might be launched from an optically-and geometricallythick super-Eddington accretion disc.Given also the young age and possible transient and intermittent activities of this RLAGN, the jets might be transient and launched from an accretion disc in a state transition, either from low/hard to high/soft or from a super-Eddington one to the standard thin disc.(v) The molecular outflow has an outflow rate of ∼ 1200 +300 −300 − 2600 +600 −600 M ⊙ yr −1 , comparable with the powerful QSO systems.If the radiative-mode AGN feedback dominates the outflow, the outflow kinetic power of ∼ (8 − 16) × 10 43 erg s −1 and momentum boost factor of 3 − 7 adopting L AGN,bol ∼ 0.9 × 10 46 erg s −1 and of ∼ 11 − 24 adopting L AGN,bol ∼ 03 × 10 46 erg s −1 suggest that the outflow could be energyconserving or be driven through direct radiation pressure on dust.If the outflow is a result of the powerful jets, i.e., the jetmode feedback, the jet-ISM coupling efficiency is ≲ 0.03 and could be as low as ≲ 0.001.During its short active timescale, the jets can provide a total energy injection of ∼ 2 × 10 57 erg with a jet-ISM coupling efficiency of 0.01, which is comparable with the outflow kinetic energy of ∼ 3 × 10 57 erg.Based on available observations, which mode dominates the outflow remains mysterious.(vi) The outflowing molecular gas in the NW galaxy is compact with a radius of ∼ 500 pc in projection and appears to be perpendicular to the jet axis and lopsided towards the blueshifted side of the NW galaxy.Although it is very massive (M H 2 ∼ (1.1 − 2.3) × 10 9 M ⊙ ), this mass is only ∼ 10% of the total cold gas deposited in the NW galaxy, suggesting that the AGN feedback is negative but not responsible for the quenching in the short-term.However, considering the possible intermittent nature of the powerful jets, the jet-mode feedback may quench the host galaxy through the long-term, episodic bursts of jets that will eventually clear out the cold gas.

Figure 2 .
Figure 2. The zoom-in investigation of the SE galaxy with white contours for CO(6-5) and black contours for 237 GHz continuum laid on the CO(6-5) mom0 map based on Cycle 6 observations.We divide the mom0 map of the SE galaxy into two regions: #Main and #Tidal.The #Tidal region is named after its possible association with the tidal bridge (see §7.1.1).There is a continuum "sub-component" ∼ 1 kpc to the west of #Main and coincided with the location of the SE hotspot.The line profiles extracted from #Main and #Tidal regions using the elliptical apertures indicated in the figure are shown in the right column.The structure in #Main shows a bimodal line profile characteristic of a rotating structure.The ellipse in the bottom left corresponds to the beam size.

Figure 3 .
Figure 3.The best fitting models of UV-to-radio SED for the Dragonfly galaxy.The upper figure shows the fitting up to 44 GHz of the radio core.The lower one is up to 1.4 GHz without the 44 GHz flux density.In each figure, the upper panel shows the observed and modeled flux densities against the observed wavelength and the lower panel shows the residuals of the model in each corresponding band.The reduced chi-square defined by CIGALE is χ 2 red,CIGALE = χ 2 /(data points − 1), where the degree of freedom is not involved.We re-calculated χ 2 red considering the number of fitted free parameters.

Figure 4 .
Figure 4. From left to right, the panels show the velocity field (mom1), the velocity dispersion (mom2) maps generated from the combined C46 dataset by excluding pixels under 3σ, and the position velocity diagrams (PVDs) extracted from the observed and modeled data cube.The upper row shows the results from the NW galaxy and the lower row for the SE galaxy.In PVDs, the color code (with black contours) indicates the observational data, the cyan contours indicate the modeled data, and the yellow dots indicate the modeled rotational velocity.∆V LOS is the observed velocity field V LOS corrected for systemic velocity (∆V LOS = ±||V LOS | − |V sys ||).The black arrows in mom1 maps indicate the slits of major and minor axes used to extract the PVD.The slits here only indicate the extraction direction and the slit width is approximately the beam size.

Figure 6 .
Figure 6.Molecular outflow properties traced by CO(6-5) of the RLAGN in the NW galaxy are shown by the five-pointed stars.Left: mass outflow rate as a function of L AGN .Middle: outflow kinetic power as a function of L AGN .Right: momentum boost factor as a function of outflow velocity, while the coral stars indicate the calculation adopting L AGN,bol ∼ 0.9 × 10 46 erg s −1 , and the blue one adopting L AGN,bol ∼ 3 × 10 46 erg s −1 .Through all panels, dots indicate the molecular outflow of low-redshift AGNs and SFGs from Fluetsch et al. (2019), filled (open) diamonds indicate the molecular (ionized) outflow from radioloud quasar host galaxies at z = 1.4 − 2.3 from Vayner et al. (2021b), triangle indicates the upper limit of the molecular outflow in the Type 2 QSO SDSS J1430+1339 at z = 0.085 (Audibert et al. 2023), where L AGN,bol ∼ (0.7 − 6) × 10 45 erg s −1 from Venturi et al. (2023) is used, and crosses indicate the ionized outflow from AGNs at z ∼ 0.1 − 3 (Fiore et al. 2017) and recomputed based on Eqs.(7-8).In the middle panel, the dashed lines indicate the thermal-tokinetic conversion efficiency of 100, 10, and 1%, while the dotted line indicates 0.06% for the upper limit of the momentum-driven molecular outflow with ϵ kin = 6 × 10 −4 vout∼1000 km s −1 1000 km s −1(Wagner et al. 2012;Costa et al. 2014).In the right panel, the dashed and dotted lines indicated the expected momentum boost factor for energy-( Ṗout / ṖAGN = 20) and momentum-driven ( Ṗout / ṖAGN = 1) outflow, respectively.

Figure B1 .
Figure B1.From left to right: the integrated intensity, velocity field, and velocity dispersion maps based on Cycle 6 observations.The upper panels show those of the NW galaxy and the lower ones are for the SE galaxy.The horizontal black line indicates a 1 kpc physical scale.The yellow dashed line indicates the jet propagation direction based on VLA 44 GHz observation and horizontally translated by 0.025 ′′ towards the east.

Figure C1 .
Figure C1.From left to right, the panels show the HST/NICMOS F160W imaging that traces the rest-frame optical continuum overlaid with ALMA Cycle 4 continuum, Cycle 6 continuum, and Cycle 4+6 CO(6-5) line emission.The grey crosses indicate the centroid of the NW and SE galaxies.

Figure D1 .
Figure D1.The kinematic modeling results using 3D Barolo based on C46 observation for the NW galaxy.From upper to bottom, the panels are integrated intensity, velocity field, and velocity dispersion maps.From left to right, the panels show the maps directly built from the observational data cubes, the maps built from the modeled data, and the residuals after subtracting the modeled from the observed data.The dashed line indicates the major axis, the grey points indicate the fitted rings, and the ellipse indicates the outermost ring.

Table 1 .
Summary of the Observations

Table 2 .
Measured and Derived Physical Properties from the Line Profiles

Table 4 .
Host galaxy and AGN physical properties computed using CIGALE