Morpho-kinematics of the molecular gas in a quasar host galaxy at redshift z=0.654
T.T. Thai, P. Tuan-Anh, P. Darriulat, D.T. Hoai, P.T. Nhung, P. N. Diep, N.B. Ngoc, N.T. Phuong
aa r X i v : . [ a s t r o - ph . GA ] J a n Submit to Communications in Physics, (2020), pp. 1-20DOI:10.15625/0868-3166/0/0/0
MORPHO-KINEMATICS OF THE MOLECULAR GAS IN A QUASAR HOSTGALAXY AT REDSHIFT z =0.654 T.T.THAI , , P. TUAN-ANH , P. DARRIULAT , D. T. HOAI , P. T. NHUNG , P. N. DIEP , N.B. NGOC , AND
N. T. PHUONG Vietnam National Space Center, Vietnam Academy of Science and Technology18 Hoang Quoc Viet, Cau Giay, Hanoi, Vietnam Graduate University of Science and Technology18 Hoang Quoc Viet, Cau Giay, Hanoi, Vietnam † E-mail: [email protected]
Received January 5, 2021Accepted for publication ??Published ??
Abstract.
We present a new study of archival ALMA observations of the CO(2-1) line emissionof the host galaxy of quasar RX J1131 at redshift z=0.654, lensed by a foreground galaxy. Asimple lens model is shown to well reproduce the optical images obtained by the Hubble SpaceTelescope. Clear evidence for rotation of the gas contained in the galaxy is obtained and a simplerotating disc model is shown to give an excellent overall description of the morpho-kinematics ofthe source. The possible presence of a companion galaxy suggested by some previous authors isnot confirmed. Detailed comparison between model and observations gives evidence for a morecomplex dynamics than implied by the model. Doppler velocity dispersion within the beam size inthe image plane is found to account for the observed line width.
Keywords: galaxies: evolution – galaxies: ISM – radio lines: galaxies . Classification numbers: 98.65-r. I . INTRODUCTIONI.1 . General features
RX J1131-1231 (simply called RX J1131 in the following), is a distant quasar, at redshift z s ∼ ∼ ∼ ©2020 Vietnam Academy of Science and Technology T.T.Thai, P. Tuan-Anh & P. Darriulat with a mass of ∼ M ⊙ ; it rotates extremely fast, reaching near half the light velocity [3].The quasar and its host galaxy are gravitationally lensed by a galaxy in the foreground, at redshift z L ∼ Fig. 1.
Left: dependence on the redshifts of the source ( z S in abscissa) and of the lens ( z L in ordinate) of the ratio between their respective angular diameter distances d aS / d aL [1,2]).The relative size of the lens with respect to the source is proportional to d aS / d aL . The starsshow the locations of quasars RX J1131 (black, P18) and RX J0911 (red, [8], [9]). Thesizes of the host galaxies are compared to the size of the lens caustic in the central andright panels respectively. tion (SED), and archival VLA observations (Program ID: AW741; PI: Wucknitz) analysed byLeung et al. 2017 [11], referred to as L17 in the following, have shown resolved continuum emis-sion from the jets and the core of the foreground elliptical galaxy as well as emission toward thebackground quasar.Thorough analyses of high angular resolution HST optical and NIR images [12–14] and ofKeck Adaptive Optics images [15] have produced a detailed description of the lensing propertiesin the neighbourhood of the quasar. They reveal a typical long axis quad configuration [16,17], thequasar being located within the eastern cusp of the lens caustic. As emission from the lens galaxyis simultaneously detected, the parameters of the lensing potential can be accurately evaluated.However, they probe only the vicinity of the cusp of the caustic curve. As the emission of thequasar host galaxy covers the whole caustic curve and extends even farther out, one cannot takeit as granted that the simple lens model obtained from the study of the quasar images reliablyapplies to the whole host galaxy. This is at variance with the gravitational lensing of quasar hoststhat are farther away and cover only part of the caustic in addition to being intrinsically smaller.The central region of the caustic corresponds to images close to an Einstein ring configuration,which dominates the picture in the case of RX J1131. We illustrate this feature in Figure 1, whichshows the location of RX J1131 in the plane z L vs z S , lens vs source redshifts, together with that ofother multiple imaged systems [12], and compares it with the case of a typical farther away quasarhost galaxy, RX J0911 [8, 9]. orpho-kinematics of quasar host galaxy RXJ 1131 3 I.2 . Millimeter observations
The present work uses archival ALMA observations to study the emission of the CO(2-1)molecular line of the quasar host galaxy, with a beam of ∼ × providing unprece-dented image quality. The data have been analysed in much detail by the team who proposed theobservation ( [18], referred to as P18 in the following) and who kindly sent us data files sum-marizing the results of their analysis. The CO(2-1) line emission of RX J1131 was observed byL17 at the Plateau de Bure Interferometer (PdBI) with an angular beam size (FWHM) of 4.4 × , a spectral resolution of ∼ − and a noise level of ∼ − per chan-nel. In addition to their CO(2-1) PdBI observations, L17 used the Combined Array for Researchin Millimeter-wave Astronomy (CARMA) to detect the CO(3-2) line emission but the signal tonoise ratio is at only 5- σ significance.ALMA observations of both 2 mm continuum and CO(2-1) line emission have been anal-ysed by P18 with an angular resolution of ∼ σ ) line emission extendedover a complete Einstein ring. The map of velocity dispersion displays values covering from ∼ ∼
50 km s − . The authors of P18 note that peaks in the intensity map and the region of highvelocity dispersion are coincident and probably reveal on-going star-formation; these peaks arenot strictly coincident with the quasar emission.As mentioned above, one cannot take it as granted that the simple lens model obtainedfrom the study of the quasar images reliably applies to the whole host galaxy. For this reason P18reconstruct the source brightness distribution using the ALMA CO(2-1) data exclusively [19–21].The result (Figure 1, central panel) is consistent with a large rotating disc inclined by 54 ◦ withrespect to the plane of the sky, having rotation velocity reaching over 400 km s − .The model used by P18 to describe the lens includes an external shear, ellipticity of the mainlens and a contribution from the small satellite galaxy revealed by the HST image north of the mainlens (central panel of Figure 2). The contribution of the latter is known to be very small [13] andthe contribution of ellipticity, while different from that of an external shear, does not affect stronglythe general picture: excellent descriptions of the lensing of the quasar point source are obtainedwith either external shear alone or ellipticity alone and their combination in the P18 model causesa rotation of only 4 ◦ of the caustic with respect to a model ignoring ellipticity. For this reason ourapproach in the present article is to use a simple lens potential with only external shear and noellipticity to describe the lensing of both the quasar point source and its host galaxy. This has theadvantage of producing a lens equation that can be solved analytically and lends itself to simpleand transparent interpretations.The aim of the present article is to shed a new light on the results obtained by P18 byevaluating uncertainties attached to their main results. To this end we use a different method ofdata reduction, resulting in a better angular resolution but an accordingly higher noise level; we T.T.Thai, P. Tuan-Anh & P. Darriulat work with the clean image in the plane of the sky rather than in the uv plane as P18 do; we usea simpler description of the lensing mechanism based exclusively on the analysis of the quasarpoint source images. Altogether, this simpler approach has the advantage of transparency and oflending itself easily to interpretation. It has no pretention for being better than the approach usedby P18; on the contrary, working in the uv plane allows for a more reliable treatment of noise thanworking on the clean sky plane image. But we show that it gives as proper a description of the mainresults, which justifies its use. Details of the analysis can be found in Thai T. T. 2020 [22]; here,the emphasis is on presenting and discussing the main results. Sections II and III describe the lensmodel and present its results both for the quasar point source (Section II ) and for the molecular gasemission (Section III ). They are compared with the results obtained by other authors and someissues, such as the possible presence of a companion of the quasar host galaxy, are addressedin this context. Section IV compares the data with the prediction of a simple rotating thin discmodel and offers a detailed discussion of the relative contributions of turbulence and rotation tothe observed Doppler velocity spectra. A summary and conclusions are presented in Section V . Fig. 2.
Images of RX J1131. Left: ALMA 2 mm continuum (P18). Centre: HST visible(CASTLES). Right: CHANDRA, 0.3 to 8 keV X-rays [3]. II . The quasar point source: a simple lens modelII.1 . Lens equation
The knowledge of the position and luminosity of the point images of the quasar, togetherwith the direct observation of point-like emission from the lens, allow for an accurate evaluationof the effective lensing potential that describes the lens in the vicinity of the quasar. The HSTimage positions measured with respect to the lens G are measured with a mean precision of ∼ ψ proportional to the integral of the gravity potential along theline of sight between source and observer. Convenient forms include an elliptical lens and/or anexternal shear [16, 17, 23]. Claeskens et al. 2006 [13], using such potentials, found that both giveexcellent results, the position angles of the minor axis of the ellipse and of the external shear beingidentical, ∼ ◦ east of north. This shows that introducing a shear or an ellipticity is an ad hocway to model the anisotropy of the lens mass distribution, mostly due to the presence of a massivecluster of galaxies distant by a few arcminutes in the north-eastern direction [24]. Accordingly, orpho-kinematics of quasar host galaxy RXJ 1131 5 we choose to use an effective potential of the form ψ = r r + γ r cos 2 ( ϕ − ϕ ) . (1)where r and ϕ (measured counter-clockwise from west) are polar coordinates of the point imagewith origin at the centre of the lens galaxy. The first term describes a spherical main lens ofEinstein radius r . The second term represents a shear of strength γ at position angle ϕ . Writingthat the gradient of the potential cancels, and calling ( r s , ϕ s ) the polar coordinates of the pointsource, one obtains the lens equation: r s e i ϕ s = re i ϕ ( − r − ∂ψ∂ r − ir − ∂ψ∂ϕ ) . For the potential ofRelation (1), writing separately the real and imaginary parts, one obtains: r s cos ( ϕ s − ϕ ) = r ( − γ cos 2 ( ϕ − ϕ )) − r = Ar s sin ( ϕ s − ϕ ) = r γ sin 2 ( ϕ − ϕ ) = B . (2)Relations (2) can be used to simply obtain the position of a point source from that of a point image: r s = ( A + B ) and ϕ s = ϕ + tan − ( B / A ) . (3)The first of these relations can be rewritten as r s = ( r − r ) + r γ − r γ ( r − r ) cos 2 ( ϕ − ϕ ) implying that r s , r γ and | r − r | form a triangle with an angle 2 ( ϕ − ϕ ) facing r s . Imag-ing a point source is done by eliminating r from Relations (2). One then obtains an equa-tion in ϕ giving four images when the source is inside the caustic and two when it is outside.The image magnification is obtained by differentiating the lens equation. Fig. 3. a) Best fit results comparing observations (blue) and model (red). The semi-majorand minor axes of the caustic are 589 and 446 mas respectively; those of the critical curveare 2.14 and 1.62 arcsec. b) Observed brightness distribution (mJy beam − ) in the imageplane (black for P18 and red for our data). c) Correlation between P18 (abscissa) and our(ordinate) brightness measurements. d) Doppler velocity spectra for P18 (black) and our(red) data after application of a 0.45 mJy cut on large pixels (250 ×
250 mas ). II.2 . Results
Seven parameters are adjusted by optimizing the match between the observed quadrupleHST point images and the prediction of the model: three parameters ( r , γ , ϕ ) define the lensingpotential; two account for a possible offset of the lens centre ( ∆ x , ∆ y ) with respect to the emissionof the lens galaxy (G); and two ( r s , ϕ s ) locate the point source with respect to the lens centre: its T.T.Thai, P. Tuan-Anh & P. Darriulat coordinates with respect to G are therefore x s + ∆ x and y s + ∆ y with x s = r s cos ϕ s and y s = r s sin ϕ s . Weminimize the value of the root mean square deviation δ between model and observation. The bestfit gives δ =14 mas and is illustrated in Figure 3a. The best fit values of the model parameters are r =1.84 ± γ =0.138 ± ϕ =106 ◦ ± ◦ , x s = − ± y s = − ± ∆ x = − ±
167 mas and ∆ y = − ±
20 mas in excellent agreement with the results obtainedby Claeskens et al. 2006 [13]. The quoted uncertainties are arbitrarily defined as doubling thevalue of δ . We find strong correlations between the model parameters due to the fact that whatis measured accurately is the relative position of the source with respect to the cusp of the caustic,not with respect to its centre. Magnifications cannot be reliably calculated because of the effect ofmicrolensing [15, 25] and are not used here.In summary a lens model using the effective potential of the form given by Relation (1)gives an excellent description of the astrometry of the HST images. Agreement has been obtainedwith the results of earlier analyses and a good understanding of the uncertainties attached to themodel parameters and of their correlations has been reached. However, these results probe onlythe environment of the eastern cusp of the caustic and the validity of the model at larger distances,in the region covered by the emission of the host galaxy, cannot be taken as granted. III . The host galaxy: de-lensing the observed emission of the CO(2-1) lineIII.1 . Data reduction
We use ALMA observations, project number 2013.1.01207.S (PI: Paraficz Danuta), col-lected on July 19th 2015 using the normal 12-m array with 37 antennas covering baselines be-tween 27.5 m and 1.6 km. Details of the observations and of the data reduction are given inThai T. T. 2020 [22]. Imaging was performed using the standard CLEAN algorithm applied tothe calibrated visibilities. With the aim to understand the effect of a different data reduction thanused by P18 on the results obtained, we used robust weighting rather than natural weighting asadopted by P 18. This means a better angular resolution (the beam is 70% in area compared withthe P18 beam) but a larger noise level. This dictated our choice of 0.5 as robust parameter, areasonable compromise between angular resolution and noise. We recall that a rigorous treatmentof the noise is only possible in the uv plane. Here, as we work in the image plane, we can onlyobtain an approximate estimate of its level. However, we have been careful to take this caveatin due account whenever relevant to the argument being made. Continuum emission is found toagree precisely with the results obtained by P18 and is not discussed in the present article. Weimaged the CO(2-1) data in the form of a cube of 640 ×
640 pixels, each 70 ×
70 mas , covering asquare of ± − G − and of 121Doppler velocity bins, 8.417 km s − each, covering an interval of ±
509 km s − , centred on thered-shifted ( z =0.654) frequency of the CO(2-1) line emission. The beam size is 380 ×
290 mas with position angle of 66 ◦ east of north; the noise rms level is 0.38 mJy beam − per channel. Inthe remaining of the article, we use coordinates centred at the best-fit lens centre, 60 mas southand 50 mas east of G, with the y axis pointing 16 ◦ east of north and the x axis pointing 16 ◦ northof west, perpendicularly to the external shear (namely ϕ =90 ◦ in Relations 1 to 3). In this newframe, using the axes of the caustic and critical curve as axes of coordinates, the quasar is locatedat x s = − y s = − orpho-kinematics of quasar host galaxy RXJ 1131 7 III.2 . Observed emission of the CO(2-1) line
Comparing the CO(2-1) data reduced above with the P18 data reveals the differences inbeam size (380 ×
290 mas instead of 440 ×
360 mas ) and noise level (3.0 instead of 1.9 mJyarcsec − ). The comparison is made in a square of 6.25 arcsec side centred on the lens centre,containing 125 ×
125 pixels, each 50 ×
50 mas . Eight different data sets are considered separately,each covering an 84.17 km s − interval of Doppler velocity. While the brightness distributionsof the two sets are similar when expressed in mJy beam − (Figure 3b), they are scaled relativeto each other by the ratio of the beam area when expressed in mJy pixel − , mJy arcsec − . Thecorrelation between the two sets of brightness measurements is illustrated in Figure 3c. Applyinga common cut of 10 µ Jy per pixel (4 mJy arcsec − ), which suppresses much of the noise, givesa good agreement between the brightness measurements of the two data sets: on average, theasymmetries (difference divided by sum) between the brightness integrated over each of the 8velocity intervals cancel and have a root mean square deviation of 5%. The effect of the differentnoise levels is attenuated when using larger pixels: Figure 3d compares the Doppler velocityspectra obtained by applying a cut of 0.45 mJy per pixel of 250 ×
250 mas (7.2 mJy arcsec − ).The differences between P18 and our brightness measurements give a measure of the uncertaintiesresulting from differences in data reduction. III.3 . De-lensing
We reconstruct the CO(2-1) emission in the source plane using the simple lens model de-fined in Section
II.2 with r =1.84 arcsec and γ =0.138. We consider separately 8 Doppler velocityintervals, each 84.17 km s − wide, covering between −
340 km s − and +333 km s − . We use twodifferent methods to de-lens the observed images. One is the simple direct de-lensing describedby Relations (3) in Section II.1 , which has the advantage of simplicity and of inviting transparentinterpretations. The method has two drawbacks. It lacks control of the effect of beam convolution,the de-lensed source brightness being smeared by the beam in a way that depends on the locationof the image pixel. And it introduces de-lensing noise, implying the application of a strong cut onimage brightness in order to stand aside from it. We have been careful to ensure that our resultswere robust with respect to both. For this reason one usually prefers to start from a model of thesource brightness, image it, convolve it with the beam and compare the result with the observedimage (as we do with the second method). In practice, application of the first method requiresprobing each image pixel over its whole area of 50 ×
50 mas : we use 1000 random points perpixel containing brightness in excess of 6 µ Jy pixel − and take a proper weighted average in eachsource pixel of the de-lensed values obtained for the brightness. Taking such a proper weightedaverage is not trivial. We need to know, for each image pixel, if it is imaged by a lens producing2 images (outside the caustic) or 4 images (inside the caustic). In the first case a weighted aver-age of the de-lensed brightness gives the source brightness outside the caustic and in the secondcase another weighted average of the de-lensed brightness gives the source brightness inside thecaustic. These two weighted averages need to be evaluated separately. Results are displayed inFigure 4. The second method proceeds in the opposite direction, from source to image. Imagingis done using a matrix of elements f i jkl equal to the brightness obtained in image pixel ( k , l ) bylensing source pixel ( i , j ) of unit brightness ( f i jkl is a pure number). The matrix has been calcu-lated once for all and includes the effect of beam convolution. We use 25 ×
25 source pixels of
T.T.Thai, P. Tuan-Anh & P. Darriulat ×
120 mas each, making a square covering 3 × . Optimization is made by minimiz-ing the value of < δ O > , the mean square deviation between observed and modelled brightnessin the image pixels. We use large pixels in the image plane, 250 ×
250 mas . We simply loop overthe source pixels containing brightness in excess of a threshold of ∼ − , vary theirbrightness by a small quantity ± ǫ and calculate the new value of < δ O > . The iteration uses theP18 source brightness as a first approximation. We retain as new value of the brightness the valuegiving the smaller value of < δ O > . We repeat the procedure until all pixels contain a brightness S giving a smaller value of < δ O > than for S ± ǫ . Convergence is achieved after 60 to 140 iterationsdepending on the velocity interval. Fig. 4.
Upper row: maps of the intensity, integrated over the whole velocity range, ofthe source emission as obtained by P18 (left), by direct de-lensing (centre) and by χ minimization (right). The location of the quasar is indicated by a cross. Units are Jykm s − arcsec − . Ellipses show the projection of the model disc. Middle row: projectionsof the source brightness (Jy km s − arcsec − ) on the major axis of its elliptical projectionon the sky plane ( x ′ axis in Figure 6); only pixels located inside this ellipse are included.Lower row: maps of the mean Doppler velocity. Figure 4 maps the source brightness integrated over the whole velocity range and the meanDoppler velocity, together with those obtained by P18 and by direct de-lensing. The three intensitymaps of the source emission are consistent with the elliptical projection on the sky plane of a thincircular disc centred on the quasar ( x = − y = − orpho-kinematics of quasar host galaxy RXJ 1131 9 of the ellipse of projected emission is oriented some 14 ◦ north of the x axis, meaning 30 ◦ north ofwest (L17 quote a value of 31 ◦ ). We evaluate the lengths of the major and minor axes to be ∼ ∼ ∼
19 and ∼
11 kpc respectively) corresponding to an inclination of the disc withrespect to the plane of the sky of cos − (1.6/2.7)=54 ◦ as obtained by P18. We note however that thelong axis of the P18 caustic is only ∼ ◦ south of east instead of 16 ◦ in our simple lens model. Alsoshown in Figure 4 are projections of the intensity on the major axis of the ellipse. All three distribu-tions, rather than peaking at the quasar location, display a small depletion in its vicinity. The mapsof the mean Doppler velocity display a strong gradient along the major axis, as expected from arotating thin disc. L17 have claimed evidence for the presence of a companion of the quasar hostgalaxy in the red-most velocity interval covering from 249 to 333 km s − . This companion is sup-posed to match approximately image F of Brewer & Lewis 2008 [26], some 100 to 200 mas south-west of the western cusp of the caustic. We find indeed an enhancement of emission in the red-mostvelocity interval, centred at x ∼ y ∼ ∼ ∼ Fig. 5.
Maps of the brightness integrated over the red-most velocity interval (249 to 333km s − ) are shown in the three leftmost panels for the source and in the two rightmostpanels for the image. The source maps are, from left to right, for P18, for direct de-lensing and for our best-fit result. The image maps show, again from left to right, theobserved images and the best-fit results. Colour scales are in Jy km s − arcsec − . IV . A simple rotating disc modelIV.1 . Geometry
The preceding sections have shown that the general morpho-kinematics of the de-lensedimages is a robust result of the analysis of P18. Using a simpler lens does not strongly affect the brightness distribution in the source plane and preserves the Doppler velocity distribution,typical of a thin rotating disc inclined with respect to the plane of the sky. In order to have areference with which one can make quantitative comparisons, it is instructive to construct theimage produced by a rotating disc of uniform brightness having morpho-kinematics matching thedistributions displayed in Figure 4. The geometry is illustrated in the left panel of Figure 6.Defining the position of a point in the disc by its polar coordinates ( R , θ ) with θ =0 along x ′ ,intersection of the disc with the plane of the sky containing the quasar, we see from Figure 6 that x ′ = R cos θ y ′ = R sin θ cos 54 ◦ z ′ = R sin θ sin 54 ◦ V x = − V ( R ) sin θ V y = V ( R ) cos θ cos 54 ◦ V z = V ( R ) cos θ sin 54 ◦ . (4)In our system of coordinates, with x pointing 16 ◦ north of west, the x ′ axis points to a direction of Fig. 6.
Left: geometry in a system of coordinates ( x ′ , y ′ , z ′ ) centred on the quasar with x ′ along the trace of the disc on the plane of the sky and z ′ perpendicular to the plane of thesky. Centre and right: the image of a disc of uniform brightness lensed by the simple lenspotential (centre) is compared with observation (right). Units are arbitrary for the modeland Jy km s − arcsec − with a cut at 0.67 Jy km s − arcsec − for the observed data. ◦ − ◦ =14 ◦ ; centring the disc on the quasar of coordinates ( x , y )=( − − x + . = x ′ cos 14 ◦ − y ′ sin 14 ◦ y + . = x ′ sin 14 ◦ + y ′ cos 14 ◦ . Using the above equations and the results obtained in the preceding section, we image a disc ofuniform brightness centred on the quasar, inclined by 54 ◦ with respect to the plane of the skyand having a radius R disc =1.35 arcsec. The inclination of the disc simply divides the disc planebrightness by cos54 ◦ to give the projected brightness. The images are convolved with the beam:the central panel of Figure 6 shows that the morphology of the image brightness produced by auniform disc is confined within an approximately elliptical annular band centred on the lens/quasarregion. We define two ellipses, E+ and E − , delimiting the outer and respectively inner edges of theuniform disc image as shown in Figure 6. The right panel of Figure 6 displays the image brightnessof the observed data. It is well contained within ellipses E+ and E − and populates the middle ofthe band delimited by them. In order to use coordinates adapted to this morphology, we define a orpho-kinematics of quasar host galaxy RXJ 1131 11 parameter λ such that λ = − − and +0.5 on E+, namely λ = ± ± . Precisely, a pointon the sky plane having Cartesian coordinates ( x = r cos ω , y = r sin ω ) and polar coordinates( r , ω ) is imaged as ( λ , ω ) with λ = [ r − ( r + + r − ) / ] / ( r + − r − ) ; here, r + and r − are the points ofposition angle ω on ellipses E+ and E − respectively. In the remaining of the section we work inthis new system of coordinates, ( λ , ω , V z ). IV.2 . The data cube
We construct a data cube in ( λ , ω , V z ) coordinates to be compared with model predictions.We set its limits as − < λ < +0.5, 0 < ω < ◦ and − < V z <
333 km s − ; we segment it in 20 λ bins, each 0.05 wide, in 18 ω bins, each 20 ◦ wide and in 16 V z bins, each ∼
42 km s − wide. Weconstruct the new data cube starting from the reduced data presented in Section III , distributedin 125 ×
125 pixels in the sky plane, each 50 ×
50 mas in area, and 80 Doppler velocity bins,each 8.417 km s − wide. Each pixel is probed in 100 points randomly distributed over its area.Working in the sky plane rather than in the uv plane prevents in practice a rigorous treatmentof the noise. Instead, we compare in Figure 7 the distributions obtained by applying no brightnesscut on the original data cube elements to those obtained by applying a cut ( ∼ σ ) such that theflux integrated over the whole data cube is the same [22]. The upper row of Figure 7 comparesthe projections of the data cube on each of the coordinates with and without application of the cut.We repeat the exercise with the data reduced by P18. Applying a cut (again ∼ σ ) producingthe same integrated flux as when no cut is applied gives the results shown in the second rowof Figure 7. Finally, the third row compares the present data with the P18 data, both with nobrightness cut being applied.In summary, the differences in noise level and angular resolution between the analysesof the present work and of P18, apart from a global rescaling factor ( ∼ ∼
13 mJy per histogram bin.
IV.3 . The disc model
We parameterize the rotation curve as V ( R ) = V ( e R / R ∗ − )( e R / R ∗ + ) and the disc brightness is takenuniform over a disc of mean radius R disc smeared radially by a Gaussian of dispersion σ disc .Thischoice has been made after having explored results obtained using different possible forms. Whilebeam convolution is properly accounted for by the model, no noise contribution is included. Wedecided to do so after having produced fits accounting for a Gaussian noise contribution mimickingthat present in the observed data cube. The results were essentially unaffected. Optimization ofthe values of the four parameters is made by minimizing a χ function describing the quality of thematch between model and observations. Rather than evaluating the value of χ as a sum over the20 × × χ is evaluated using as uncertainty a commonarbitrary value of 10 mJy per histogram bin and is divided by the number of degrees of freedom.The best fit is obtained for V =405 km s − , R ∗ =0.22 arcsec (1.6 kpc), R disc =1.10 arcsec (7.7 kpc)and σ disc =0.32 arcsec (2.2 kpc). It corresponds to a steeper rise of the rotation curve at the originthan implied by P18 and L17. The results are illustrated in Figure 8. Fig. 7.
Projections of the new data cube (see text) on λ (left), ω (centre) and V z (right).The upper row compares the present data without cut (black) or with a 1.6- σ cut (red)corresponding to 11.7 µ Jy per original data cube element (50 × × km s − ).The second row compares the P18 data without cut (black) or with a 1.5- σ cut (red)corresponding to 7.1 µ Jy per original data cube element. The third row compares thepresent data (black) with the P18 data (red).
The best fit value of χ is ∼
3, meaning that the mean deviation between model and of ob-servation is ∼
17 mJy per histogram bin. As could be expected, the model is unable to account forthe observed brightness inhomogeneity, particularly well revealed by the central panel of Figure 8,which we discuss in the next section. The results obtained for the rotation curve and for the discbrightness are essentially uncorrelated; but the results obtained for each independently display avery strong correlation between V and R ∗ on the one hand and between R disc and σ disc on theother. This is illustrated in Figure 9 which maps the value of χ in the R ∗ vs V and σ disc vs R disc planes respectively. IV.4 . Brightness inhomogeneity
Maps of the observed and modelled brightness are compared in Figure 10 in the λ vs ω , λ vs V z and V z vs ω planes. In general, the agreement between observed and modelled mapsis remarkable given the crudeness of the model. The V z vs ω maps display an oscillation thatreveals very clearly the disc rotation. The strong asymmetry of the Doppler velocity spectrum,present in both model and observed brightness distributions, is the direct result of the lens caustic orpho-kinematics of quasar host galaxy RXJ 1131 13 Fig. 8.
Comparison between observations (black) and best fit model (red and blue) areshown as projections of the data cube on λ (left), on ω (centre) and on V z (right). The red(blue) histograms are with (without) allowance for the presence of an emission hot spot(see Section IV.4 ). Fig. 9.
Left: dependence of χ on the model parameters: V (abscissa) and R ∗ (ordinate)in the left-most panels and R disc (abscissa) and σ disc (ordinate) in the right-most panels.In each pair of panels, the results obtained without and with allowance for an excess ofemission in the vicinity of the quasar are ordered from left to right. In each panel the othermodel parameters are set at their best-fit values. The values of F E listed in the inserts referto the excess of emission discussed in Section IV.4 . being entirely located in the red-shifted region of the rotating disc. These maps provide impor-tant information on the morpho-kinematics without requiring de-lensing of the observed images.Differences between model and observation are illustrated in the lower panels of the figure. Theyare dominantly in the two red-most V z bins and confined to a narrow interval of the polar angle ω ,approximately between 200 ◦ and 240 ◦ , where they take the form of an excess for small negativevalues of λ and a depletion for small positive values. We define two rectangles in the λ vs ω plane that delimit these regions: an excess E (200 ◦ < ω < ◦ , − < λ < − ◦ < ω < ◦ ,0.15 < λ < ∼ ∼
5. The companion image of D is small and centred at ( x , y ) ∼ (0.75, 0.25) arcsec; its magnificationspans values between 0.5 and 1. In contrast, the source of E is located inside the caustic and Fig. 10.
Maps of the brightness projected on the ( λ vs ω ), ( λ vs V z ) and ( V z vs ω )planes (respectively from left to right) and integrated over the third coordinate, V z , ω and λ respectively. Bins are 0.05 wide in λ , 20 ◦ wide in ω and 42 km s − wide in V z .Upper panels display observed images with no brightness cut. Central panels display theimages obtained from lensing the best fit disc model. Contours of the disc model (centralrow) are superimposed on observations (upper row). Lower panels map the differencebetween observations and model. The rectangles in the left panel show the regions ofstrong inhomogeneity discussed in Section 4.4, the critical curve (black) is also plotted. touches its edge; its magnification spans accordingly a very broad range of values starting at ∼ ∼
20. It produces three additional images, whichspan magnifications ranging between 1.5 and 5.5; one of these, centred at ( x , y ) ∼ ( −
1, 2) arcsecis in a region where the lower-left panel of Figure 10 shows that the observed brightness exceedsslightly that of the model. The fact that de-lensing E produces a source limited by the causticcurve is not an artefact of the lensing mechanism. In fact, the source that produces E overlaps thecaustic but only the part of it that is inside contributes to image E. The part that is outside producesinstead additional images, one located near ( x , y ) ∼ ( −
1, 2) and the other close to E.Interpreting the inhomogeneity as the conjunction of an excess and a depletion is arbitrary;it assumes that on average observed and modelled fluxes are equal, which they have no reasonto be; it is indeed simpler and more natural to blame the inhomogeneity on a single hot spot, inwhich case the flux predicted by the model needs to be scaled down, or as a single depletion, inwhich case the flux predicted by the model needs to be scaled up. We tried both interpretations by orpho-kinematics of quasar host galaxy RXJ 1131 15 including in the model an enhancement (depletion) of emission by a factor F E ( F D ) in the sourceregion of E (D) shown in Figure 11. In both cases the best fit values of the model parameters, V , R ∗ , R disc and σ disc , were essentially unaffected but allowing for a depletion did not producelower values of χ ; in contrast, as illustrated in the right panels of Figure 11, allowing for anexcess resulted in a significant improvement of the quality of the fit. The best fit values of themodel parameters are now V =435 km s − , R ∗ =0.26 arcsec (1.8 kpc), R disc =1.10 arcsec (7.7 kpc), σ disc =0.32 arcsec (2.25 kpc) and F E =2.5 (see Figure 9).The above results call for a word of caution. Qualitatively, the presence of an excess ofemission in the vicinity of the quasar is likely to be a robust result of the analysis; it had been notedby P18 who suggest that it is most likely associated with a site of on-going star-formation. How-ever, its location on top of the caustic makes it difficult to specify reliably its morpho-kinematicalproperties. The lens model used in the present work is very crude and so is also that used byP18: none of these allows for local inhomogeneity of the lensing potential, they describe it in verygeneral terms; the vicinity of the caustic is particularly difficult to model reliably, small variationsof the lensing potential having strongly amplified effects in the construction of the images. As weshall see in the next section, the overall remarkable agreement between model and observationsdisplayed in Figure 10 does not sustain a significantly more detailed scrutiny. IV.5 . Rotation and turbulence
The evaluation of the rotation curve made by P18 and L17 consists in defining a band brack-eting the major axis of the projection of the disc on the plane of the sky ( x ′ axis in Figure 6). Wechoose for it a width of ± ∼
19 kpc), divided in nine segments,each having a length of 0.3 arcsec ( ∼ θ of the sine of the inclination angle in Relation (4). However, includingwarping in the model by writing V z = V ( R ) cos θ sin ϕ with ϕ depending simply on θ and R , givesonly a modest improvement of the match between model and observations. This suggests that amore complex dynamics than described by the simple model is at stake.Beam convolution causes an effective smearing of the disc region contributing to each of thenine segments, making the radial distribution of the mean Doppler velocity measured in the centralsegments less steep and causing an increase of the velocity dispersion. This important result is awarning: measuring the velocity dispersion requires a careful evaluation of the contribution ofrotation within the finite angular resolution of relevance. Making the angular acceptance smallerwill decrease the direct contribution of rotation, independent from beam convolution, but willnot decrease the contribution resulting from the smearing caused by the beam. This contributionis important: using the central segment as an example, the line width predicted by the modelincreases from 60 km s − to nearly 100 km s − when accounting for beam convolution, namely abeam contribution of over 70 km s − (both contributions add up in quadrature). Fig. 11.
Left panels display the brightness distribution in the source plane obtained byde-lensing depletion D (left) and excess E (centre-left) respectively. Right panels displaythe dependence of χ on the factors F D and F E measuring the amplitudes of the depletion(centre-right) and excess (right) respectively when the other model parameters are fixedat their best-fit values. Fig. 12.
Dependence of < V z > on x ′ (major axis of the projection of the disc on the planeof the sky) for the data (black) and the model (red) is shown in the left panel and therotation curve in the right panel together with the L17 (blue curve) and P18 (red crosses)results. The differences between observed and predicted velocity dispersions may receive contribu-tions from the intrinsic line width (turbulence), which is absent from the model, or from differentrotation contributions, as could be caused by disc warping. We recall that P18 claim that the obser-vation of a high dispersion in the vicinity of the quasar demonstrates that the region of enhancedemission studied in the preceding section is the seat of increased gas turbulence. However, thedifferences observed in the central segments of the band used to evaluate the rotation curve do notconsist of a simple broadening of the line, as would be expected from turbulence, but reveal ratheran excess of red-shifted emission.Maps of the intensity, mean Doppler velocity and Doppler velocity dispersion in the λ vs ω plane reveal important differences between model and observations [22]: predicted mean Dopplervelocities tend to be too small in the 90 ◦ < θ < ◦ quadrant and too large in the 270 ◦ < θ < ◦ quadrant of the disc plane; velocity dispersions reach some 120 km s − for the model and some orpho-kinematics of quasar host galaxy RXJ 1131 17 Fig. 13.
Line profiles in the region | λ | < ◦ < ω < ◦ . Left: observationswith cuts at 11.7 µ Jy (black) or 16 µ Jy (red); the lines show Gaussian fits with σ =85km s − and 56 km s − respectively. Gaussian fits to the predictions of the disc modelwithout (centre) and with (right) beam convolution give σ values of 18 km s − and 58km s − respectively.
160 km s − for the observations. We failed to obtain a substantial improvement of the quality ofthe match between model and observations by including a simple description of disc warping inthe model: this confirms the need for a more complex dynamics than implied by the simple discmodel. From a number of comparisons between the line widths predicted by the model with andwithout beam convolution, the contribution is found to be in the range between 50 and 70 km s − .An independent evaluation is obtained by selecting the central region of the λ vs ω plane(left column of Figure 10) defined as | λ | < | ω − ◦ | < ◦ . We segment it in 10 × λ and 10 ◦ in ω . For each pixel ( i , j ) we evaluate the mean value < V z > i j of theDoppler velocity and show the distribution of V z − < V z > i j in Figure 13. We obtain this way anevaluation of the line profile. The observed line widths ( σ ) are 85 and 56 km s − for cuts at 11.7(1.5- σ ) and 16 µ Jy (2- σ ) per 50 × × km s − respectively. As the model includes nointrinsic line broadening, the line broadening shown by the model is entirely due to the contributionof rotation. Gaussian fits give σ values of 18 and 58 km s − without and with beam convolutionrespectively, implying a contribution of ∼
50 km s − from the smearing effect of the beam size.The line width measured with the 16 µ Jy cut does not require any contribution of turbulencewhile that measured with an 11.7 µ Jy cut allows for a contribution reaching ∼
60 km s − : withouta precise understanding of the contribution of noise, it is therefore difficult to obtain a reliableevaluation of the line broadening caused by turbulence.In summary, a simple rotating disc model gives a global picture that describes well thegeneral trend of the observed kinematics. However, such a model fails to reproduce quantitativelythe details of the Doppler velocity distribution: it reveals a more complex dynamics. Attempts todescribe the mismatch as the result of a simple disc warping have failed. An important result ofour analysis is that the smearing in the image plane caused by the beam size contributes between50 and 70 km s − to the dispersion of the Doppler velocity. It combines in quadrature with thecontribution of rotation within the angular acceptance being probed by the pixel size of relevance.Taken together, these effects make it difficult to evaluate reliably the contribution of turbulence tothe line width. But they prevent claiming that such contribution is important: within experimental uncertainties the line width can be accounted for by the dispersion of rotation velocities within thebeam size. V . Summary and conclusions
The present study of the emission of the CO(2-1) molecular line by the host galaxy of quasarRX J1131 uses ALMA observations of unprecedented quality; they have been previously analysedin much detail by their proponents, [18], P18. We have shown that the HST images of the quasarand of the lens galaxy are very well reproduced by a simple lensing potential, sphere+externalshear. We obtained parameters that are in agreement with the results of previous authors withinerrors. We discussed the uncertainties attached to the model parameters and their correlations andexplained them by remarking that these results probe only the environment of the eastern cusp ofthe caustic. At this stage, the validity of the model at larger distances, in the region covered by theemission of the host galaxy, cannot be taken as granted. But we demonstrated its validity from ouranalysis of the ALMA observations of the CO(2-1) line emission. We reduced the raw data andproduced cleaned images, at variance with P18 who perform their analysis in the uv plane. Weobtained morpho-kinematics properties of the gas that are in good agreement with the results ofP18. We used two different methods for de-lensing these images, which produce source brightnessdistributions in agreement with each other as well as with the P18 results. We have shown that thegeneral morphology of the image is dictated by the properties of the lens and confined within aband bracketing the critical curve. This leaves in practice no freedom to conceive lensing potentialsthat would be significantly different from those that we and other authors have been using. Itguarantees the robustness of the results obtained by P18 within observational uncertainties.We used polar coordinates better adapted to this geometry than Cartesians to reveal veryclearly the evidence for rotation. We looked for possible deviations from the predictions of asimple rotating disc model, in particular for the possible presence of a companion as was predictedby L17 using Plateau de Bure observations of much less good angular resolution than the ALMAobservations; but our analysis did not support such a presence. We paid special attention to the red-most velocity interval, which displays a particularly complex morphology. We found evidence forenhanced emission revealing a significant deviation from the prediction of the simple disc model.However the associated hot spot in the source plane overlaps the caustic, implying importantuncertainties on its de-lensed properties.We discussed the rotation curve in some detail and argued that it probably rises faster thanassumed by P18 in the vicinity of the quasar; we noted that the angular resolution is such thatthe dispersion of the rotation velocity within the beam size in the image plane causes an importanteffective broadening of the line width, which we evaluated in the range between 50 and 70 km s − .This is sufficient to account for the observed line width when applying a ∼ σ cut on the data.However, uncertainties attached to the effect of noise prevent from giving a reliable evaluation ofthe contribution of turbulence. This is at variance with the conclusion of P18 who attribute the highDoppler velocity dispersion observed in the vicinity of the AGN to gas turbulence exclusively.P18 obtained from their analysis a number of results concerning the physical properties ofthe galaxy, which they compared with those of other similar galaxies. The nature of our studyprevents us from adding much to their conclusions. We simply remark that their estimate ofthe total dynamical mass enclosed within 5 kpc, (1.46 ± × M ⊙ , may be affected by thesteepness of the rotation curve: at fixed radius, the dynamical mass is proportional to the square of orpho-kinematics of quasar host galaxy RXJ 1131 19 the rotation velocity. The analysis presented in the preceding sections implies therefore a possiblyhigher value of the gas mass, obtained by P18 and accordingly affect the value of the CO − H conversion factor. ACKNOWLEDGMENT
We express our deep gratitude to Professors Frederic Courbin and Matus Rybak who kindlyprovided us with copies of data files summarizing the results of the P18 analysis. This paper usesarchival ALMA data from project 2013.1.01207.S (PI: Paraficz Danuta). ALMA is a partnershipof ESO (representing its member states), NSF (USA), NINS (Japan), NRC(Canada), NSC/ASIAA(Taiwan), and KASI (South Korea), in cooperation with Chile. The Joint ALMA Observatory isoperated by ESO, AUI/NRAO and NAOJ. We are deeply indebted to the ALMA partnership,whose open access policy means invaluable support and encouragement for Vietnamese astro-physics. Financial support from the World Laboratory and VNSC is gratefully acknowledged.This research is funded by the Vietnam National Foundation for Science and Technology Devel-opment (NAFOSTED) under grant number 103.99-2019.368.
REFERENCES [1] Adam R. et al., 2016, A&A 596, A108[2] Ade P.A.R. et al., 2016, A&A 594, A13[3] Reis R.C., Reynolds M.T., Miller J.M. & Walton D.J., 2014, Nature, 507, 207-209[4] Tewes, M., Courbin F., Meylan G. et al., 2013, A&A, 556, A22.[5] Bonvin, V., Courbin, F., Suyu, S.H. et al., 2017, MNRAS, 465, 4914[6] Chartas G., Kochanek C.S. Dai, X. et al., 2012, ApJ, 757, 137[7] Sugai H., Kawai A., Shimono A. et al., 2007, ApJ, 660, 2[8] Anh P.T., Boone F., Hoai D.T. et al., 2013, A&A 552 L12[9] Anh P.T. 2014, PhD thesis, Vietnam Institute of Physics and Toulouse University[10] Stacey H. R., McKean J. P., Robertson, N. C., et al., 2018, MNRAS, 476, 5075[11] Leung T. K. D., Riechers D. A., & Pavesi, R., 2017, ApJ, 836, 180[12] CASTLES http://cfa.harvard.edu/castles/Individual/RXJ1131.html[13] Claeskens, J.-F., Sluse, D., Riaud, P., & Surdej, J., 2006, A&A, 451, 865[14] Birrer, S., Amara, A., & Refregier, A., 2016, Cosmology Astropart. Phys., 8, 020[15] Chen, G.C.-F., Suyu, S. H., Wong, K. C., et al., SHARP - III, 2016, MNRAS, 462, 3457[16] Blandford, R.D., Kochanek, C.S., Kovner, I. and Narayan, R., 1989, Science 245/4920, 824[17] Saha P. & Williams, L.L.R., 2003, A.J., 125/6, 2769[18] Paraficz, D., Rybak, M., McKean, J.P. et al., 2018, A&A, 613, A34[19] Rybak, M., McKean, J. P., Vegetti, S., Andreani, P., & White, S. D. M. 2015a, MNRAS, 451, L40[20] Rybak M., Vegetti S., McKean J. P., Andreani P. & White S. D. M. 2015b, MNRAS, 453, L26[21] Vegetti S., & Koopmans L. V. E. 2009, MNRAS, 392, 945[22] Thai T.T., 2020,