Molecular Gas Kinematics and Star Formation Properties of the Strongly-lensed Quasar Host Galaxy RXS J1131–1231

Observations of the $\mathrm{CO}(J=3\to 2)$ rotational line (${\nu }_{\mathrm{rest}}=345.79599$ GHz) toward RXJ1131 at ${z}_{{\rm{s}},\mathrm{QSO}}=0.658$ were carried out with the Combined Array for Research in Millimeter-wave Astronomy (CARMA; Program ID: cf0098; PI: D. Riechers) in the D array configuration. The line frequency is redshifted to ${\nu }_{\mathrm{obs}}=209.10443$ GHz at the quasar redshift. Observations were carried out on 2014 February 02 under poor 1.5 mm weather conditions and on 2014 February 17 under good 1.5 mm weather conditions. This resulted in a total on-source time of 2.94 hr after flagging poor visibility data. The correlator setup provides a bandwidth of 3.75 GHz in each sideband and a spectral resolution of 12.5 MHz (∼17.9 km s⁻¹). The line was placed in the lower sideband with the local oscillator tuned to ${\nu }_{\mathrm{LO}}\sim 216\,\mathrm{GHz}$. The radio quasars J1127−189 (first track) and 3C273 (second track) were observed every 15 minutes for pointing, amplitude, and phase calibration. Mars was observed as the primary absolute flux calibrator and 3C279 was observed as the bandpass calibrator for both tracks.

Given that the phase calibrator used for the first track was faint and was observed under poor weather conditions and that the phase calibrator used for the second track was far from our target source, the phase calibration is subpar, with an rms scatter ∼50° over a baseline length of ∼135 m. We thus conservatively estimate a calibration accuracy of ∼40% based on the flux scale uncertainties, the gain variations over time, and the phase scatter on the calibrated data. We therefore treat the derived $\mathrm{CO}(J=3\to 2)$ line intensity with caution and ensure that our physical interpretation of this system and the conclusion of this paper do not rely on this quantity.

The miriad package was used to calibrate the visibility data. The calibrated visibility data were imaged and deconvolved using the CLEAN algorithm with "natural" weighting. This yields a synthesized clean beam size of 32 × 19 at a position angle (PA) of 8° for the lower sideband image cube. The final rms noise is σ = 13.3 mJy beam⁻¹ over a channel width of 25 MHz. An rms noise of σ = 0.83 mJy beam⁻¹ is reached by averaging over the line-free channels in both sidebands.

Observations of the CO(J = 2 → 1) rotational line (${\nu }_{\mathrm{rest}}=230.53800\,\mathrm{GHz}$) toward RXJ1131 were carried out using the IRAM Plateau de Bure Interferometer (PdBI; Program ID: S14BX; PI: D. Riechers). Based on the CARMA $\mathrm{CO}(J=3\to 2)$ line redshift of ${z}_{\mathrm{CO}(3-2)}=0.655$, the CO(J = 2 → 1) line is redshifted to ${\nu }_{\mathrm{obs}}=139.256$ GHz. Two observing runs were carried out on 2014 December 06 and 2015 February 05 under good weather conditions in the C and D array configurations, respectively. This resulted in 3.75 hr of cumulative six antenna-equivalent on-source time after discarding unusable visibility data. The 2 mm receivers were used to cover the redshifted CO(J = 2 → 1) line and the underlying continuum emission, employing a correlator setup that provides an effective bandwidth of 3.6 GHz (dual polarization) and a native spectral resolution of 1.95 MHz (∼4.2 km s⁻¹). The nearby quasars B1127−145 and B1124−186 were observed every 22 minutes for pointing, secondary amplitude, and phase calibration, and B1055+018 was observed as the bandpass calibrator for both tracks. MWC349 and 3C279 were observed as primary flux calibrators for the C and D array observations, respectively, yielding calibration accuracy better than 15%.

The gildas package was used to calibrate and analyze the visibility data. The calibrated visibility data were imaged and deconvolved using the CLEAN algorithm with "natural" weighting. This yields a synthesized clean beam size of 444 × 195 (PA = 13°). The final rms noise is σ = 1.45 mJy beam⁻¹ over 10 MHz (21.5 km s⁻¹). The continuum image at ${\nu }_{\mathrm{cont}}\sim 139$ GHz is produced by averaging over 3.16 GHz of line-free bandwidth. This yields an rms noise of 0.082 mJy beam⁻¹.

Our analysis also uses archival data of the 4.885 GHz radio continuum obtained with the Karl G. Jansky Very Large Array (VLA) (Program ID: AW741; PI: Wucknitz). Observations were carried out on 2008 December 29 under excellent weather conditions in the A array configuration for a total of ∼7 hr of on-source time. The C-band receivers were used with a continuum mode setup, providing a bandwidth of 50 MHz for the two IF bands with full polarization. The nearby radio quasar J1130−149 was observed every 10 minutes for pointing, amplitude, and phase calibration. J1331+305 was observed as the primary flux calibrator, and J0319+415 was observed as the bandpass calibrator, yielding ∼10% calibration accuracy. We used aips to calibrate the visibility data. The calibrated visibility data were imaged and deconvolved using the CLEAN algorithm with robust = 0, which was chosen to obtain a higher resolution image given high S/N. This yields a synthesized clean beam size of 049 × 035 (PA = 018) and a final rms noise of σ = 13 μJy beam⁻¹.

We obtained a Hubble Space Telescope (HST) image taken with the ACS using the F555W filter (V-band) from the Hubble Legacy Archive. The details of the observations can be found in C06. We apply an astrometric correction to the optical image by adopting the VLA 5 GHz map as the reference coordinate frame. We shift the latter to the east by 05963 in R.A. and +08372 in decl., which is consistent with the typical astrometric precision (1''–2'') of images from the Hubble Legacy Archive.¹ This astrometric correction is critical to avoid artificial spatial offsets between different emitting regions and to carry out our lens modeling, in which the absolute position of the foreground lensing galaxy is based on its coordinates in the high-resolution optical image. The VLA image is calibrated using a well-monitored phase calibrator, with an absolute positional accuracy of ∼2 mas. For this reason, the absolute alignment between the VLA image and other interferometric images reported in this paper are expected to have an astrometric precision better than 01, modulo uncertainties related to the SNR and phase instability.

We detect CO(J = 2 → 1) line emission toward the background source RXJ1131 in the PdBI data at ≳27σ significance. Based on this measurement, we refine the redshift of RXJ1131 to ${z}_{\mathrm{CO}}=0.6541\pm 0.0002$.² The emission is spatially and kinematically resolved with a highly asymmetric double-horned line profile as shown in Figure 1. Fitting a double-Gaussian results in peak flux densities of 75.3 ± 2.6 mJy and 24.0 ± 2.0 mJy, and an FWHM of 179 ± 9 km s⁻¹ and 255 ± 28 km s⁻¹ for the two components, respectively. The peaks are separated by ${\rm{\Delta }}{v}_{\mathrm{sep}}=400\pm 12$ km s⁻¹. The total integrated line flux is 20.3 ± 0.6 Jy km s⁻¹.

