Morpho-kinematics of the molecular gas in a quasar host galaxy at redshift z = 0:654

Nội dung text: Morpho-kinematics of the molecular gas in a quasar host galaxy at redshift z = 0:654

1. Communications in Physics, Vol.31, No. 2 (2021), pp. 149-168 DOI:10.15625/0868-3166/15599 MORPHO-KINEMATICS OF THE MOLECULAR GAS IN A QUASAR HOST GALAXY AT REDSHIFT z = 0.654 T. T.THAI1,2, P. TUAN-ANH1,†, P. DARRIULAT1, D. T. HOAI1, P. T. NHUNG1, P. N. DIEP1, N. B. NGOC1 AND N. T. PHUONG1 1Vietnam National Space Center, Vietnam Academy of Science and Technology 18 Hoang Quoc Viet, Cau Giay, Hanoi, Vietnam 2Graduate University of Science and Technology 18 Hoang Quoc Viet, Cau Giay, Hanoi, Vietnam E-mail: †ptanh@vnsc.org.vn Received 14 October 2020 Accepted for publication 29 December 2020 Published 27 January 2021 Abstract. We present a new study of archival ALMA observations of the CO(2-1) line emission of the host galaxy of quasar RX J1131 at redshift z=0.654, lensed by a foreground galaxy. A simple lens model is shown to well reproduce the optical images obtained by the Hubble Space Telescope. Clear evidence for rotation of the gas contained in the galaxy is obtained and a simple rotating disc model is shown to give an excellent overall description of the morpho-kinematics of the source. The possible presence of a companion galaxy suggested by some previous authors is not confirmed. Detailed comparison between model and observations gives evidence for a more complex dynamics than implied by the model. Doppler velocity dispersion within the beam size in the 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. INTRODUCTION I.1. General features RX J1131-1231 (simply called RX J1131 in the following), is a distant quasar, at redshift zs ∼0.654, corresponding to a distance of ∼1.45 Gpc or a time of ∼7.5 Gyr after the Big Bang, about half way from us, and some 4 Gyr later than the time of maximal star formation [1, 2]. At such distance, 1 arcsec spans 7.03 kpc. It hosts a Super Massive Black Hole (SMBH) in its centre ©2021 Vietnam Academy of Science and Technology
2. 150 T. T. THAI et al. 8 with a mass of ∼2 10 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 zL ∼0.295. They are the object of numerous studies, in particular aiming at a better understanding of the cosmological parameters governing the expansion of the Universe ( [4, 5] and references therein). Microlensing caused by stars transiting across the line of sight to the quasar has been used to study the structure of the lens halo ( [6, 7] and references therein). Infrared observations obtained by Herschel [10] have measured the spectral energy distribu- Fig. 1. Left: dependence on the redshifts of the source (zS in abscissa) and of the lens (zL in ordinate) of the ratio between their respective angular diameter distances daS/daL [1,2]). The relative size of the lens with respect to the source is proportional to daS/daL. The stars show the locations of quasars RX J1131 (black, P18) and RX J0911 (red, [8], [9]). The sizes of the host galaxies are compared to the size of the lens caustic in the central and right panels respectively. tion (SED), and archival VLA observations (Program ID: AW741; PI: Wucknitz) analysed by Leung 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 the background quasar. Thorough analyses of high angular resolution HST optical and NIR images [12–14] and of Keck Adaptive Optics images [15] have produced a detailed description of the lensing properties in the neighbourhood of the quasar. They reveal a typical long axis quad configuration [16,17], the quasar being located within the eastern cusp of the lens caustic. As emission from the lens galaxy is 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 the quasar host galaxy covers the whole caustic curve and extends even farther out, one cannot take it as granted that the simple lens model obtained from the study of the quasar images reliably applies to the whole host galaxy. This is at variance with the gravitational lensing of quasar hosts that 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, which shows the location of RX J1131 in the plane zL vs zS, lens vs source redshifts, together with that of other multiple imaged systems [12], and compares it with the case of a typical farther away quasar host galaxy, RX J0911 [8, 9].
3. MORPHO-KINEMATICS OF THE MOLECULAR GAS IN A QUASAR HOST GALAXY 151 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 ∼0.4×0.3 arcsec2 providing unprece- dented image quality. The data have been analysed in much detail by the team who proposed the observation ( [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 by L17 at the Plateau de Bure Interferometer (PdBI) with an angular beam size (FWHM) of 4.4×2.0 arcsec2, a spectral resolution of ∼21.5 km s−1 and a noise level of ∼1.45 mJy beam−1 per chan- nel. In addition to their CO(2-1) PdBI observations, L17 used the Combined Array for Research in Millimeter-wave Astronomy (CARMA) to detect the CO(3-2) line emission but the signal to noise 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 ∼0.3 arcsec, an order of magnitude better than for the PdBI data of L17. These are the data used in the present article, about which P18 kindly sent us useful documentation that complements the published article (private communications by Professors Frederic Courbin and Matus Rybak). The continuum image obtained by P18 shows four clearly separated compact components, three coincident with the lensed optical point images and one associated with the lens galaxy (Figure 2 left). In contrast with continuum emission, CO(2-1) line emission detects no signal from the lens galaxy. The velocity integrated intensity map clearly shows (76-σ) line emission extended over a complete Einstein ring. The map of velocity dispersion displays values covering from ∼10 to ∼50 km s−1. The authors of P18 note that peaks in the intensity map and the region of high velocity dispersion are coincident and probably reveal on-going star-formation; these peaks are not strictly coincident with the quasar emission. As mentioned above, one cannot take it as granted that the simple lens model obtained from the study of the quasar images reliably applies to the whole host galaxy. For this reason P18 reconstruct 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◦ with respect to the plane of the sky, having rotation velocity reaching over 400 km s−1. The model used by P18 to describe the lens includes an external shear, ellipticity of the main lens and a contribution from the small satellite galaxy revealed by the HST image north of the main lens (central panel of Figure 2). The contribution of the latter is known to be very small [13] and the contribution of ellipticity, while different from that of an external shear, does not affect strongly the general picture: excellent descriptions of the lensing of the quasar point source are obtained with either external shear alone or ellipticity alone and their combination in the P18 model causes a rotation of only 4◦ of the caustic with respect to a model ignoring ellipticity. For this reason our approach in the present article is to use a simple lens potential with only external shear and no ellipticity to describe the lensing of both the quasar point source and its host galaxy. This has the advantage of producing a lens equation that can be solved analytically and lends itself to simple and transparent interpretations. The aim of the present article is to shed a new light on the results obtained by P18 by evaluating uncertainties attached to their main results. To this end we use a different method of data reduction, resulting in a better angular resolution but an accordingly higher noise level; we
4. 152 T. T. THAI et al. work with the clean image in the plane of the sky rather than in the uv plane as P18 do; we use a simpler description of the lensing mechanism based exclusively on the analysis of the quasar point source images. Altogether, this simpler approach has the advantage of transparency and of lending itself easily to interpretation. It has no pretention for being better than the approach used by P18; on the contrary, working in the uv plane allows for a more reliable treatment of noise than working on the clean sky plane image. But we show that it gives as proper a description of the main results, 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 lens model and present its results both for the quasar point source (Sec. II) and for the molecular gas emission (Sec. III). They are compared with the results obtained by other authors and some issues, such as the possible presence of a companion of the quasar host galaxy, are addressed in this context. Sec. IV compares the data with the prediction of a simple rotating thin disc model and offers a detailed discussion of the relative contributions of turbulence and rotation to the observed Doppler velocity spectra. A summary and conclusions are presented in Sec. 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 MODEL II.1. Lens equation The knowledge of the position and luminosity of the point images of the quasar, together with the direct observation of point-like emission from the lens, allow for an accurate evaluation of the effective lensing potential that describes the lens in the vicinity of the quasar. The HST image positions measured with respect to the lens G are measured with a mean precision of ∼8 mas. The deflection induced by the lens as a function of the sky coordinates of the image can be described by an effective potential ψ proportional to the integral of the gravity potential along the line of sight between source and observer. Convenient forms include an elliptical lens and/or an external shear [16, 17, 23]. Claeskens et al. 2006 [13], using such potentials, found that both give excellent results, the position angles of the minor axis of the ellipse and of the external shear being identical, ∼16◦ east of north. This shows that introducing a shear or an ellipticity is an ad hoc way to model the anisotropy of the lens mass distribution, mostly due to the presence of a massive cluster of galaxies distant by a few arcminutes in the north-eastern direction [24]. Accordingly,
5. MORPHO-KINEMATICS OF THE MOLECULAR GAS IN A QUASAR HOST GALAXY 153 we choose to use an effective potential of the form 1 ψ = r r + γ r2 cos2(ϕ − ϕ ). (1) 0 2 0 0 where r and ϕ (measured counter-clockwise from west) are polar coordinates of the point image with origin at the centre of the lens galaxy. The first term describes a spherical main lens of Einstein radius r0. The second term represents a shear of strength γ0 at position angle ϕ0. Writing that the gradient of the potential cancels, and calling (rs, ϕs) the polar coordinates of the point iϕs iϕ −1 ∂ψ −2 ∂ψ source, one obtains the lens equation: rse = re (1 − r ∂r − ir ∂ϕ ). For the potential of Relation (1), writing separately the real and imaginary parts, one obtains: rs cos(ϕs − ϕ) = r(1 − γ0 cos2(ϕ − ϕ0)) − r0 = A rs sin(ϕs − ϕ) = rγ0 sin2(ϕ − ϕ0) = B. (2) Relations (2) can be used to simply obtain the position of a point source from that of a point image: 2 2 1 −1 rs = (A + B ) 2 and ϕs = ϕ + tan (B/A). (3) 2 2 2 2 The first of these relations can be rewritten as rs = (r − r0) + r γ0 − 2rγ0(r − r0)cos2(ϕ − ϕ0) implying that rs, rγ0 and |r − r0| form a triangle with an angle 2(ϕ − ϕ0) facing rs. 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-major and minor axes of the caustic are 589 and 446 mas respectively; those of the critical curve are 2.14 and 1.62 arcsec. b) Observed brightness distribution (mJy beam−1) in the image plane (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 mas2). II.2. Results Seven parameters are adjusted by optimizing the match between the observed quadruple HST point images and the prediction of the model: three parameters (r0,γ0,ϕ0) define the lensing potential; two account for a possible offset of the lens centre (∆x, ∆y) with respect to the emission of the lens galaxy (G); and two (rs, ϕs) locate the point source with respect to the lens centre: its
6. 154 T. T. THAI et al. coordinates with respect to G are therefore xs+∆x and ys+∆y with xs = rs cosϕs and ys=rs sinϕs. We minimize the value of the root mean square deviation δ between model and observation. The best fit gives δ=14 mas and is illustrated in Figure 3a. The best fit values of the model parameters are ◦ ◦ r0=1.84±0.02 arcsec, γ0=0.138±0.007, ϕ0=106 ±1 , xs=−0.47±0.03 arcsec, ys=−0.14±0.01 arcsec, ∆x=−50±167 mas and ∆y=−60±20 mas in excellent agreement with the results obtained by Claeskens et al. 2006 [13]. The quoted uncertainties are arbitrarily defined as doubling the value of δ 2. We find strong correlations between the model parameters due to the fact that what is 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 of microlensing [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 obtained with the results of earlier analyses and a good understanding of the uncertainties attached to the model parameters and of their correlations has been reached. However, these results probe only the 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) LINE III.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 in Thai T. T. 2020 [22]. Imaging was performed using the standard CLEAN algorithm applied to the calibrated visibilities. With the aim to understand the effect of a different data reduction than used by P18 on the results obtained, we used robust weighting rather than natural weighting as adopted by P 18. This means a better angular resolution (the beam is 70% in area compared with the P18 beam) but a larger noise level. This dictated our choice of 0.5 as robust parameter, a reasonable compromise between angular resolution and noise. We recall that a rigorous treatment of the noise is only possible in the uv plane. Here, as we work in the image plane, we can only obtain an approximate estimate of its level. However, we have been careful to take this caveat in due account whenever relevant to the argument being made. Continuum emission is found to agree precisely with the results obtained by P18 and is not discussed in the present article. We imaged the CO(2-1) data in the form of a cube of 640×640 pixels, each 70×70 mas2, covering a square of ±22.4 arcsec centred on the continuum emission of the lens galaxy − G − and of 121 Doppler velocity bins, 8.417 km s−1 each, covering an interval of ±509 km s−1, centred on the red-shifted (z=0.654) frequency of the CO(2-1) line emission. The beam size is 380×290 mas2 with position angle of 66◦ east of north; the noise rms level is 0.38 mJy beam−1 per channel. In the remaining of the article, we use coordinates centred at the best-fit lens centre, 60 mas south and 50 mas east of G, with the y axis pointing 16◦ east of north and the x axis pointing 16◦ north ◦ of west, perpendicularly to the external shear (namely ϕ0=90 in Relations 1 to 3). In this new frame, using the axes of the caustic and critical curve as axes of coordinates, the quasar is located at xs=−0.49 arcsec and ys=−0.005 arcsec.
7. MORPHO-KINEMATICS OF THE MOLECULAR GAS IN A QUASAR HOST GALAXY 155 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 in beam size (380×290 mas2 instead of 440×360 mas2) and noise level (3.0 instead of 1.9 mJy arcsec−2). The comparison is made in a square of 6.25 arcsec side centred on the lens centre, containing 125×125 pixels, each 50×50 mas2. Eight different data sets are considered separately, each covering an 84.17 km s−1 interval of Doppler velocity. While the brightness distributions of the two sets are similar when expressed in mJy beam−1 (Figure 3b), they are scaled relative to each other by the ratio of the beam area when expressed in mJy pixel−1, mJy arcsec−2. The correlation between the two sets of brightness measurements is illustrated in Figure 3c. Applying a common cut of 10 µJy per pixel (4 mJy arcsec−2), which suppresses much of the noise, gives a good agreement between the brightness measurements of the two data sets: on average, the asymmetries (difference divided by sum) between the brightness integrated over each of the 8 velocity intervals cancel and have a root mean square deviation of 5%. The effect of the different noise levels is attenuated when using larger pixels: Figure 3d compares the Doppler velocity spectra obtained by applying a cut of 0.45 mJy per pixel of 250×250 mas2 (7.2 mJy arcsec−2). The differences between P18 and our brightness measurements give a measure of the uncertainties resulting 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 Subsec. II.2 with r0=1.84 arcsec and γ0=0.138. We consider separately 8 Doppler velocity intervals, each 84.17 km s−1 wide, covering between −340 km s−1 and +333 km s−1. We use two different methods to de-lens the observed images. One is the simple direct de-lensing described by Relations (3) in Subsec. II.1, which has the advantage of simplicity and of inviting transparent interpretations. 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 location of the image pixel. And it introduces de-lensing noise, implying the application of a strong cut on image brightness in order to stand aside from it. We have been careful to ensure that our results were robust with respect to both. For this reason one usually prefers to start from a model of the source brightness, image it, convolve it with the beam and compare the result with the observed image (as we do with the second method). In practice, application of the first method requires probing each image pixel over its whole area of 50×50 mas2: we use 1000 random points per pixel containing brightness in excess of 6 µJy pixel−1 and take a proper weighted average in each source pixel of the de-lensed values obtained for the brightness. Taking such a proper weighted average is not trivial. We need to know, for each image pixel, if it is imaged by a lens producing 2 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 second case another weighted average of the de-lensed brightness gives the source brightness inside the caustic. These two weighted averages need to be evaluated separately. Results are displayed in Figure 4. The second method proceeds in the opposite direction, from source to image. Imaging is done using a matrix of elements fi jkl equal to the brightness obtained in image pixel (k,l) by lensing source pixel (i, j) of unit brightness ( fi 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
8. 156 T. T. THAI et al. 120×120 mas2 each, making a square covering 3×3 arcsec2. Optimization is made by minimiz- ing the value of , the mean square deviation between observed and modelled brightness in the image pixels. We use large pixels in the image plane, 250×250 mas2. We simply loop over the source pixels containing brightness in excess of a threshold of ∼0.3 mJy arcsec−2, vary their brightness by a small quantity ± and calculate the new value of . The iteration uses the P18 source brightness as a first approximation. We retain as new value of the brightness the value giving the smaller value of . We repeat the procedure until all pixels contain a brightness S giving a smaller value of than for S±. Convergence is achieved after 60 to 140 iterations depending on the velocity interval. Fig. 4. Upper row: maps of the intensity, integrated over the whole velocity range, of the source emission as obtained by P18 (left), by direct de-lensing (centre) and by χ2 minimization (right). The location of the quasar is indicated by a cross. Units are Jy km s−1 arcsec−2. Ellipses show the projection of the model disc. Middle row: projections of the source brightness (Jy km s−1 arcsec−1) on the major axis of its elliptical projection on the sky plane (x0 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 mean Doppler velocity, together with those obtained by P18 and by direct de-lensing. The three intensity maps of the source emission are consistent with the elliptical projection on the sky plane of a thin circular disc centred on the quasar (x=−0.49 arcsec, y=−0.005 arcsec). We find that the major axis
9. MORPHO-KINEMATICS OF THE MOLECULAR GAS IN A QUASAR HOST GALAXY 157 of the ellipse of projected emission is oriented some 14◦ north of the x axis, meaning 30◦ north of west (L17 quote a value of 31◦). We evaluate the lengths of the major and minor axes to be ∼2.7 and ∼1.6 arcsec (∼19 and ∼11 kpc respectively) corresponding to an inclination of the disc with respect to the plane of the sky of cos−1(1.6/2.7)=54◦ as obtained by P18. We note however that the long axis of the P18 caustic is only ∼12◦ south of east instead of 16◦ in our simple lens model. Also shown 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 maps of the mean Doppler velocity display a strong gradient along the major axis, as expected from a rotating thin disc. L17 have claimed evidence for the presence of a companion of the quasar host galaxy in the red-most velocity interval covering from 249 to 333 km s−1. 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-most velocity interval, centred at x ∼0.8 arcsec and y ∼0.2 arcsec (Figure 5). It is most probably what L17 are referring to, but it does not match precisely image F of Brewer & Lewis 2008 [26], being north rather than south of the western cusp of the caustic. Lensing this lump produces two images, shown as A and B in Figure 5. Contrary to image B, image A stands out of the general image mor- phology; it has a magnification of ∼0.5 while image B has a magnification of ∼4.3. While much weaker than image B, image A contributes as much as image B to the de-lensed source brightness and is therefore responsible for the appearance as an isolated lump in the source plane: the case for a lump of separate emission is much weaker when made in the image plane than in the source plane. The fact that the Doppler velocity of the lump matches well that implied in this region by the galaxy rotation curve argues against the L17 interpretation as a companion galaxy. It is more natural to interpret it as a lump of enhanced emission on the disc. Fig. 5. Maps of the brightness integrated over the red-most velocity interval (249 to 333 km s−1) are shown in the three leftmost panels for the source and in the two rightmost panels 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, the observed images and the best-fit results. Colour scales are in Jy km s−1 arcsec−2. IV. A SIMPLE ROTATING DISC MODEL IV.1. Geometry The preceding sections have shown that the general morpho-kinematics of the de-lensed images is a robust result of the analysis of P18. Using a simpler lens does not strongly affect
10. 158 T. T. THAI et al. 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 a reference with which one can make quantitative comparisons, it is instructive to construct the image produced by a rotating disc of uniform brightness having morpho-kinematics matching the distributions 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 x0, intersection of the disc with the plane of the sky containing the quasar, we see from Figure 6 that x0 = Rcosθ y0 = Rsinθ cos54◦ z0 = Rsinθ sin54◦ ◦ ◦ Vx = −V(R)sinθ Vy = V(R)cosθ cos54 Vz = V(R)cosθ sin54 . (4) In our system of coordinates, with x pointing 16◦ north of west, the x0 axis points to a direction of Fig. 6. Left: geometry in a system of coordinates (x0,y0,z0) centred on the quasar with x0 along the trace of the disc on the plane of the sky and z0 perpendicular to the plane of the sky. Centre and right: the image of a disc of uniform brightness lensed by the simple lens potential (centre) is compared with observation (right). Units are arbitrary for the model and Jy km s−1 arcsec−2 with a cut at 0.67 Jy km s−1 arcsec−2 for the observed data. 30◦ −16◦ =14◦ ; centring the disc on the quasar of coordinates (x,y)=(−0.49, −0.005) arcsec, we obtain: x + 0.49 = x0 cos14◦ − y0 sin14◦ y + 0.005 = x0 sin14◦ + y0 cos14◦. Using the above equations and the results obtained in the preceding section, we image a disc of uniform brightness centred on the quasar, inclined by 54◦ with respect to the plane of the sky and having a radius Rdisc=1.35 arcsec. The inclination of the disc simply divides the disc plane brightness 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 a uniform disc is confined within an approximately elliptical annular band centred on the lens/quasar region. We define two ellipses, E+ and E−, delimiting the outer and respectively inner edges of the uniform disc image as shown in Figure 6. The right panel of Figure 6 displays the image brightness of the observed data. It is well contained within ellipses E+ and E− and populates the middle of the band delimited by them. In order to use coordinates adapted to this morphology, we define a
11. MORPHO-KINEMATICS OF THE MOLECULAR GAS IN A QUASAR HOST GALAXY 159 parameter λ such that λ=−0.5 on E− and +0.5 on E+, namely λ=±0.5 on E±. Precisely, a point on the sky plane having Cartesian coordinates (x = r cosω, y = r sinω) and polar coordinates (r, ω) is imaged as (λ, ω) with λ = [r − (r+ + r−)/2]/(r+ − r−); here, r+ and r− are the points of position angle ω on ellipses E+ and E− respectively. In the remaining of the section we work in this new system of coordinates, (λ, ω, Vz). IV.2. The data cube We construct a data cube in (λ, ω, Vz) coordinates to be compared with model predictions. ◦ −1 We set its limits as −0.5<λ<+0.5, 0<ω<360 and −340<Vz<333 km s ; we segment it in 20 ◦ −1 λ bins, each 0.05 wide, in 18 ω bins, each 20 wide and in 16 Vz bins, each ∼42 km s wide. We construct the new data cube starting from the reduced data presented in Sec. III, distributed in 125×125 pixels in the sky plane, each 50×50 mas2 in area, and 80 Doppler velocity bins, each 8.417 km s−1 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 treatment of the noise. Instead, we compare in Figure 7 the distributions obtained by applying no brightness cut on the original data cube elements to those obtained by applying a cut (∼1.5-σ) such that the flux integrated over the whole data cube is the same [22]. The upper row of Figure 7 compares the 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 ∼1.5-σ) producing the same integrated flux as when no cut is applied gives the results shown in the second row of Figure 7. Finally, the third row compares the present data with the P18 data, both with no brightness cut being applied. In summary, the differences in noise level and angular resolution between the analyses of the present work and of P18, apart from a global rescaling factor (∼0.8) that is irrelevant to the following arguments, are unimportant. Differences between the four data cubes illustrated in Figure 7 help with an evaluation of the uncertainties attached to these measurements, on average ∼13 mJy per histogram bin. IV.3. The disc model R/R∗ We parameterize the rotation curve as V(R) = V (e −1) and the disc brightness is taken 0 (eR/R∗+1) uniform over a disc of mean radius Rdisc smeared radially by a Gaussian of dispersion σdisc.This choice has been made after having explored results obtained using different possible forms. While beam convolution is properly accounted for by the model, no noise contribution is included. We decided to do so after having produced fits accounting for a Gaussian noise contribution mimicking that present in the observed data cube. The results were essentially unaffected. Optimization of the values of the four parameters is made by minimizing a χ2 function describing the quality of the match between model and observations. Rather than evaluating the value of χ2 as a sum over the 20×18×16=5760 data cube elements, we find that restricting the sum to the 20+18+16=54 bins of the histograms displayed in Figure 7 is better adapted to the precision allowed by the quality of the data and the crudity of the model. The χ2 is evaluated using as uncertainty a common arbitrary value of 10 mJy per histogram bin and is divided by the number of degrees of freedom. −1 ∗ The best fit is obtained for V0=405 km s , R =0.22 arcsec (1.6 kpc), Rdisc=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 origin than implied by P18 and L17. The results are illustrated in Figure 8.
12. 160 T. T. THAI et al. Fig. 7. Projections of the new data cube (see text) on λ (left), ω (centre) and Vz (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×50×8.417 mas2 km s−1). 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 the present data (black) with the P18 data (red). The best fit value of χ2 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 for the 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 disc brightness are essentially uncorrelated; but the results obtained for each independently display a ∗ very strong correlation between V0 and R on the one hand and between Rdisc and σdisc on the 2 ∗ other. This is illustrated in Figure 9 which maps the value of χ in the R vs V0 and σdisc vs Rdisc planes respectively. IV.4. Brightness inhomogeneity Maps of the observed and modelled brightness are compared in Figure 10 in the λ vs ω, λ vs Vz and Vz vs ω planes. In general, the agreement between observed and modelled maps is remarkable given the crudeness of the model. The Vz vs ω maps display an oscillation that reveals 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
13. MORPHO-KINEMATICS OF THE MOLECULAR GAS IN A QUASAR HOST GALAXY 161 Fig. 8. Comparison between observations (black) and best fit model (red and blue) are shown as projections of the data cube on λ (left), on ω (centre) and on Vz (right). The red (blue) histograms are with (without) allowance for the presence of an emission hot spot (see Subsec. IV.4). 2 ∗ Fig. 9. Left: dependence of χ on the model parameters: V0 (abscissa) and R (ordinate) in the left-most panels and Rdisc (abscissa) and σdisc (ordinate) in the right-most panels. In each pair of panels, the results obtained without and with allowance for an excess of emission in the vicinity of the quasar are ordered from left to right. In each panel the other model parameters are set at their best-fit values. The values of FE listed in the inserts refer to the excess of emission discussed in Subsec. 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. They are dominantly in the two red-most Vz 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 negative values of λ and a depletion for small positive values. We define two rectangles in the λ vs ω plane that delimit these regions: an excess E (200◦<ω<235◦,−0.18<λ<−0.05) and a depletion D (207◦<ω<237◦,0.15<λ<0.27). The left panels of Figure 11 show the associated regions of the source plane. The source of D is outside the caustic and its magnification spans values between ∼3 and ∼5. The companion image of D is small and centred at (x, y)∼(0.75, 0.25) arcsec; its magnification spans values between 0.5 and 1. In contrast, the source of E is located inside the caustic and
14. 162 T. T. THAI et al. Fig. 10. Maps of the brightness projected on the (λ vs ω), (λ vs Vz) and (Vz vs ω) planes (respectively from left to right) and integrated over the third coordinate, Vz, ω ◦ −1 and λ respectively. Bins are 0.05 wide in λ, 20 wide in ω and 42 km s wide in Vz. Upper panels display observed images with no brightness cut. Central panels display the images obtained from lensing the best fit disc model. Contours of the disc model (central row) are superimposed on observations (upper row). Lower panels map the difference between observations and model. The rectangles in the left panel show the regions of strong 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 ∼4 and extending to infinity with a mean value of ∼20. It produces three additional images, which span magnifications ranging between 1.5 and 5.5; one of these, centred at (x, y)∼(−1, 2) arcsec is in a region where the lower-left panel of Figure 10 shows that the observed brightness exceeds slightly that of the model. The fact that de-lensing E produces a source limited by the caustic curve is not an artefact of the lensing mechanism. In fact, the source that produces E overlaps the caustic but only the part of it that is inside contributes to image E. The part that is outside produces instead 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 reason to be; it is indeed simpler and more natural to blame the inhomogeneity on a single hot spot, in which case the flux predicted by the model needs to be scaled down, or as a single depletion, in which case the flux predicted by the model needs to be scaled up. We tried both interpretations by
15. MORPHO-KINEMATICS OF THE MOLECULAR GAS IN A QUASAR HOST GALAXY 163 including in the model an enhancement (depletion) of emission by a factor FE (FD) in the source region of E (D) shown in Figure 11. In both cases the best fit values of the model parameters, ∗ V0, R , Rdisc and σdisc, were essentially unaffected but allowing for a depletion did not produce lower values of χ2; in contrast, as illustrated in the right panels of Figure 11, allowing for an excess resulted in a significant improvement of the quality of the fit. The best fit values of the −1 ∗ model parameters are now V0=435 km s , R =0.26 arcsec (1.8 kpc), Rdisc=1.10 arcsec (7.7 kpc), σdisc=0.32 arcsec (2.25 kpc) and FE =2.5 (see Figure 9). The above results call for a word of caution. Qualitatively, the presence of an excess of emission in the vicinity of the quasar is likely to be a robust result of the analysis; it had been noted by 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-kinematical properties. The lens model used in the present work is very crude and so is also that used by P18: none of these allows for local inhomogeneity of the lensing potential, they describe it in very general terms; the vicinity of the caustic is particularly difficult to model reliably, small variations of the lensing potential having strongly amplified effects in the construction of the images. As we shall see in the next section, the overall remarkable agreement between model and observations displayed 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 (x0 axis in Figure 6). We choose for it a width of ±1 kpc and a length of 2.7 arcsec (∼19 kpc), divided in nine segments, each having a length of 0.3 arcsec (∼2 kpc). In each of these segments we compare the observed and modelled velocity spectra. The modelled spectra are obtained by de-lensing the images pro- duced by lensing the model disc source and convolved with the beam. Results are illustrated in the left panel of Figure 12. Qualitatively, the general trend is well reproduced by the model but significant differences are observed in the central segments: the data display larger Doppler ve- locities on the red side and lower Doppler velocities on the blue side than implied by the model. Moreover, in the central segment, the line width predicted by the model is much smaller than that observed in the data [22]. A natural interpretation of such an effect is disc warping causing an effective dependence on θ of the sine of the inclination angle in Relation (4). However, including warping in the model by writing Vz = V(R)cosθ sinϕ with ϕ depending simply on θ and R, gives only a modest improvement of the match between model and observations. This suggests that a more 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 the nine segments, making the radial distribution of the mean Doppler velocity measured in the central segments less steep and causing an increase of the velocity dispersion. This important result is a warning: measuring the velocity dispersion requires a careful evaluation of the contribution of rotation within the finite angular resolution of relevance. Making the angular acceptance smaller will decrease the direct contribution of rotation, independent from beam convolution, but will not decrease the contribution resulting from the smearing caused by the beam. This contribution is important: using the central segment as an example, the line width predicted by the model increases from 60 km s−1 to nearly 100 km s−1 when accounting for beam convolution, namely a beam contribution of over 70 km s−1 (both contributions add up in quadrature).
16. 164 T. T. THAI et al. Fig. 11. Left panels display the brightness distribution in the source plane obtained by de-lensing depletion D (left) and excess E (centre-left) respectively. Right panels display 2 the dependence of χ on the factors FD and FE measuring the amplitudes of the depletion (centre-right) and excess (right) respectively when the other model parameters are fixed at their best-fit values. 0 Fig. 12. Dependence of on x (major axis of the projection of the disc on the plane of the sky) for the data (black) and the model (red) is shown in the left panel and the rotation 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 different rotation 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 enhanced emission studied in the preceding section is the seat of increased gas turbulence. However, the differences observed in the central segments of the band used to evaluate the rotation curve do not consist of a simple broadening of the line, as would be expected from turbulence, but reveal rather an 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 Doppler velocities tend to be too small in the 90◦<θ<180◦ quadrant and too large in the 270◦<θ<360◦ quadrant of the disc plane; velocity dispersions reach some 120 km s−1 for the model and some
17. MORPHO-KINEMATICS OF THE MOLECULAR GAS IN A QUASAR HOST GALAXY 165 Fig. 13. Line profiles in the region |λ| i j of the Doppler velocity and show the distribution of Vz− i j in Figure 13. We obtain this way an evaluation of the line profile. The observed line widths (σ) are 85 and 56 km s−1 for cuts at 11.7 (1.5-σ) and 16 µJy (2-σ) per 50×50×8.417 mas2 km s−1 respectively. As the model includes no intrinsic line broadening, the line broadening shown by the model is entirely due to the contribution of rotation. Gaussian fits give σ values of 18 and 58 km s−1 without and with beam convolution respectively, implying a contribution of ∼50 km s−1 from the smearing effect of the beam size. The line width measured with the 16 µJy cut does not require any contribution of turbulence while that measured with an 11.7 µJy cut allows for a contribution reaching ∼60 km s−1: without a precise understanding of the contribution of noise, it is therefore difficult to obtain a reliable evaluation of the line broadening caused by turbulence. In summary, a simple rotating disc model gives a global picture that describes well the general trend of the observed kinematics. However, such a model fails to reproduce quantitatively the details of the Doppler velocity distribution: it reveals a more complex dynamics. Attempts to describe the mismatch as the result of a simple disc warping have failed. An important result of our analysis is that the smearing in the image plane caused by the beam size contributes between 50 and 70 km s−1 to the dispersion of the Doppler velocity. It combines in quadrature with the contribution 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 to the line width. But they prevent claiming that such contribution is important: within experimental
18. 166 T. T. THAI et al. uncertainties the line width can be accounted for by the dispersion of rotation velocities within the beam size. V. SUMMARY AND CONCLUSIONS The present study of the emission of the CO(2-1) molecular line by the host galaxy of quasar RX J1131 uses ALMA observations of unprecedented quality; they have been previously analysed in much detail by their proponents, [18], P18. We have shown that the HST images of the quasar and of the lens galaxy are very well reproduced by a simple lensing potential, sphere+external shear. We obtained parameters that are in agreement with the results of previous authors within errors. We discussed the uncertainties attached to the model parameters and their correlations and explained them by remarking that these results probe only the environment of the eastern cusp of the caustic. At this stage, the validity of the model at larger distances, in the region covered by the emission of the host galaxy, cannot be taken as granted. But we demonstrated its validity from our analysis of the ALMA observations of the CO(2-1) line emission. We reduced the raw data and produced cleaned images, at variance with P18 who perform their analysis in the uv plane. We obtained morpho-kinematics properties of the gas that are in good agreement with the results of P18. We used two different methods for de-lensing these images, which produce source brightness distributions in agreement with each other as well as with the P18 results. We have shown that the general morphology of the image is dictated by the properties of the lens and confined within a band bracketing the critical curve. This leaves in practice no freedom to conceive lensing potentials that would be significantly different from those that we and other authors have been using. It guarantees the robustness of the results obtained by P18 within observational uncertainties. We used polar coordinates better adapted to this geometry than Cartesians to reveal very clearly the evidence for rotation. We looked for possible deviations from the predictions of a simple rotating disc model, in particular for the possible presence of a companion as was predicted by L17 using Plateau de Bure observations of much less good angular resolution than the ALMA observations; 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 for enhanced 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 important uncertainties on its de-lensed properties. We discussed the rotation curve in some detail and argued that it probably rises faster than assumed by P18 in the vicinity of the quasar; we noted that the angular resolution is such that the dispersion of the rotation velocity within the beam size in the image plane causes an important effective broadening of the line width, which we evaluated in the range between 50 and 70 km s−1 . This is sufficient to account for the observed line width when applying a ∼2-σ cut on the data. However, uncertainties attached to the effect of noise prevent from giving a reliable evaluation of the contribution of turbulence. This is at variance with the conclusion of P18 who attribute the high Doppler 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 of the galaxy, which they compared with those of other similar galaxies. The nature of our study prevents us from adding much to their conclusions. We simply remark that their estimate of 11 the total dynamical mass enclosed within 5 kpc, (1.46±0.31)×10 M, may be affected by the steepness of the rotation curve: at fixed radius, the dynamical mass is proportional to the square of
19. MORPHO-KINEMATICS OF THE MOLECULAR GAS IN A QUASAR HOST GALAXY 167 the rotation velocity. The analysis presented in the preceding sections implies therefore a possibly higher value of the gas mass, obtained by P18 and accordingly affect the value of the CO−H2 conversion factor. ACKNOWLEDGMENT We express our deep gratitude to Professors Frederic Courbin and Matus Rybak who kindly provided us with copies of data files summarizing the results of the P18 analysis. This paper uses archival ALMA data from project 2013.1.01207.S (PI: Paraficz Danuta). ALMA is a partnership of 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 is operated 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] R. Adam, N. Aghanim, M. Ashdown, J. Aumont, C. Baccigalupi, M. Ballardini, A. Banday, R. Barreiro, N. Bar- tolo, S. Basak and R. Battye, A&A 596 (2016) A108. [2] P. A. Ade, N. Aghanim, M. Arnaud, M. Ashdown, J. Aumont, C. Baccigalupi, A. Banday, R. Barreiro, J. Bartlett, N. Bartolo et al., A&A 594 (2016) A13. [3] R. C. Reis, M. T. Reynolds, J. M. Miller and D. J. Walton, Nature 507 (2014) 207. [4] M. Tewes, F. Courbin, G. Meylan, C. Kochanek, E. Eulaers, N. Cantale, A. Mosquera, P. Magain, H. Van Winckel, D. Sluse and G. Cataldi, A&A 556 (2013) A22. [5] V. Bonvin, F. Courbin, S. H. Suyu, P. J. Marshall, C. E. Rusu, D. Sluse, M. Tewes, K. C. Wong, T. Collett, C. D. Fassnacht and T. Treu, Monthly Notices of the Royal Astronomical Society 465 (2017) 4914. [6] G. Chartas, C. S. Kochanek, X. Dai, D. Moore, A. M. Mosquera and J. A. Blackburne, ApJ 757 (2012) 137. [7] H. Sugai, A. Kawai, A. Shimono, T. Hattori, G. Kosugi, N. Kashikawa, K. T. Inoue and M. Chiba, ApJ 660 (2007) 1016. [8] P. Anh, F. Boone, D. Hoai, P. Nhung, A. Weiß, J.-P. Kneib, A. Beelen and P. Salome,´ A&A 552 (2013) L12. [9] A. P.T., PhD thesis,Vietnam Institute of Physics and Touluse University (2014) . [10] H. R. Stacey, J. P. McKean, N. C. Robertson, R. J. Ivison, K. G. Isaak, D. R. G. Schleicher, P. P. van der Werf, W. A. Baan, A. Berciano Alba, M. A. Garrett and A. F. Loenen, Mon. Notices Royal Astron. Soc. 476 (2018) 5075. [11] T. D. Leung, D. A. Riechers and R. Pavesi, ApJ 836 (2017) 180. [12] CASTLES, Summary of RXJ1131-1231 observations, RXJ1131.html, 2015, Mon Oct 26 19:54:39 EDT 2015. [13] J.-F. Claeskens, D. Sluse, P. Riaud and J. Surdej, A&A 451 (2006) 865. [14] S. Birrer, A. Amara and A. Refregier, J. Cosmol. Astropart. Phys. 2016 (2016) 020. [15] G. C.-F. Chen, S. H. Suyu, K. C. Wong, C. D. Fassnacht, T. Chiueh, A. Halkola, I. S. Hu, M. W. Auger, L. V. Koopmans, D. J. Lagattuta et al., Mon. Notices Royal Astron. Soc. 462 (2016) 3457. [16] R. Blandford, C. Kochanek, I. Kovner and R. Narayan, Science 245 (1989) 824. [17] P. Saha and L. L. R. Williams, AJ 125 (2003) 2769. [18] D. Paraficz, M. Rybak and J. P. McKean, A&A A34 (2018) 613. [19] M. Rybak, J. P. McKean, S. Vegetti, P. Andreani and S. D. M. White, Mon. Notices Royal Astron. Soc. 451 (2015) L40. [20] M. Rybak, S. Vegetti, J. P. McKean, P. Andreani and S. D. M. White, Mon. Notices Royal Astron. Soc. 453 (2015) L26. [21] S. Vegetti and L. V. E. Koopmans, Mon. Notices Royal Astron. Soc. 945 (2009) 392.
20. 168 T. T. THAI et al. [22] T. T. Thai, Master thesis, Graduate University of Science and Technology (2020) . [23] P. T. Nhung, P. T. Anh, F. Boone, P. Darriulat, P. N. Diep, P. N. Dong, N. Van Hiep and N. T. Thao, Res. Astron. Astrophys. 13 (2013) 803. [24] S. H. Suyu, M. W. Auger and S. Hilbert, ApJ 70 (2013) 766. [25] D. Sluse, J.-F. Claeskens and B. Altieri, A&A 539 (2006) 449. [26] B. B. J. and G. F. Lewis, Mon. Notices Royal Astron. Soc. 39 (2008) 390.