Zoom In Zoom Out Reset image size

We construct a zeroth-order moment map, red/blue channel maps, and first and second moment maps, as shown in Figure 2, using the uv-continuum-subtracted data cube over a velocity range of ${\rm{\Delta }}v\sim 750$ km s⁻¹. The higher-order moment maps are produced using unbinned channel maps with 3σ clipping. The peak flux density is 8.12 ± 0.30 Jy km s⁻¹ beam⁻¹ in the intensity map. Observed properties of the CO(J = 2 → 1) emission line are summarized in Table 1.

Zoom In Zoom Out Reset image size

Table 1.  Observed Properties of RXJ1131 and Its Companion

Parameter	Unit	Value
${z}_{\mathrm{CO}(2-1)}$		0.6541 ± 0.0002
${I}_{\mathrm{CO}(2-1)}$	Jy km s−1	20.3 ± 0.6
${S}_{\mathrm{CO}(2-1)}$ a	Jy km s−1 beam−1	8.12 ± 0.30
${\mathrm{FWHM}}_{\mathrm{CO}(2-1)}$ b	km s−1	179 ± 9, 255 ± 28
${\mathrm{FWHM}}_{\mathrm{CO}(2-1)}$ c	km s−1	220 ± 72
${I}_{\mathrm{CO}(3-2)}$	Jy km s−1	35.7 ± 6.9

Notes.

^aPeak flux density in the intensity map. ^bFrom fitting a double Gaussian to the observed CO(J = 2 → 1) spectrum (Figure 1). ^cFrom fitting a double Gaussian with a common FWHM to the de-lensed CO(J = 2 → 1) spectrum (Figure 5).

Download table as:  ASCIITypeset image

The deconvolved source size FWHM obtained from fitting a single two-dimensional Gaussian to the integrated line emission in the image plane yields 51 ± 07 × 37 ± 07, which is consistent with that obtained by visibility-plane fitting within the uncertainties. Since the spatial distribution of the observed CO emission is unlikely to be fully described by a simple Gaussian and appears to be a superposition of at least two components (top left panel of Figure 2), we also fit two Gaussians to the intensity map. This yields deconvolved source sizes of 38 ± 04 × 19 ± 04 and 36 ± 03 × 15 ± 03, separated by ∼22 in R.A. and ∼17 in decl. The deconvolved source sizes of both models suggest that the gravitationally lensed CO emission is more extended than the optical "Einstein ring," which has a diameter of ∼36 (i.e., the "Einstein ring" formed by CO emission is likely to have a larger diameter compared to the optical one). This is consistent with the centroid position of the redshifted emission, which is along the quasar arc seen in the optical image, and the blueshifted emission, which is offset further to the SE (top right of Figure 2). Therefore, the CO-emitting region in RXJ1131 is likely to be more extended than its stellar and quasar emission.

We also place an upper limit on HNC(J = 2 → 1) line emission in the foreground galaxy at $z\sim 0.295$. Assuming a typical linewidth of 300 km s⁻¹, this corresponds to a 3σ limit of 0.35 Jy km s⁻¹ beam⁻¹.

We detect $\mathrm{CO}(J=3\to 2)$ line emission toward RXJ1131 in the CARMA data at ≳5σ significance. The $\mathrm{CO}(J=3\to 2)$ spectrum appears to be consistent with a double-peaked profile, as shown in Figure 3, where we over-plot spectra of the CO(J = 2 → 1) and $\mathrm{CO}(J=3\to 2)$ lines. We extract the peak fluxes and their corresponding uncertainties for the blue and red wing independently. We find a peak line flux of 5.13 ± 1.43 Jy km s⁻¹ beam⁻¹ for the blue wing, indicating a ≳3σ detection for this component alone, and a peak line flux of 11.45 ± 1.99 Jy km s⁻¹ beam⁻¹ for the red wing, indicating a ∼ 6σ detection. We measure a line intensity of 35.7 ± 6.9 Jy km s⁻¹ (Table 1) by summing up fluxes over the FWZI linewidth used to infer the CO(J = 2 → 1) line intensity.

Zoom In Zoom Out Reset image size

Assuming that the spatial extent of CO(J = 2 → 1) and $\mathrm{CO}(J=3\to 2)$ is similar and therefore the emission is magnified by the same amount, the measured line intensities correspond to a brightness temperature ratio of ${r}_{32}\,={T}_{\mathrm{CO}(J=3\to 2)}/{T}_{\mathrm{CO}}(J=2\to 1)=0.78\pm 0.37$. The quoted error bar is derived by adding the uncertainties associated with the CO line intensities and those from absolute flux calibrations in quadrature. This brightness temperature ratio is consistent with thermalized excitation within the uncertainties, as commonly observed in nuclear regions of nearby ULIRGs and high-z quasars (e.g., Weiß et al. 2007; Riechers et al. 2011b; Carilli & Walter 2013), but also with the lower excitation seen in normal star-forming disks (e.g., Dannerbauer et al. 2009; Carilli & Walter 2013; Daddi et al. 2015).

No 1.4 mm continuum emission is detected at the position of $\mathrm{CO}(J=3\to 2)$ down to a 3σ limit of 2.49 mJy beam⁻¹. This is consistent with the spectrum shown in Figure 3.

We detect 2.2 mm continuum emission at an integrated flux density of 1.2 ± 0.2 mJy, with a peak flux of ${S}_{\nu }=799\pm 88$ μJy beam⁻¹ centered on the lensing galaxy (Figure 4). Slightly extended emission along the lensing arc is also detected. This suggests that we detect emission in both the foreground and the background galaxy and that the emission is marginally resolved along its major axis. We subtract a point-source model in the visibility-plane to remove the unresolved part of the emission, which we here assume to be dominated by the foreground galaxy. The emission in the residual map coincides spatially with the lensing arc. We measure a flux density of ${S}_{\nu }=0.39\pm 0.08$ mJy for this residual component. This flux density is consistent with the difference between the integrated and the peak flux density measured in the original continuum map (∼0.4 mJy). We therefore adopt ${S}_{\nu }=0.39\pm 0.12$ mJy as the best estimate for the 2 mm continuum flux of the background galaxy (RXJ1131). We here quote a conservative error bar, which is derived by adding the uncertainty associated with the flux density of the point-source model ($\delta {S}_{\nu }=0.088$ mJy) with that of the peak flux in the residual map (0.08 mJy) in quadrature. We caution that this does not account for the systematic uncertainties of the de-blending procedure, where we have assigned 100% of the point-source flux to the foreground galaxy. We report the peak flux in the original map (${S}_{\nu }=799\pm 88$ μJy beam⁻¹) for the foreground galaxy, which is the best estimate possible at the resolution of our observations, but we acknowledge that a non-negligible contribution from the background source to the peak flux cannot be ruled out.

Zoom In Zoom Out Reset image size

The VLA C-band continuum image in Figure 4 shows resolved emission from the jets and the core of the foreground elliptical galaxy as well as emission toward the background quasar. Multiple peaks are seen along the arc with their centroids coincident with the optical emission from the quasar. We extract the flux densities for the lensing arc and the radio core in Table 2. We find a spectral index of ${\alpha }_{6\,\mathrm{cm}}^{2\,\mathrm{mm}}=-0.02\,\pm \,0.07$ for the foreground galaxy and ${\alpha }_{6\,\mathrm{cm}}^{2\,\mathrm{mm}}=-0.35\pm 0.21$ for the background galaxy by fitting power laws (${S}_{\nu }\propto {\nu }^{\alpha }$) to their continuum fluxes at 5 GHz and 2 mm. The spectral slope derived for the background source is flatter than the typical slope of pure synchrotron emission (α ∼ −0.7; e.g., Andreani et al. 1993). This likely suggests that at least a fraction of the observed 2 mm emission arises from thermal dust emission. This spectral slope would be even shallower if the background source contributes to the unresolved fraction of the 2 mm flux. In this case, the 2 mm flux of the foreground galaxy would be lower than the value reported here and lead to a slope steeper than ${\alpha }_{6\,\mathrm{cm}}^{2\,\mathrm{mm}}\sim -0.02$, which is flatter than that typical of elliptical galaxies. Assuming a spectral slope of α ∼ −0.7 to account for synchrotron radiation in RXJ1131, we expect a flux density of ${S}_{2\mathrm{mm}}=0.122\pm 0.004$ mJy at 2 mm. The flux excess of ${S}_{2\mathrm{mm}}=0.27\pm 0.08$ mJy therefore likely arises due to thermal dust emission.

Table 2.  Photometry Data

Wavelength	Frequency	Flux Density	Instrument
(μm)	(GHz)	(mJy)
Combined/Unresolved
1.25	239834	1.009 ± 0.090	CTIO/J-Band
1.65	181692	1.448 ± 0.121	CTIO/H-Band
2.17	138153	2.064 ± 0.160	CTIO/Ks-Band
3.4	88174.2	7.027 ± 0.142	WISE/W1
3.6	83275.7	5.618 ± 0.002	Spitzer/IRAC
4.5	66620.5	7.803 ± 0.002	Spitzer/IRAC
4.6	65172.3	8.872 ± 0.163	WISE/W2
5.8	51688.4	10.720 ± 0.005	Spitzer/IRAC
8.0	37474.1	14.470 ± 0.004	Spitzer/IRAC
12	24982.7	21.960 ± 0.425	WISE/W3
12	24982.7	<400	IRAS
22	13626.9	55.110 ± 1.878	WISE/W4
24	12491.4	70.204 ± 0.026	Spitzer/MIPS
25	11991.7	<500	IRAS
60	4996.54	<600	IRAS
100	2997.92	<1000	IRAS
250	1199.17	289.4 ± 9.6	Herschel/SPIRE
350	856.55	168.2 ± 8.6	Herschel/SPIRE
500	599.585	56.8 ± 8.8	Herschel/SPIRE
1387.93	216	<2.49	CARMA
2152.82	139.256	1.23 ± 0.22	PdBI
Foreground Lensing Galaxy (deblended bands)
0.555	540167	0.056 ± 0.006	HST-ACS/V-Band
0.814	368295	0.238 ± 0.013	HST-ACS/I-Band
1.6	187370	0.539 ± 0.041	HST-NICMOS(NIC2)/H-Band
3.6	83275.7	0.585 ± 0.003a	Spitzer/IRAC
4.5	66620.5	1.794 ± 0.003a	Spitzer/IRAC
5.8	51688.4	3.163 ± 0.006a	Spitzer/IRAC
8.0	37474.1	4.589 ± 0.006a	Spitzer/IRAC
2152.82	139.256	0.799 ± 0.082	PdBI
61414	4.8815	0.866 ± 0.027	VLA
Background Galaxy RXJ1131 (deblended bands)
0.555	540167	0.009 ± 0.004b	HST-ACS/V-Band
0.814	368295	0.041 ± 0.005b	HST-ACS/I-Band
1.6	187370	0.133 ± 0.004b	HST-NICMOS(NIC2)/H-Band
3.6	83275.7	5.034 ± 0.002	Spitzer/IRAC
4.5	66620.5	6.009 ± 0.002	Spitzer/IRAC
5.8	51688.4	7.557 ± 0.003	Spitzer/IRAC
8.0	37474.1	9.881 ± 0.004	Spitzer/IRAC
2152.82	139.256	0.39 ± 0.12c	PdBI
61414	4.8815	1.273 ± 0.042	VLA

Notes. The IRAC photometry for channel 1 (3.6 μm) is extracted directly from the image and from the Spitzer Heritage Archive for channels 2–4 (4.5, 5.8, and 8.0 μm). The flux uncertainties quoted for radio and millimeter observations (PdBI, CARMA, and VLA) do not include those from absolute flux calibration. All upper limits are 3σ.

^aFlux obtained using aperture photometry after subtracting the emission of RXJ1131 from the total emission. ^bA contribution from the quasar has been removed (see C06), and thus the flux density corresponds to the host galaxy only. ^cFlux extracted from the residual map after subtracting a point-source model. For SED modeling, we use ${S}_{2\mathrm{mm}}=0.27\pm 0.08$ mJy to exclude synchrotron emission (see Section 3.3).

References. The HST photometry is adopted from C06.

Download table as:  ASCIITypeset image

We compile mid-IR (MIR) to far-IR broadband photometry from various catalogs available on the NASA/IPAC Infrared Science Archive (IRSA) in Table 2 with aperture corrections when warranted. These data were obtained using the Cerro Tololo Inter-American Observatory (CTIO) for the Two Micron All Sky Survey (2MASS; Skrutskie et al. 2006), the Wide-field Infrared Survey Explorer (WISE; Wright et al. 2010), the Infrared Astronomical Satellite (IRAS; Neugebauer et al. 1984), and the Multiband Imaging Photometer (MIPS; Rieke et al. 2004) and Mid-infrared Infrared Array Camera (IRAC; Fazio et al. 2004) on the Spitzer Space Telescope . We retrieve PBCD (level 2) Spitzer/IRAC images from the Spitzer Heritage Archive and perform aperture photometry on the channel 1 image to extract the flux density at 3.6 μm since it is not available from the IRSA archive.

The emission in the IRAC images is slightly extended. We thus use an HST image (∼007 resolution) to determine the origin of their centroids, all of which are found to be centered at the position corresponding to the lensed emission from the background galaxy. To recover the diffuse background emission, we subtract a point-source model centered on the lensing galaxy, using the average FWHM found by fitting a Gaussian profile to several field stars with the imexam routine of IRAF. We perform aperture photometry on the residual image to obtain decomposed flux measurements of the background galaxy. The photometry for the foreground galaxy is then obtained by subtracting the background emission from the observed total flux. The resulting photometry in Table 2 is obtained after performing an aperture correction described in the IRAC Instrument Handbook³ to correct for the fact that the imaging was calibrated using a 12^'' aperture, which is larger than the aperture (58) we used to perform aperture photometry.

We fit a power-law spectrum to the decomposed IRAC photometry to disentangle the background and foreground emission from the total flux observed in the MIPS 24 μm band. The spectral indices corresponding to the best-fitting curves are α = −1.8 and α = −0.85 for the lensing galaxy and RXJ1131, respectively. The latter is consistent with the mean 3.6–8 μm spectral slope of α = −1.07 ± 0.53 found for unobscured AGN (Stern et al. 2005). An extrapolation of the fit to 24 μm yields 33.96 ± 0.01 mJy and 25.19 ± 0.03 mJy for the foreground galaxy and RXJ1131, respectively. The uncertainties are the standard deviations of the extrapolated fluxes obtained from two independent Monte Carlo simulations, each of 500 iterations. We incorporate the decomposed 24 μm data in our SED fitting to provide some constraints on the Wien tail beyond the dust peak of the SED of RXJ1131. Details of the SED modeling are presented in Section 4.5.

Extraction of the Herschel/SPIRE photometry at 250, 350, and 500 μm was carried out using sussextractor within the Herschel Interactive Processing Environment (HIPE; Ott 2010) on Level 2 maps obtained from the Herschel Science Archive. These maps were processed by the SPIRE pipeline version 13.0 within HIPE. The sussextractor task estimates the flux density from an image convolved with a kernel derived from the SPIRE beam. The flux densities measured by sussextractor were confirmed by using the Timeline Fitter, which performs photometry by fitting a 2D elliptical Gaussian to the Level 1 data at the source position given by the output of sussextractor. The fluxes obtained from both methods are consistent within the uncertainties.

At the angular resolution of the CO(J = 2 → 1) data, the source is resolved $\gtrsim 2$ resolution elements. Given the extent of the lensed emission (see Figure 2), this implies that we do not resolve structures (e.g., knots and arcs) of the lensed emission in our CO(J = 2 → 1) data. Nevertheless, the high spectral resolution of these data provides kinematic information on spatial scales smaller than the beam (see Figure 2). Hence, we reconstruct the intrinsic line profile and source-plane velocity structure by carrying out a parametric lens modeling over different channel slices of the interferometric data using our lensing code uvmcmcfit (Bussmann et al. 2015a; see Bussmann et al. 2015b for details of the code). Our approach follows a similar strategy as Riechers et al. (2008), who reconstruct a source-plane velocity gradient and constrain the gas dynamics in the z > 4 quasar host galaxy of PSS J2322+1944, which is also lensed into an Einstein ring configuration. To ensure adequate SNRs for lens modeling, we bin the frequency channels by a factor of five to produce seven independent ${\rm{\Delta }}v\sim \,105$ km s⁻¹ channels (dashed line in Figure 5) that cover the full linewidth of ∼750 km s⁻¹.

Zoom In Zoom Out Reset image size

We model the lens mass distribution using a singular isothermal ellipsoid (SIE) profile, which is described by five free parameters: the positional offset in R.A. and decl. relative to an arbitrary chosen fixed coordinate in the image, the Einstein radius, the axial ratio, and the PA. Positional offset between the foreground galaxy and the pre-defined coordinate is initialized using the VLA radio continuum map. We impose a uniform prior of ±005 in both ΔR.A. and Δdecl., motivated by the astrometry uncertainties in the VLA image as well as the uncertainties provided by previous SIE lens model (C06). We initialize the Einstein radius based on the model parameters reported by C06 and impose a uniform prior using ±3σ of their uncertainties. The sources are modeled using elliptical Gaussian profiles, which are parameterized by six free parameters: the positional offset in R.A. and decl. relative to the lens, the intrinsic flux density, the effective radius, the axial ratio, and the PA. The position of each source is allowed to vary between ±15 (i.e., within the Einstein radius) and the effective radius is allowed to vary from 001–2''.

Our code uses an Markov Chain Monte Carlo (MCMC) approach to sample the posterior probability distribution function (PDF) of the model parameters. In each model, we require a target acceptance rate of ∼0.25–0.5 and check for chain convergence by inspecting trace plots and by requiring that the samples are obtained beyond at least an autocorrelation time. We thus employ ∼50,000 samples as the initial "burn-in" phase to stabilize the Markov chains (which we then discard) and use the final ∼5000 steps, sampled by 128 walkers, to identify the posterior. Here, we identify the best-fit model and the quoted uncertainties using the median and the 68% confidence intervals in the marginal PDFs.

We first obtain a preliminary lens model for each channel slice independently, where their lens parameters are allowed to vary and are initialized according to the aforementioned way. We obtain the final model by repeating the modeling over each slice but fixing their lens parameters to the overall median in the preliminary models, as listed in Table 3. This ensures that all models share the same lens profile. The magnification factors in Table 4 are determined by taking the ratio between the image-plane flux and the source-plane flux of each model.

Table 3.  Lens Parameters Constrained by Models of Seven Velocity Channels

Parameters		Median Values
Offset in R.A.	('')	0.004 ± 0.027
Offset in decl.	('')	0.003 ± 0.027
Axial Ratio		0.56 ± 0.16
Position Angle	(deg)	103 ± 22
Einstein Radiusa	('')	1.833 ± 0.002

Notes. Parameters describing the foreground lens are obtained based on the median in the preliminary models (see the text for details). All angular offsets are with respect to α = 11^h31^m5144, δ = −12°31'583 (J2000).

^aThis corresponds to a mass of $M(\theta$ < θ_E) = (7.42 ± 0.02) × 10¹¹ M_⊙ within the Einstein radius.

Download table as:  ASCIITypeset image

Table 4.  Magnification Factors of Various Kinematic Components in CO(J = 2 → 1)

Velocity Range (km s−1)	Source 1 ${\mu }_{{\rm{L}}}$	Source 2 ${\mu }_{{\rm{L}}}$
346–260	6.7 ± 2.5	7.2 ± 5.6
238–153	7.6 ± 1.6	⋯
131–45	8.7 ± 2.0	⋯
24 to −62	4.1 ± 0.9	⋯
−84 to −170	4.2 ± 0.6	⋯
−191 to −277	4.3 ± 2.4	⋯
−300 to −385	3.1 ± 0.9	⋯
weighted average	4.4	⋯
median	5.5	⋯

Note. First column corresponds to the rest-frame velocity ranges taken from the center of an unbinned channel (see Figure 5). Each row corresponds to a (binned) channel slice used for lens modeling. Source 1 is RXJ1131 and source 2 is its companion.

Download table as:  ASCIITypeset image

Our model parameters in Table 3, describing the mass distribution of the lensing galaxy, are consistent (within the uncertainties) with that of the SIE model presented by C06. We find a mass of $M(\theta \lt {\theta }_{{\rm{E}}})=(7.47\pm 0.02)\times {10}^{11}$ M_⊙ within the Einstein radius.

4.1.1. Interpretation of the Source-plane Morphology

The reconstructed source locations, as represented by magenta ellipses in Figure 6, demonstrate an intrinsic velocity gradient across the source plane, which is consistent with a kinematically ordered disk-like galaxy. Additional support to the disk conjecture can be found in the double-horned line profile (Figure 1) and the observed (image plane) velocity field (Figure 2). Furthermore, C06 also find that the reconstructed source-plane emission in optical-NIR is best-reproduced using a n = 1 Sérsic profile. We thus interpret RXJ1131 as a disk galaxy.

Zoom In Zoom Out Reset image size

A better fit is found for the lens model of the red-most channel if we add a second source component (see the top left panel in Figure 6). This is consistent with previous results reported by Brewer & Lewis (2008, hereafter B08), who find an optically faint companion (component F in their paper) ∼2.4 kpc in projection from the AGN host galaxy in V-band, and with C06, who find evidence for an interacting galaxy near RXJ1131. Spatially, the red velocity component of the CO emission is also consistent with this component F. It is therefore likely that we detect CO(J = 2 → 1) emission in a companion galaxy.

We decompose the total line flux into two components: one from RXJ1131 and the other from its companion. Since the companion is only detected in the red-most channel, we derive its intrinsic gas mass using the best-fit flux densities and magnification factors obtained from the models of this channel. Assuming a brightness temperature ratio of ${r}_{21}=1$ between CO(J = 2 → 1) and CO(J = 1 → 0) lines and a CO luminosity-to-H₂ mass conversion factor of ${\alpha }_{\mathrm{CO}}$  = 0.8 ${M}_{\odot }\,{({\rm{K}}\mathrm{km}{\mathrm{pc}}^{2})}^{-1}$, we find a molecular gas mass of ${M}_{\mathrm{gas}}=(1.92\pm 0.09)\times {10}^{9}$ M_⊙. For the molecular gas mass in RXJ1131, we derive its intrinsic line flux over the FWZI linewidth using the respective magnification factors listed in Table 4, which to first order takes into account the effect of differential lensing. This yields ${I}_{\mathrm{CO}(J=2\to 1)}=2.93\,\pm 0.70$ Jy km s⁻¹, where the uncertainty includes those on the magnification factors. Adopting the same brightness temperature ratio and ${\alpha }_{\mathrm{CO}}$ as used for the companion, this corresponds to a gas mass of ${M}_{\mathrm{gas}}=(1.38\pm 0.33)\times {10}^{10}$ M_⊙, which implies a gas mass ratio of ∼7:1 between RXJ1131 and its companion.

The spatial resolution of the data in hand is a few arcseconds, which implies that despite the high SNR and spectral resolution, constraints on the intrinsic sizes of the lensed galaxies are modest, and thus the magnification factors may be under-predicted (see, e.g., Bussmann et al. 2015b; Dye et al. 2015; Rybak et al. 2015).

4.1.2. Spatial Extent and Differential Lensing

In the image-plane integrated line map shown in Figure 2, the redshifted component is cospatial with the Einstein ring that is seen in the optical image, with most of its apparent flux originating along the lensing arc, whereas the centroid of the blueshifted emission is offset to the SE of the lensing arc. This suggests that the CO-emitting region in RXJ1131 is extended. To further illustrate this, we show the channel maps of 21.5 km s⁻¹ width and a spatial spectra map of 15 resolution in Figures 7 and 8, respectively. These figures show that redshifted emission is present to the west, peaking toward the lensing arc (black crosses in Figure 7), and shifts to the east with decreasing velocity (blue wing). This is consistent with the source-plane positions in our models and is suggestive of an extended CO-emitting region.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Previous studies of RXJ1131 find evidence for differential lensing across the HST V-, I-, and H-bands, where the magnification factor varies from 10.9 to 7.8 (C06). This indicates that the emission from different stellar populations within the host galaxy have various spatial extents and positions with respect to the caustic. The best-fit lens models obtained here for different CO channels show that differential lensing also plays an important role in the observed CO(J = 2 → 1) emission, with a magnification factor $({\mu }_{{\rm{L}}})$ that varies from 8.7 to 3.1 across different kinematic components (Table 4). The asymmetry in the line profile (Figures 1 and 5) is therefore predominantly a result of the redshifted CO-emitting gas being more strongly magnified than the blueshifted component. A secondary reason is likely due to the inclusion of the emission of the companion in the most redshifted velocity channels. The variation in ${\mu }_{{\rm{L}}}$ found across channels is consistent with the source-plane positions relative to the caustics in Figure 6, where the red wing emission mainly originates near the cusp of the caustic and the blue wing emission is located beyond the caustics. In fact, the intrinsic line flux of the redshifted and blueshifted emission in RXJ1131 (after subtracting a contribution from the companion) are ${I}_{\mathrm{CO}(J=2\to 1)}=1.26\pm 0.23$ Jy km s⁻¹ and 1.25 ± 0.23 Jy km s⁻¹, respectively, implying an intrinsically symmetric line profile (Figure 5). This is consistent with the symmetric source-plane velocity gradient in our lens model (Figures 6 and 9).

Zoom In Zoom Out Reset image size

Fitting two Gaussians with a common FWHM to the "intrinsic" CO(J = 2 → 1) line profile of RXJ1131 (after correcting for lensing using the magnification factors for various channels and separating the emission from RXJ1131 and its companion), we find a roughly symmetric double-horned profile with a flux ratio of 1.2 ± 0.4 between the peaks, which are separated by ${\rm{\Delta }}{v}_{\mathrm{sep}}=387\pm 45$ km s⁻¹, and each with an FWHM of 220 ± 72 km s⁻¹. The peak separation obtained from this "intrinsic" line profile is slightly lower than that obtained from the observed spectrum (i.e., without lensing corrections). This discrepancy is likely a result of differential lensing, which causes the line peak of the red wing to shift toward higher velocity channels, and thereby biasing the centroid of one Gaussian to a higher velocity than otherwise. To facilitate a comparison (Section 5.1.1) with previous works, which were observed at lower spectral resolution, we also fit a single Gaussian to the intrinsic line profiles. This yields FWHMs of 600 ± 160 km s⁻¹ for RXJ1131 and 73 ± 43 km s⁻¹ for the companion galaxy.

A clear velocity gradient and a high velocity dispersion (≳400 km s⁻¹) near the central region are seen in Figure 2. While beam smearing is inevitably the dominant factor in the observed velocity dispersion at the spatial resolution of these data, the exceedingly high velocity dispersion may hint at potential perturbations from the AGN, or internal turbulence due to interactions with the companion, and/or instability due to the large gas content. Therefore, in this scenario, RXJ1131 is consistent with a disrupted disk galaxy hosting an optically bright quasar and in the process of merging.

As discussed in Section 4.1.1, we interpret RXJ1131 as a disk galaxy since it displays a kinematically-ordered velocity gradient in the source-plane velocity map of the CO emission, a symmetric double-horned line profile (Figures 5, 6 and, left panel of Figure 9), and a disk-like morphology in the source-plane reconstruction of the optical-NIR emission (C06). We extract a one-dimensional position–velocity (PV) profile by assuming that the source-plane centroids of different velocity components obtained from dynamical lens modeling are dominated by the tangential component of the true velocity vector of a rotating disk, i.e., each velocity component would be seen as lying along the major axis of a rotating disk if observed with sufficiently high angular resolution (see the right panel of Figure 9). In this process, the positions for each velocity component (plotted as data points in the right panel of Figure 9) are extracted along the best-fitted major axis, which is along a PA of 121°.

We then attempt to characterize the molecular gas kinematics using an empirically motivated disk model (e.g., Courteau 1997; Puech et al. 2008; Miller et al. 2011):

$$\begin{eqnarray}&&V={V}_{0}+\displaystyle \frac{2}{\pi }{V}_{a}\arctan \left(\displaystyle \frac{R}{{R}_{t}}\right),\end{eqnarray} \tag{ 1 }$$

where V is the observed velocity, V₀ is the velocity at dynamical center, V_a is the asymptotic velocity, and R_t is the "turnover" radius at which the rising part of the curve begins to flatten. We perform non-linear least squares fitting using an orthogonal distance regression to find the best-fit parameters, taking into account the uncertainties in both velocity (channel width) and distance offset. We also place an upper limit on R_t < 15 kpc to keep this parameter physical (e.g., Puech et al. 2008; Miller et al. 2011). The parameter uncertainties are inferred based on a Monte Carlo simulation of 500 iterations, where the input parameters are perturbed according to random Gaussian distributions with standard deviations corresponding to their uncertainties. Using this model, we find V_a = 988 ± 618 km s⁻¹, R_t = 10.9 ± 7.8 kpc, and V₀ = 0 ± 9 km s⁻¹. However, since the emission is not resolved along the flat regime of the rotation curve, the asymptotic velocity and the "turnover" radius are poorly constrained. In particular, V_a and R_t are highly correlated with a Pearson coefficient R = 0.998, and −0.400 between V_a and V₀.

The asymptotic velocity (V_a)—an extrapolation of the model out to radius beyond the disk scale-length and half-light radius—is not equivalent to the maximum observed velocity (${V}_{\max }$), which is commonly used in literature to parameterize disk rotation. The arctangent model is most commonly used in studies of the Tully–Fisher relation, where an extrapolation to ${V}_{2.2}$ (velocity at 2.2 disk scale-length or ∼1.375 half-light radius, or $\sim 0.7\,{R}_{\mathrm{opt}}$ ⁴ ) is typically adopted as the rotation velocity (${V}_{\max }$ in their terminology), since this corresponds to the radius at which the velocity of a pure exponential disk peaks (Courteau & Rix 1997). Here, we adopt the maximum observed velocity ${V}_{\mathrm{rot}}=303\pm 55$ km s⁻¹ at 6 ± 3 kpc from the dynamical center as a proxy to the rotation velocity. This radius corresponds to ∼0.6 R_e, where R_e is the half-light radius ∼10.3 kpc inferred from the HST I-band lens model (C06; converted to our cosmology). We note that the source-plane half-light radius varies substantially with wavelength. In particular, the half-light radius is found to be ∼ 4 kpc and ∼7 kpc in V-band (B08) and H-band (C06), respectively. The CO gas is thus of similar spatial extent as in the H- and I-bands.

In the rest-frame, emission in the observed-frame H-band corresponds to NIR emission (∼1 μm), tracing radiation from the accretion disk surrounding the central AGN and also from old and evolved stellar populations; I-band corresponds to roughly the optical V-band, tracing stellar radiation from existing, less massive (i.e., longer-lasting) stars; V-band corresponds to roughly U-band, tracing radiation from massive young stars in the host galaxy. Hence, the relative compactness observed in the V-band may be explained in part due to the fact that the emission in this band is more susceptible to dust extinction than in other bands and/or dominated by a central starburst caused by higher concentrations of star-forming gas toward the central regions—owing to gravitational perturbations induced from interactions with the companion (e.g., Di Matteo et al. 2005). This would be consistent with the picture that old stars form first and constitute the bulge component of a spiral galaxy, and that nuclear starbursts (in the inner few kiloparsecs) can be triggered at a later time as the progenitor disk galaxy interacts with other galaxies to form a larger bulge.

Assuming the gas to be virialized, the dynamical mass can be approximated by ${M}_{\mathrm{dyn}}\sim {\sigma }^{2}R/G$, where σ is the velocity dispersion, or the rotational velocity in the case of a rotating disk model (i.e., $\sigma ={V}_{\mathrm{rot}}\,\sin i)$. Using a rotational velocity ${V}_{\mathrm{rot}}\,\sin i=303$ km s⁻¹ (see Section 4.3), we find a dynamical mass of ${M}_{\mathrm{dyn}}$ ${\sin }^{2}i$ (<6 kpc) = 1.3 × 10¹¹ M_⊙ enclosed within the CO-emitting region in RXJ1131. If we instead consider the CO(J = 2 → 1) line peak separation $({\rm{\Delta }}{v}_{\mathrm{sep}}/2\sim 200$ km s⁻¹) as the rotation velocity, we find ${M}_{\mathrm{dyn}}$ ${\sin }^{2}i$ (<6 kpc) = 5.8 × 10¹⁰ M_⊙. We derive an inclination angle of 564 from the morphological axial ratio of a/b ∼ 18/325, which we estimate from the source-plane image reconstructed by C06 (Figure 3 in their paper). This corresponds to an inclination-corrected dynamical mass of 8.3 × 10¹⁰ M_⊙ < ${M}_{\mathrm{dyn}}$ < 25 × 10¹⁰ M_⊙. Our estimate should be considered at best an upper limit since the gas in RXJ1131 is unlikely to be virialized. In the following sections, we use the lower limit (8.3 ± 1.9) × 10¹⁰ M_⊙ as the dynamical mass since it is derived in a manner similar to what is commonly used in the literature (e.g., Solomon et al. 1997, hereafter S97; Downes & Solomon 1998, hereafter DS98; Greve et al. 2005, hereafter G05).

Using the velocity dispersion (σ = 30 km s⁻¹) obtained by fitting a single Gaussian to the de-lensed line profile of the companion and a half-light radius of ${R}_{\mathrm{CO}}=4.2\pm 2.8$ kpc from the best-fit lens model, we find a dynamical mass of M_dyn = (3.5 ± 2.3) × 10⁹ M_⊙ for the companion, assuming an inclination angle of i = 30°. The uncertainty here only includes that of the CO source size. On the other hand, we find ${M}_{\mathrm{dyn}}$ ${\sin }^{2}30^\circ =5.8\,\times \,{10}^{8}$ M_⊙ if we adopt the better-constrained V-band source size of ∼700 pc (B08). Since the V-band based dynamical mass measurement is substantially lower than the gas mass, the V-band emitting region may appear to be much smaller than its true extent due to dust obscuration.

The CO-based dynamical mass estimates correspond to a mass ratio of ∼24:1 between RXJ1131 and the companion, with a gas mass ratio of ∼7:1 derived in Section 4.1.1. We thus classify the system as a gas-rich,"wet" minor merger.

We fit dust SED models to the 24 μm−2.2 mm photometry using a modified-blackbody (MBB) function with a power-law attached to the Wien side to account for the MIR excess due to emission of warm and small dust grains. The IRAS 60 and 100 μm upper limits are included to constrain the dust peak. Here, we use a flux density of ${S}_{2\mathrm{mm}}=0.27\pm 0.08$ mJy derived in Section 3.3 instead of the deblended flux listed in Table 2 to exclude a potential contribution due to synchrotron emission (see Section 3.3) in the dust SED modeling. An uncertainty from absolute flux calibration of ∼15% is added in quadrature to the PdBI 2 mm continuum photometry in our fitting procedure.

The fit is performed using the code mbb_emcee (e.g., Riechers et al. 2013; Dowell et al. 2014), which samples the posterior distributions using an MCMC approach and uses instrumental response curves to perform color correction. The model is described by five free parameters: the rest-frame characteristic dust temperature (T_d), the emissivity index (β), the power-law index (α), the flux normalization at 500 μm (${f}_{\mathrm{norm}}$), and the observed-frame wavelength at which the emission becomes optically thick (${\lambda }_{0}$). We impose a uniform prior with an upper limit of 100 K on T_d (see e.g., Sajina et al. 2012), a Gaussian prior centered around 1.9 with a standard deviation of 0.3 on β, and a uniform prior with an upper limit of 1000 μm on ${\lambda }_{0}$. We check for chain convergence by requiring that the autocorrelation length of each parameter is less than the number of steps taken for the burn-in phase (which are then discarded). Here we report the statistical means and the 68% confidence intervals in the marginal PDFs as the best-fit parameters, as listed in Table 5. The best-fit models are shown in Figure 10 along with the broadband photometry that is listed in Table 2.

Zoom In Zoom Out Reset image size

Table 5.  SED Fitting Results

Parameters		With 24 μm	Without 24 μm
Td	(K)	54${}_{-10}^{+8}$	55${}_{-20}^{+21}$
β		1.6${}_{-0.4}^{+0.5}$	2.2${}_{-0.3}^{+0.3}$
α		1.6${}_{-0.6}^{+0.5}$	8.5${}_{-6.2}^{+7.0}$
${\lambda }_{0}$ a	(μm)	559${}_{-324}^{+278}$	365${}_{-120}^{+111}$
${\lambda }_{\mathrm{peak}}$ b	(μm)	159${}_{-40}^{+19}$	155${}_{-43}^{+38}$
${f}_{\mathrm{norm},500\mu {\rm{m}}}$ c	(mJy)	55${}_{-13}^{+13}$	59${}_{-6}^{+6}$
${L}_{\mathrm{FIR}}$ d	(1012 L⊙)	3.81${}_{-1.97}^{+1.92}$	4.24${}_{-2.00}^{+2.17}$
${M}_{{\rm{d}}}$ e	(108 M⊙)	16${}_{-12}^{+5}$	14${}_{-7}^{+5}$

Notes. Errors reported here are ±1σ. ${L}_{\mathrm{FIR}}$ and ${M}_{{\rm{d}}}$ are not corrected for lensing.

^aObserved-frame wavelength, where ${\tau }_{\nu }=1$. ^bObserved-frame wavelength of the SED peak. ^cObserved-frame flux density at 500 μm. ^dRest-frame 42.5–122.5 μm luminosity. ^eDerived assuming an absorption mass coefficient of κ = 2.64 m² kg⁻¹ at λ = 125.0 μm (Dunne et al. 2003).

Download table as:  ASCIITypeset image

In the first model, we attempt to constrain the power-law index by including the 24 μm data. Based on the resulting posterior PDFs, we find an apparent IR luminosity (rest-frame 8–1000 μm) of 8.22${}_{-2.98}^{+2.75}$ × 10¹² L_⊙, a far-IR luminosity (rest-frame 42.5–122.5 μm) of 3.81${}_{-1.97}^{+1.92}$ × 10¹² L_⊙, and a dust mass of 16${}_{-12}^{+5}$ × 10⁸ M_⊙, none of which are corrected for lensing magnification. For the mass absorption coefficient, we adopt κ = 2.64 m² kg⁻¹ at rest-frame 125.0 μm (Dunne et al. 2003). The dust mass uncertainty does not include that of the absorption coefficient.

A fit including the MIR 24 μm photometry is likely an upper limit on the far-IR luminosity due solely to star formation in the AGN host galaxy. If we instead fit for a model excluding this constraint, two major consequences are immediately apparent. First, the power-law index is poorly constrained (see Table 5). Second, the steep power-law implies only a small contribution from the power-law regime to the total IR luminosity as compared to the graybody component. Thus, the far-IR luminosity in this model should, in principle, correspond to a lower limit on the cold dust emission. Using the best-fit parameters for this model, we find a total IR luminosity ${L}_{\mathrm{IR}}$(rest-frame 8–1000 μm) of 8.67${}_{-5.27}^{+5.27}$ × 10¹² L_⊙, a far-IR luminosity ${L}_{\mathrm{FIR}}$ of 4.24${}_{-2.00}^{+2.17}$ × 10¹² L_⊙, and a dust mass ${M}_{\mathrm{dust}}$ of 14${}_{-7}^{+5}$ × 10⁸ M_⊙, all of which are not lensing-corrected. Taken at face value, this implies an FIR-to-IR luminosity ratio of ∼49 ± 38%.

The dust temperature from both models is similar to that of ULIRGs at 0.6 < z < 1.0 (54 ± 5 K; Combes et al. 2013, hereafter C13). The far-IR luminosity is comparable in both models, which is not surprising given the lack of constraints in the MIR. For the subsequent analysis, we adopt the physical quantities from the first model (i.e., with constraints at 24 μm). The choice of SED model does not affect the derived SFR given the similar far-IR luminosity, and their dust masses are consistent within the uncertainties. We correct for lensing using the median magnification factor $({\mu }_{{\rm{L}}}=5.5)$ from the CO lens models. This yields a ${L}_{\mathrm{FIR}}$ of (6.9 ± 3.6) × 10¹¹ L_⊙ and an intrinsic total IR luminosity of ∼1.5 × 10¹² $(5.5/{\mu }_{{\rm{L}}})$ L_⊙, implying that RXJ1131 can be classified as a ULIRG. Assuming a Salpeter (1955) initial mass function (IMF), we find an SFR${}_{\mathrm{FIR}}$ of 120 ± 63 M_⊙ yr⁻¹ using a standard conversion (Kennicutt 1998).

We derive the stellar mass of RXJ1131 by fitting SED models to the rest-frame UV-to-millimeter photometry using the high-z version of the magphys code (da Cunha et al. 2008, 2015). Two sets of stellar templates modeled using either the Bruzual & Charlot (2003) or the unpublished Bruzual (2007) stellar population synthesis code are provided in the magphys package. We adopt the former set. To minimize contaminations from the quasar, we only fit to the HST, Herschel, and PdBI data, where both the HST and the PdBI 2 mm photometry are deblended from the AGN (see bottom section of Table 2). The input photometry are corrected for lensing using their respective magnification factors to account for differential lensing (light blue circles in Figure 10). We thus find a stellar mass of M_* = 2.95${}_{-0.86}^{+1.32}$ × 10¹⁰ M_⊙, which is the median value of the posterior probability distribution and the uncertainties are derived from the 16th and 84th percentiles. We note that the models are over-fitted with a best-fit of ${\chi }^{2}=0.41$, which is unsurprising due to sparse sampling of the SED compare to the number of free parameters. The resulting dust mass and IR luminosity are consistent with those obtained from the MBB+power-law models within the uncertainties, albeit some differences in the assumptions behind the two methods. The consistency may be attributed to the large uncertainties arising from the lack of photometric constraints on the models and the fact that the best-fit parameters from the MBB method are similar to those of the magphys method.

In this section, we derive the gas properties of the merging system RXJ1131 based on CO(J = 2 → 1) and compare them with those reported by C13 ⁵ —the largest sample of CO-detected ULIRGs at similar redshift (0.6 < z < 1.0). Their results are based on spatially unresolved CO(J = 2 → 1) and CO(J = 4 → 3) line observations with the IRAM 30-m single-dish telescope.

5.1.1. Linewidths and Sizes

5.1.2. Gas Mass Fractions and Gas-to-dust Ratio (GDR)

5.1.3. Star Formation Efficiency and Specific SFR

Sluse et al. (2007) report two sets of AGN lines observed in RXJ1131. The first set of lines is at z ∼ 0.654, including the narrow component of the Balmer lines, the [O iii] 4959, 5007 Å lines, and the Mg ii 2798 Å absorption line; the second set is at ${z}_{{\rm{s}},\mathrm{QSO}}\,\sim \,0.658$, including the broad component of the Balmer lines and the Mg ii 2798 Å emission line. Using the CO line center redshift as the systemic redshift, we find that the redshift of the first set is fully consistent with the systemic redshift. This supports previous claims that [O iii] lines, tracing the narrow line region (NLR), are good proxies to the true systemic redshift (e.g., Vrtilek 1985; Nelson 2000). On the other hand, the second set of lines is redshifted by ∼715 km s⁻¹.

Velocity offsets between broad line region (BLR) and NLR lines have been reported in the literature. Richards et al. (2002) find a median offset of ∼100 ± 270 km s⁻¹ between [Mg ii] and [O iii] lines in a sample of >3800 quasars, and Bonning et al. (2007) report a mean offset of ∼100 ± 210 km s⁻¹ between the broad component of Hβ and [O iii] lines in a sample of ∼2600 quasars at 0.1 < z < 0.8, where only $\lesssim 20$ of them (i.e., <1%) are found to have offsets >800 km s⁻¹ and ∼1% are found to have offsets >500 km s⁻¹.⁷ Thus, large velocity offsets between BLR and NLR lines comparable to that of RXJ1131 are uncommon but have been observed in some cases.

The observed velocity offset between the BLR and NLR lines may be explained by a recoiling black hole (BH), where the BLR is moving at high velocity relative to the bulk of its host galaxy (Madau & Quataert 2004; Bonning et al. 2007; Loeb 2007). Depending on the initial conditions of the black hole pair (e.g., black hole mass ratio, spin–orbit orientation, spin magnitude), numerical relativity simulations have shown that recoil velocities can reach up to ${v}_{\mathrm{kick}}\sim 4000$ km s⁻¹ for spinning BHs, with typical recoil velocities of ${v}_{\mathrm{kick}}\,\sim 100\mbox{--}500$ km s⁻¹ (e.g., Libeskind et al. 2006; Campanelli et al. 2007). Several sources have been proposed as recoiling BH candidates (Komossa et al. 2008; Civano et al. 2010; Steinhardt et al. 2012). However, Decarli et al. (2014) have recently refuted such scenarios for one of the candidates—SDSS J0927+2943—by finding that the redshift of its BLR lines is indeed consistent with its CO systemic redshift. This is in contrast with RXJ1131, where our CO observations confirm the redshifted BLR lines compared to the CO systemic redshift. Since this scenario requires a coalesced BH, it would imply that RXJ1131 is a product of a previous merger, which is not implausible and might also explain the highly spinning BH in RX1131 (a ∼ 0.9; Reis et al. 2014).

Alternative scenarios, e.g., outflow/inflow of gas in the BLR, viewing angle toward the accretion disk, and obscuration in the clumpy accretion disk are more commonly invoked to explain velocity offsets between BLR and systemic redshift. Since the BLR lines of RXJ1131 show positive velocity offsets with respect to its systemic redshift, it may imply that the observed BLR line emission is dominated by the gas that is flowing into the central BH, or by the receding component of the accretion disk, owing to the viewing angle or the obscuration in the accretion disk. Sluse et al. (2007) report a covering factor of 20% for the accretion disk in RXJ1131 based on its broad Mg ii 2798 Å absorption line at z = 0.654. Additionally, the centroids of the BLR lines in RXJ1131 may be biased toward longer wavelengths due to microlensing (e.g., Sluse et al. 2007, 2012), which may have magnified the redshifted component of the compact BLR more strongly than its blueshifted component.

We find a ${M}_{\mathrm{BH}}$/${M}_{\mathrm{bulge}}$ ratio of >0.27${}_{-0.08}^{+0.11}$% using the black hole mass of ${M}_{\mathrm{BH}}\sim 8\times {10}^{7}$ M_⊙ (Sluse et al. 2012) and the stellar mass derived in Section 4.5 as an upper limit to the bulge mass. This ratio is consistent with those of other intermediate-z radio-loud AGNs (McLure et al. 2006) but is higher than those of nearby AGNs (Häring & Rix 2004). Our results therefore support the emerging picture that quasars grow faster and/or earlier than their host galaxies at higher redshifts (e.g., Walter et al. 2004; McLure et al. 2006; Peng et al. 2006; Riechers et al. 2008). The elevated ${M}_{\mathrm{BH}}$/${M}_{\mathrm{bulge}}$ ratio of RXJ1131 compared to local AGNs suggests that the bulk of the black hole mass of RXJ1131 is largely in place while its stellar bulge is still assembling.