Modelling high resolution ALMA observations of strongly lensed highly star forming galaxies detected by Herschel

TL;DR

This paper models high-resolution ALMA images of six strongly lensed galaxies detected by Herschel, revealing their intrinsic properties and lens mass profiles, and comparing modeling approaches.

Contribution

It introduces a detailed modeling approach for ALMA data of lensed galaxies, recovering lens and source properties and comparing image plane and visibility modeling methods.

Findings

All galaxies have high star formation rate to dust mass ratios.

Most galaxies show disturbed morphologies.

Lens mass profiles are close to isothermal in most cases.

Abstract

We have modelled high resolution ALMA imaging of six strong gravitationally lensed galaxies detected by the Herschel Space Observatory. Our modelling recovers mass properties of the lensing galaxies and, by determining magnification factors, intrinsic properties of the lensed sub-millimetre sources. We find that the lensed galaxies all have high ratios of star formation rate to dust mass, consistent with or higher than the mean ratio for high redshift sub-millimetre galaxies and low redshift ultra-luminous infra-red galaxies. Source reconstruction reveals that most galaxies exhibit disturbed morphologies. Both the cleaned image plane data and the directly observed interferometric visibilities have been modelled, enabling comparison of both approaches. In the majority of cases, the recovered lens models are consistent between methods, all six having mass density profiles that are close to isothermal. However, one system with poor signal to noise shows mildly significant differences.

Tables

Table 1. Table 1: The six lenses systems modelled in this work with their lens galaxy redshifts, z l subscript 𝑧 𝑙 z_{l} , and source redshifts, z s subscript 𝑧 𝑠 z_{s} . a Bussmann et al. ( 2012 ) . b Cox et al. ( 2011 ) . c Messias et al. ( 2014 ) . d Negrello et al. ( 2017 ) . e Nayyeri et al. ( 2016 ) . f Amvrosiadis et al. ( 2017 ) . g Marchetti et al. (in prep.).

ID	$z_l$	$z_s$
H-ATLAS J142413.9+022303	0.595a	4.243b
H-ATLAS J142935.3-002836	0.218c	1.026d
HELMS J004714.2+032454	0.478e	1.190e
HELMS J001626.0+042613	0.215f,g	2.509e
HELMS J004723.6+015751	0.365f,g	1.441e
HELMS J001615.7+032435	0.663e	2.765e

Table 2. Table 2: Observed (i.e., lensed) source flux densities in mJy. Subscripts indicate the passband central wavelength in μ 𝜇 \mu m. Fluxes f 100 subscript 𝑓 100 f_{100} to f 500 subscript 𝑓 500 f_{500} inclusive are taken from the H-ATLAS first data release (Valiante et al., 2016 ) for the two H-ATLAS sources. For the four HELMS sources, f 100 subscript 𝑓 100 f_{100} and f 160 subscript 𝑓 160 f_{160} are PACS flux densities extracted from maps acquired from the Herschel Science Archive and flux densities f 250 subscript 𝑓 250 f_{250} to f 500 subscript 𝑓 500 f_{500} are taken from Nayyeri et al. ( 2016 ) . Flux densities f 850 subscript 𝑓 850 f_{850} , f 880 SMA superscript subscript 𝑓 880 SMA f_{880}^{\rm SMA} , f 880 ALMA superscript subscript 𝑓 880 ALMA f_{880}^{\rm ALMA} and f 1280 subscript 𝑓 1280 f_{1280} are taken from Bakx et al. (2017, in prep.), Bussmann et al. ( 2013 ) , this work and Messias et al. ( 2014 ) respectively. Finally, f 60 subscript 𝑓 60 f_{60} is the 60 μ 𝜇 \,\mu m flux taken from the IRAS faint source catalogue (Moshir, Kopman & Conrow, 1992 ) .

ID	$f_{60}$	$f_{100}$	$f_{160}$	$f_{250}$	$f_{350}$	$f_{500}$	$f_{850}$	$f_{880}^{\mathrm{SMA}}$	$f_{880}^{\mathrm{ALMA}}$	$f_{1280}$
H-ATLAS J142413.9+022303	-	-	-	$112\pm 7$	$182\pm 8$	$193\pm 8$	$121\pm 8$	$90\pm 5$	$116\pm 8$	-
H-ATLAS J142935.3-002836	$190\pm 38$	$911\pm 29$	$1254\pm 34$	$802\pm 7$	$438\pm 7$	$200\pm 7$	-	-	$38\pm 3$	$5.86\pm 0.99$
HELMS J004714.2+032454	-	$82\pm 11$	$164\pm 22$	$312\pm 6$	$244\pm 7$	$168\pm 8$	-	-	$49\pm 5$	-
HELMS J001626.0+042613	-	$13\pm 10$	$53\pm 20$	$117\pm 7$	$151\pm 6$	$127\pm 7$	-	-	$39\pm 4$	-
HELMS J004723.6+015751	-	$104\pm 15$	$285\pm 32$	$398\pm 6$	$320\pm 6$	$164\pm 8$	-	-	$42\pm 5$	-
HELMS J001615.7+032435	-	$23\pm 11$	$92\pm 24$	$195\pm 6$	$221\pm 6$	$149\pm 7$	-	-	$33\pm 4$	-

Table 3. Table 3: Lens model parameters. The top half of the table gives the parameters obtained from the image plane analysis and the bottom half gives those from the visibility plane analysis. Only HELMS J004723.6+015751 and HELMS J001615.7+032435 showed significant improvement in the fit when external shear was included in the lens model, hence the remaining four were modelled without it. Parameters are: lens normalisation, κ 0 subscript 𝜅 0 \kappa_{0} in units of 10 10 ​ M ⊙ ​ kpc − 2 superscript 10 10 subscript 𝑀 direct-product superscript kpc 2 10^{10}M_{\odot}\,{\rm kpc}^{-2} ; co-ordinates of the lens model centroid with respect to the phase-tracking centre of observations (west and north correspond to positive x c subscript 𝑥 𝑐 x_{c} and y c subscript 𝑦 𝑐 y_{c} respectively); lens semi-major axis orientation, θ 𝜃 \theta , measured counter-clockwise from north; lens semi-major to semi-minor axis ratio, ϵ italic-ϵ \epsilon ; logarithmic slope of the power-law density profile, α 𝛼 \alpha ; external shear strength, γ 𝛾 \gamma ; shear direction angle, θ γ subscript 𝜃 𝛾 \theta_{\gamma} , measured counter-clockwise from north; Einstein radius, θ E subscript 𝜃 𝐸 \theta_{E} .

ID	${\kappa }_0$	$(x_c,y_c)$ (arcsec)	$\theta \left(\mathrm{deg}\right)$	$ϵ$	$\alpha$	$\gamma$	${\theta }_{\gamma }\left(\mathrm{deg}\right)$	${\theta }_E\left(\mathrm{arcsec}\right)$
Image plane
H-ATLAS J142413.9	$0.59\pm 0.01$	$(0.18\pm 0.01,0.68\pm 0.01)$	$84\pm 2$	$1.07\pm 0.02$	$2.03\pm 0.04$			$0.97\pm 0.04$
H-ATLAS J142935.3	$0.44\pm 0.01$	$(1.60\pm 0.01,0.62\pm 0.01)$	$124\pm 1$	$1.33\pm 0.02$	$1.82\pm 0.05$			$0.71\pm 0.03$
HELMS J004714.2	$0.50\pm 0.01$	$(1.56\pm 0.02,2.34\pm 0.03)$	$94\pm 2$	$1.25\pm 0.02$	$1.96\pm 0.04$			$0.59\pm 0.03$
HELMS J001626.0	$0.56\pm 0.01$	$(2.88\pm 0.02,1.67\pm 0.02)$	$36\pm 1$	$1.37\pm 0.03$	$2.14\pm 0.06$			$0.98\pm 0.07$
HELMS J004723.6	$1.18\pm 0.02$	$(2.52\pm 0.02,-0.60\pm 0.02)$	$178\pm 2$	$1.18\pm 0.01$	$1.87\pm 0.04$	$0.09\pm 0.01$	$167\pm 2$	$2.16\pm 0.10$
HELMS J001615.7	$2.21\pm 0.04$	$(0.12\pm 0.05,-0.96\pm 0.07)$	$18\pm 2$	$1.41\pm 0.02$	$2.00\pm 0.07$	$0.13\pm 0.01$	$55\pm 2$	$2.79\pm 0.20$
Visibility plane
H-ATLAS J142413.9	$0.59\pm 0.01$	$(0.18\pm 0.01,0.68\pm 0.01)$	$85\pm 2$	$1.07\pm 0.01$	$2.06\pm 0.04$			$0.97\pm 0.04$
H-ATLAS J142935.3	$0.43\pm 0.01$	$(1.60\pm 0.01,0.61\pm 0.01)$	$125\pm 1$	$1.35\pm 0.02$	$1.79\pm 0.05$			$0.70\pm 0.03$
HELMS J004714.2	$0.50\pm 0.01$	$(1.55\pm 0.02,2.34\pm 0.03)$	$93\pm 2$	$1.24\pm 0.02$	$1.91\pm 0.05$			0.$58\pm 0.03$
HELMS J001626.0	$0.58\pm 0.01$	$(2.89\pm 0.02,1.66\pm 0.02)$	$36\pm 1$	$1.38\pm 0.03$	$2.18\pm 0.06$			$0.98\pm 0.07$
HELMS J004723.6	$1.18\pm 0.03$	$(2.51\pm 0.02,-0.60\pm 0.02)$	$178\pm 2$	$1.20\pm 0.02$	$1.89\pm 0.06$	$0.08\pm 0.01$	$161\pm 2$	$2.08\pm 0.12$
HELMS J001615.7	$2.00\pm 0.07$	$(0.11\pm 0.05,-0.94\pm 0.06)$	$18\pm 2$	$1.42\pm 0.02$	$1.90\pm 0.05$	$0.10\pm 0.01$	$53\pm 2$	$2.96\pm 0.16$

Table 4. Table 4: Intrinsic source properties. Columns are the total source magnification computed using the image plane method and uv-plane method, μ tot img subscript superscript 𝜇 img tot \mu^{\rm img}_{\rm tot} and μ tot uv subscript superscript 𝜇 uv tot \mu^{\rm uv}_{\rm tot} respectively, dust mass assuming a single temperature optically thick SED, M d thick superscript subscript 𝑀 𝑑 thick M_{d}^{\rm thick} , dust mass assuming a dual temperature optically thin SED, M d thin superscript subscript 𝑀 𝑑 thin M_{d}^{\rm thin} , temperature of the optically thick SED, T thick superscript 𝑇 thick T^{\rm thick} , temperatures of the optically thin SED, T thin / K superscript 𝑇 thin 𝐾 T^{\rm thin}/K , the opacity at 100 ​ μ 100 𝜇 100\,\mu m for the optically thick SED, τ 100 subscript 𝜏 100 \tau_{100} , de-magnified luminosity (computed as the integral of the best fit SED from 3 to 1100  μ 𝜇 \mu m using the optically thin SED), L FIR subscript 𝐿 FIR L_{\rm FIR} , H 2 gas mass calculated using the scaling relation of Hughes et al. ( 2017 ) , M gas subscript 𝑀 gas M_{\rm gas} , and star formation rate (SFR) scaled from L FIR subscript 𝐿 FIR L_{\rm FIR} using the prescription given by Kennicutt & Evans ( 2012 ) with a Kroupa IMF. Dust masses are expressed as log 10 ⁡ ( M d / M ⊙ ) subscript 10 subscript 𝑀 𝑑 subscript 𝑀 direct-product \log_{10}(M_{d}/M_{\odot}) , gas masses as log 10 ⁡ ( M gas / M ⊙ ) subscript 10 subscript 𝑀 gas subscript 𝑀 direct-product \log_{10}(M_{\rm gas}/M_{\odot}) and the luminosity values are log ⁡ ( L FIR / L ⊙ ) subscript 𝐿 FIR subscript 𝐿 direct-product \log(L_{\rm FIR}/L_{\odot}) .

ID	${\mu }_{\mathrm{tot}}^{\mathrm{img}}$	${\mu }_{\mathrm{tot}}^{\mathrm{uv}}$	$M_d^{\mathrm{thick}}$	$M_d^{\mathrm{thin}}$	$T^{\mathrm{thick}}/K$	$T^{\mathrm{thin}}/K$	${\tau }_{100}$	$L_{\mathrm{FIR}}$	$M_{\mathrm{gas}}$	SFR ($M_{\odot }$/yr)
H-ATLAS J142413.9	$6.6\pm 0.5$	$6.4\pm 0.5$	8.7	9.7	59	41 / 21	5.8	$13.2\pm 0.1$	$11.8\pm 0.1$	$2200\pm 500$
H-ATLAS J142935.3	$23.6\pm 1.3$	$22.3\pm 1.3$	7.9	8.2	70	45 / 26	4.4	$12.3\pm 0.1$	$10.7\pm 0.1$	$330\pm 80$
HELMS J004714.2	$8.3\pm 0.6$	$8.7\pm 0.6$	8.7	9.2	43	51 / 22	9.2	$12.2\pm 0.1$	$11.3\pm 0.1$	$220\pm 60$
HELMS J001626.0	$4.1\pm 0.3$	$4.3\pm 0.3$	8.8	9.3	48	57 / 27	4.4	$12.8\pm 0.1$	$11.5\pm 0.1$	$980\pm 240$
HELMS J004723.6	$16.5\pm 1.0$	$15.2\pm 1.0$	8.2	8.7	52	48 / 26	5.2	$12.2\pm 0.1$	$10.9\pm 0.1$	$230\pm 60$
HELMS J001615.7	$15.9\pm 1.0$	$17.1\pm 1.0$	7.9	8.5	58	72 / 34	2.4	$12.5\pm 0.1$	$10.7\pm 0.1$	$480\pm 100$

Equations

$${\chi }^2=\sum _{j=1}^J\frac{{\sum }_{i=1}^I{\left|s_i{g}_{ij}-V_j\right|}^2}{{\sigma }_j^2},$$

$$\mathbf{s}={\mathbf{F}}^{-1}\mathbf{v}$$

$${\mathrm{F}}_{ij}$$

$${\mathrm{v}}_i$$

$$\kappa ={\kappa }_0(\tilde{r}/1kpc{)}^{1-\alpha }.$$

Full text

††pagerange: Modelling high resolution ALMA observations of strongly lensed highly star forming galaxies detected by Herschel ††thanks: Herschel is an ESA space observatory with science instruments provided by European-led Principal Investigator consortia and with important participation from NASA.–References††pubyear: 2017

1 Introduction

The most prodigious star formation rates observed in the Universe are located within strongly optically obscured galaxies at high redshift (e.g., Alexander et al., 2005; Greve et al., 2005; Tacconi et al., 2006; Pope et al., 2008). The ultra-violet radiation emitted by their hot young stars is absorbed by copious quantities of enshrouding dust and re-emitted in the mid- and far-infrared (far-IR). Observations indicate that on average they are substantially more energetic per unit mass than local star forming galaxies and have higher star formation efficiencies (e.g., Santini et al., 2014). They are also considerably more abundant than local ultra-luminous infra-red galaxies (ULIRGs) which have comparable bolometric luminosities (e.g., Chapman et al., 2005; Swinbank et al., 2010; Alaghband-Zadeh et al., 2012; Rowlands et al., 2014). Capturing these systems in the midst of a high rate of assembly is of key importance for a complete understanding of galaxy formation. Thanks to recent advances in sub-millimetre (submm) interferometric imaging capability with facilities such as the Atacama Large Millimetre/submillimeter Array (ALMA), study of these high redshift submm-bright galaxies can now be conducted with resolutions $<0.1$ arcsec, providing vastly more detail than was previously possible.

Strong gravitational lensing offers an additional increase in spatial resolution, with magnification factors often in excess of 10. This neatly complements the high lensing bias that occurs at submm wavelengths, which makes selection of strong lens systems relatively easy (Blain, 1996; Negrello et al., 2007). In this way ALMA follow-up of significant numbers of strongly lensed far-IR sources detected in large area surveys such as the Herschel Astrophysical Terahertz Large Area Survey (H-ATLAS; Eales et al., 2010), the Herschel Extragalactic Multi-tiered Extragalactic Survey (HerMES; Oliver et al., 2012) and the Herschel Stripe 82 Survey (HerS Viero et al., 2014) conducted using the Herschel Space Observatory (Pilbratt et al., 2010) and the millimetre wavelength surveys carried out by the South Pole Telescope (Carlstrom et al., 2011; Vieira et al., 2013) and the Planck satellite (Cañameras et al., 2015) are beginning to bring about rapid progress in our understanding of the early stages of galaxy formation. In particular, the improved sensitivity of these facilities allows study of less luminous galaxies than previously possible, pushing down towards the main sequence of star formation occupied by more typical star forming systems.

Not only are these surveys quickly increasing the size of current strong lens samples (e.g., Wardlow et al., 2013; Hezaveh et al., 2013; Bussmann et al., 2013; Calanog et al., 2014; Rowan-Robinson et al., 2014; Bussmann et al., 2015; Nayyeri et al., 2016; Negrello et al., 2017), they are also extending their redshift range owing to the more favourable submm K-correction than that which occurs at shorter wavelengths. Due to the scaling of the lensing cross-section with lens redshift, higher redshift sources are lensed by higher redshift lenses on average and so the extended redshift range also allows study of lens mass profiles in galaxies at an earlier epoch, to widen the time period over which structural evolution in lens galaxies can be studied. Submm lens samples therefore allow the density profile slope to be measured at earlier times when galaxies were evolving more quickly (see, for example, Dye et al., 2014; Negrello et al., 2014).

One particular measurement which has generated significant interest owing to its simplicity and because it provides an observational benchmark for simulations of large scale structure is that of the mass profile of lens galaxies on scales where baryons often dominate the mass budget (i.e., on scales of the Einstein radius; see, for example, Ruff et al., 2011; Bolton et al., 2012; Barnabé et al., 2012; Sonnenfeld et al., 2015). The physics governing the baryons is complex and this gives rise to significant uncertainties in simulations. Observational characterisation of the way in which baryons shape the central mass profile of galaxies therefore brings valuable insight to this problem.

The more accurate lens models afforded by higher resolution submm follow-up also bring about improvements in model-dependent source characteristics such as luminosity, star formation rate and gas and dust mass but also emission line ratios, source morphology and source kinematics which are subjected to differential magnification effects in the reconstructed source plane. A striking example of the degree to which enhancements to our understanding of submm sources can be made by strong lensing can be found in several studies which recently analysed ALMA follow-up imaging of the H-ATLAS discovered lens system SDP81 (see Dye et al., 2015; Swinbank et al., 2015; Rybak et al., 2015a, b; Wong, Suyu & Matsushita, 2015; Tamura et al., 2015; Hezaveh et al., 2016; Inoue et al., 2016). These studies serve to illustrate how high resolution submm imaging brings about a dramatically different interpretation of the lensed source compared to what is inferred from optical data. Whilst significant differences between optical and submm observations, such as large offsets in flux centroids, are not limited to lensed sources, (see, for e.g., Hodge et al., 2015; Chen et al., 2015), differences are expected to be more prevalent at higher redshifts when the rate of galaxy evolution and assembly was higher. At these redshifts, lensing efficiency and therefore lens magnification is high, enabling much enhanced spatial resolution for more detailed morphological study.

Techniques to reconstruct the lensed source from interferometric data naturally divide into those which directly model the visibilities in the uv-plane (e.g., Bussmann et al., 2012, 2013; Rybak et al., 2015a; Hezaveh et al., 2016) and those which model the cleaned data in the image plane (e.g., Dye et al., 2015; Inoue et al., 2016). The advantage of the latter approach is that the reconstruction is often vastly less computationally intensive but this comes at a price of not working with the purest form of the data. This can in principle cause biases in the lens modelling, especially when coverage of the uv-plane is sparse.

In this paper, we have opted to use both uv-plane and image-plane modelling, so that comparison between both methods can be made. We carry out lens modelling of ALMA imaging of six galaxy-galaxy strong lens systems originally detected by the Herschel space observatory within H-ATLAS and the HerMES Large Mode Survey (HELMS; Asboth et al., 2016; Nayyeri et al., 2016) which is an extension to the original HerMES fields.

The layout of this paper is as follows: Section 2 outlines the data. In Section 3 we describe the methodology of the lens modelling. Section 4 presents the results and we summarise the findings of this work in Section 5. Throughout this paper, we assume the following cosmological parameters; ${\rm H}_{0}=67\,{\rm km\,s}^{-1}\,{\rm Mpc}^{-1}$, $\Omega_{m}=0.32$, $\Omega_{\Lambda}=0.68$ (Planck Collaboration, 2013).

2 Data

The ALMA observations modelled in this paper are contained within the ALMA dataset ADS/JAO.ALMA#2013.1.00358.S (PI: Eales). The ALMA spectral setup used for each lens system is identical, comprising Band 7 continuum observations in four spectral windows, each of width 1875 MHz centred on the frequencies 336.5, 338.5, 348.5 and 350.5 GHz. In each spectral window, there are 128 frequency channels giving a resolution of 15.6 MHz. Forty two 12 m antennas were used with an on-source integration time of approximately 125 s. This results in an angular resolution of 0.12 arcsec and an RMS of approximately 230 $\mu$Jy/beam and 130 $\mu$Jy/beam for the H-ATLAS and HELMS sources respectively after combining all four spectral windows. In this paper, we have used the calibrated visibilities as provided in the ALMA science archive. The cleaned data used for the image plane modelling were constructed using Briggs weighting with a robustness parameter of -0.2 and were primary beam corrected. Both calibration and cleaning were carried out using version 4.3.1 of the Common Astronomy Software Applications package (McMullin et al., 2007). The image pixel scale used for the H-ATLAS and HELMS sources was 0.02 and 0.03 arcsec respectively.

When calculating intrinsic source properties, in addition to the photometry obtained from our own ALMA imaging data, we have drawn from a variety of other datasets. We have used submm photometry obtained by the Herschel space observatory using both the Spectral and Photometric Imaging Receiver (SPIRE Griffin et al., 2010) at the wavelengths $250$, $350$ and $500$ $\mu$m and the Photoconductor Array Camera and Spectrometer (PACS; Poglitsch et al., 2010) at wavelengths of 100 and $160$ $\mu$m. For the H-ATLAS sources, SPIRE and PACS photometry was taken from the H-ATLAS first data release (Valiante et al., 2016). For the HELMS sources, SPIRE fluxes were taken from Nayyeri et al. (2016, N16 hereafter) whereas PACS fluxes were extracted from imaging held in the Herschel Science Archive111http://archives.esac.esa.int/hsa/whsa. Where available, we have also used $880\,\mu$m photometry obtained with the Submillimeter Array (SMA) as detailed in Bussmann et al. (2013), $850\,\mu$m Submillimeter Common User Bolometer Array 2 fluxes as given in Bakx et al. (2017, in prep.) and ALMA Band 6 data (1280 $\mu$m) from Messias et al. (2014). Finally, the source H-ATLAS J142935.3-002836 is the Infrared Astronomical Satellite (IRAS) source IRAS 14269-0014 for which we have taken the 60$\,\mu$m flux density as given in the IRAS faint source catalogue (Moshir, Kopman & Conrow, 1992).

Table 1 lists the six systems modelled in this paper along with their lens and source redshifts. Table 2 gives their observed photometry.

3 Methodology

In this paper, we have applied the standard image plane version of the Warren & Dye (2003) semi-linear inversion (SLI) lens modelling method and a modified version which works directly in the interferometric uv-plane on the visibility data. Both use the framework derived by Suyu et al. (2006) for optimising the model Bayesian evidence. The image plane version adopts an implementation similar to that described by Nightingale & Dye (2015) which uses a randomised Voronoi tessellation in the source plane to minimise biases in the lens model parameters. The only differences are that here we have used k-means clustering for the source pixels and Markov Chain Monte Carlo (MCMC) optimisation, whereas Nightingale & Dye used h-means clustering and MultiNest (Feroz, Hobson & Bridges, 2009). The uv-plane version is described in more detail below.

3.1 Adapting the SLI method to visibility data

At the heart of the SLI method lies a pixelised source plane. Using a given lens model, an image of each pixel is formed. In the image plane version of the method, the source surface brightness distribution for a given lens model is determined by finding the linear superposition of these images which best fits the observed lensed image. Adapting this scheme to work with interferometric visibility data requires forming a model visibility dataset for each source pixel image. The linear combination of each model visibility dataset that best fits the observed visibilities then recovers the source surface brightness distribution for a given lens model, in the same manner as the image plane SLI version.

This scheme was used recently by Hezaveh et al. (2016) in application to ALMA data. In their implementation, phase calibration was included in the modelling procedure by introducing the phase offset of each antenna as a free parameter of the fit. In our implementation, the sources are too faint to provide such self-calibration hence we have instead opted to apply the phase calibration provided by external calibrators observed throughout acquisition of our science data.

In the image plane SLI method, the rectangular matrix $f_{ij}$ holds the fluxes of lensed image pixels $j$ for each source plane pixel $i$ assuming the source pixel has unit surface brightness. Analogously, in the uv-plane version, the rectangular matrix $g_{ij}$ is used instead, where each row holds the complex visibilities determined from the lensed image of the unit surface brightness source pixel. Each row of $g_{ij}$ therefore contains the Fourier transform of its corresponding row in $f_{ij}$, evaluated at the same points on the uv-plane as the observed visibilities. This is achieved by incorporating the MIRIAD software package library (Sault, Teuben & Wright, 1995) into our reconstruction code, but using a much streamlined version of the uvmodel procedure. The inputs to uvmodel are the observed visibility dataset and, in turn, the lensed images of the source plane pixels. In this way, a model visibility dataset is created with visibilities equal to $\sum_{i}\,s_{i}g_{ij}$ for each visibility $j$ given source pixel surface brightnesses $s_{i}$. With observed complex visibilities $V_{j}$, the $\chi^{2}$ statistic is therefore computed as

[TABLE]

where the summations act over $I$ total Voronoi source pixels and $J$ visibilities and it is assumed that there is no covariance between visibilities. We used a similar method as Hezaveh et al. (2016) for determining the $1\sigma$ uncertainties, $\sigma_{j}$, on the visibilities. These were computed from the rms of differences in neighbouring visibilities grouped in the uv-plane to remove sky contribution. Whereas Hezaveh et al. computed this for each baseline, our computation was applied over all baselines although our analysis excluded baselines flagged as being bad (and therefore exceptionally noisy) by the ALMA data reduction pipeline. The minimum $\chi^{2}$ solution is given by

[TABLE]

where the elements of the real quantities ${\rm\mathbf{F}}$ and ${\rm\mathbf{v}}$ are respectively

[TABLE]

Here, the superscripts $\mathbb{R}$ and $\mathbb{I}$ denote the real and imaginary components respectively and the column vector ${\rm\mathbf{s}}$ contains the real source pixel surface brightnesses.

The source is linearly regularised, introducing the real regularisation matrix ${\rm\mathbf{H}}$ as described in Warren & Dye (2003). The regularisation scheme we adopted follows that of Nightingale & Dye (2015), computing the mean gradient between a given Voronoi source pixel and its three nearest neighbours. To find the most probable lens model parameters, we used Markov Chain Monte Carlo (MCMC) optimisation to maximise the Bayesian evidence derived by Suyu et al. (2006). We performed multiple MCMC runs for a range of power-law density profile slopes which were kept fixed in each case to help simplify parameter space. The number of source pixels was kept fixed during optimisation and the regularisation weight was optimised following the procedure outlined in Dye et al. (2008).

3.2 Lens model

We used an elliptical power-law density profile with an external shear component where necessary to model the lenses in this work. We used the form introduced by Kassiola & Kovner (1993) which has a surface mass density, $\kappa$,

[TABLE]

where $\kappa_{0}$ is the normalisation surface mass density and $\alpha$ is the power-law index of the volume mass density profile. Here, the elliptical radius ${\tilde{r}}$ is defined by ${\tilde{r}}^{2}=x^{\prime 2}+y^{\prime 2}/\epsilon^{2}$ where $\epsilon$ is the lens elongation (i.e., the ratio of semi-major to semi-minor axis length). The orientation of the semi-major axis measured in a counter-clockwise sense from north is described by the parameter $\theta$ and the co-ordinates of the centre of the lens in the image plane are $(x_{c},y_{c})$. The external shear field is characterised by the shear strength, $\gamma$, and the shear direction angle measured counter-clockwise from north, $\theta_{\gamma}$. The shear direction angle is defined to be perpendicular to the direction of resulting image stretch. We only incorporated external shear in the lens model when the Bayesian evidence was improved by its inclusion. We found that only two of the six lenses in this work needed external shear. The total number of lens model parameters is thus eight when shear is included and six when not.

4 Results

Figure 1 shows the model reconstructions of each of the six lenses using both the image plane and visibility plane methods. It is apparent from the figure that whilst there are differences in the reconstructed sources between both methods, these are quite subtle. The variation in source plane pixelisation between image plane and uv-plane reconstructions likely accounts for a significant amount of this variation; the largest difference in morphology is seen in the case of H-ATLAS J142935.3-002836 but, owing to the random nature of the k-means clustering, this source also possesses the largest differences in source pixelisation. An anticipated tendency of the image plane method to reproduce possible artifacts arising from transformation from the visibility plane or cleaning procedure has not manifested itself in the reconstructions. Faint source features seen in each lens system are commonly reconstructed with both methods, giving an indication of their robustness. Additionally, the fact that the optimal regularisation weight may differ between the image and visibility plane due to correlated image plane pixels appears to have had little consequence222We adopted a uniform noise map for the image plane modelling, neglecting correlations between image pixels although we found that varying the pixel scale produced no significant changes in the reconstruction., although this effect may be at least partly responsible for the differences seen between some residual plots. (For example, H-ATLAS J142935.3-002836 and HELMS J001626.0+042613 show significant residuals at the location of image peaks in the image plane reconstruction compared to the uv-plane reconstruction.) The strongest features identified in the residual plots, such as those of H-ATLAS J142935.3-002836 and HELMS J004714.2+032454, have a significance of $\sim 2.5\sigma$.

Figure 1 also shows the dirty beam maps for each lens system. The strongest sidelobes occur in the HELMS beams approximately 1 arcsec east and west of the central beam component. These sidelobes each contain 6 per cent of the flux contained in the main beam component. To assess the impact that such sidelobes might have on the reconstructions, we carried out a simple test whereby we reconstructed the cleaned image of HELMS J001626.0+042613 with the dirty beam and the model beam. The resulting reconstructions showed differences in the source and model images which were only at the level of a few per cent, smaller than the differences between uv-plane and image-plane reconstructions. We therefore conclude that beam sidelobes in the current data play a negligible role.

The lens model parameters recovered for each of the six lenses using the image plane and visibility plane methods are given in table 3. On the whole, there is good agreement between the parameters obtained using the two methods, although there are mildly significant differences in the case of HELMS J001615.7+032435. However, this system has the lowest signal to noise ratio and the lack of detection of a counter image introduces additional uncertainty.

Figure 3 shows how source magnification varies as a fraction of ranked source surface brightness. We took the best fit lens model for each system (determined from the image plane modelling although the results are very similar from the uv-plane modelling – see table 4) and computed the average source magnification factor of 100 different source plane pixelisations. This was computed for different fractions of the total source flux density by working down a list of source pixels ranked by flux density (i.e. the product of source pixel area and reconstructed surface brightness). The plots show how sensitive the inferred magnification is to different interferometric configurations which probe different scales and surface brightness limits. The two systems HELMS J004723.6+015751 and HELMS J001615.7+032435 exhibit the largest variation in magnification since their sources are located in the vicinity of a caustic cusp where magnification gradients are significantly stronger.

4.1 Intrinsic source properties

We have computed intrinsic properties of the background sources in each lens system. To do this, we de-magnified the available submm photometry (see table 2) by the total source magnification factors derived from the image plane reconstructions, $\mu^{\rm img}_{\rm tot}$, as given in table 4. These are consistent with the magnifications from the uv-plane reconstructions in the sense that all differences in magnification propagate to differences in intrinsic source properties that are significantly smaller than the uncertainties arising from the SED fitting. Using the source redshifts given in table 1, we then fitted the rest-frame photometry with both a single temperature optically thick spectral energy distribution (SED) and a dual temperature optically thin SED. This SED choice gives an estimate of the upper and lower values in the range of possible dust masses, which we computed using the method outlined in Dunne et al. (2011). Here, we used the observed ALMA 880 $\mu$m flux density and a dust mass absorption coefficient computed by extrapolating the 850 $\mu$m value of $\kappa_{850}=0.077$ m2kg*-1* (James et al., 2002) to the rest-frame wavelength corresponding to the observer-frame wavelength of 880 $\mu$m (see Dunne et al., 2000, for more details). Computing dust masses in this way minimises the propagation of errors in dust temperature.

When fitting the optically thin SED, the temperature and normalisation of both components were varied. For the optically thick SED, temperature, normalisation and the opacity at 100 $\mu$m, $\tau_{100}$, were varied in the fit. In all cases, the emissivity index was fixed to 2.0 (see, for example Smith et al., 2013). The best fit SED parameters and the corresponding de-magnified luminosity of the source computed by integrating the best fit optically thin SED from 3-1100 $\mu$m are given in table 4. Finally, we computed the star formation rate of the source with the conversion from luminosity given by Kennicutt & Evans (2012) which uses a Kroupa (Kroupa, 2001) initial mass function (IMF).

4.1.1 Object notes

• H-ATLAS J142413.9+022303 - Keck K-band imaging of this system (see Calanog et al., 2014) reveals two compact galaxies interior to the Einstein ring, each consistent with an early-type morphology. Follow-up spectroscopy by Bussmann et al. (2012, B12 hereafter) gives a redshift of $z=0.595$ but due to lack of spatial resolution, it is unclear if this corresponds to solely the brighter primary galaxy or whether both galaxies have the same redshift. In this work, we have used a single power-law profile, finding that this gives a perfectly acceptable fit to the data. The lens profile centre, which is a free parameter of the fit, aligns within 0.05 arcsec of the centre of the brighter of the two galaxies. Adding a second mass to the lens model does not provide a significant improvement to the fit and makes a negligible difference to the inferred intrinsic source properties reported herein.

B12 found that a source model comprising two sersic profiles gives a significantly better fit than a single sersic profile source model. At a qualitative level, this is consistent with the irregular morphology of the reconstructed source we have obtained in the current work. B12 also estimated the de-magnified luminosity of the CO(1-0) line emitted by the source and found this to be a factor of 2.4 greater than that inferred from the line dispersion (which correlates with line luminosity; see, for example Harris et al., 2012). This discrepancy is significantly lessened to 1.4 using our magnification factor which is 80 per cent higher than that determined by B12.

The lensed source in this system has a very high star formation rate (SFR) of 2200 $M_{\odot}$/yr (see below for more discussion). This compares to the value of $\simeq 5000\,M_{\odot}$/yr reported by Bussmann et al. (2013), although this becomes $\simeq 2800\,M_{\odot}$/yr using our magnification factor instead.

• H-ATLAS J142935.3-002836 - This lens system has been previously investigated in detail by Messias et al. (2014, M14 hereafter) who analysed a broad range of multi-wavelength imaging, including ALMA Band 3 and Band 6 data (with central wavelengths of 3.1 mm and 1.3 mm respectively and maximum resolutions of 1.4 arcsec and 0.6 arcsec respectively). Optical imaging acquired with the Keck telescope (see Calanog et al., 2014) indicates that the lens is an edge-on spiral and optical spectroscopy by M14 from the Gemini-South telescope gives a lens redshift of 0.218.

The power-law lens model determined by M14 using image plane modelling of their submm/mm data has parameters $\kappa_{0}=(0.40\pm 0.01)\times 10^{10}M_{\odot}\,{\rm kpc}^{-2}$, $\alpha=2.08\pm 0.08$, $\epsilon=1.46\pm 0.04$, $\theta=136\pm 1\,{\rm deg}$ and $\theta_{E}=0.62\pm 0.08$ arcsec compared to the parameters $\kappa_{0}=(0.43\pm 0.01)\times 10^{10}M_{\odot}\,{\rm kpc}^{-2}$, $\alpha=1.79\pm 0.05$, $\epsilon=1.35\pm 0.02$, $\theta=125\pm 1\,{\rm deg}$ and $\theta_{E}=0.70\pm 0.03$ arcsec obtained directly from our much higher resolution ALMA visibility data. Whilst the models are similar, there are some significant discrepancies in certain parameters. One likely cause of this might stem from degeneracies between the triplet $\kappa_{0}$, $\alpha$ and $\epsilon$ which can give rise to substantial differences if any systematics are present (for example, arising from the fixed source plane grid used in the modelling method of M14; see Nightingale & Dye, 2015, for more details).

Our reconstructed ALMA Band 7 source has the same linear structure as that found by M14 in the submm/mm wavebands, aligned with approximately the same orientation along the lens fold caustic. Regarding the source magnification factor, our value of 24 is consistent with the values quoted in M14333In M14, magnifications were computed over different fractions of the source plane area containing 10, 50 and 100 per cent of the total source plane flux. M14 computed a 50 per cent magnification of 14 and a 10 per cent magnification of 26. To be consistent with the definition used by M14 would require a source plane fraction somewhere between these two values.. In our reconstruction, there is a hint of morphological disturbance at the southern end of the source. This is exactly where M14 find that a second optically detected source intersects in what they interpret as a possible merger.

This source has an extremely high SFR to dust mass ratio, the highest in our sample. The source lies $>3\sigma$ away from the mean in the distribution of SFR to dust mass ratios of high redshift submm galaxies (SMGs) and lower redshift ULIRGs determined by Rowlands et al. (2014) as Figure 5 shows.

• HELMS J004714.2+032454 - This is a double image system which is very well fit with a single power-law density profile and no external shear. The source exhibits a long faint structure extending to the south-east and this is readily seen in the lensed image.

The SPIRE and ALMA photometry alone continues to rise towards shorter wavelengths, the peak of the SED being constrained purely by the PACS photometry. The relatively high 100$\,\mu$m PACS flux is suggestive of a warmer dust component and this is reflected in a significantly better fit by the dual temperature SED compared to the single temperature template, although both SEDs give a comparable dust mass.

De-magnifying the far-IR luminosity given in N16 using our magnification factor of 8.3 gives $\log(L_{\rm FIR}/L_{\odot})=12.1\pm 0.1$, slightly less than our determination but consistent within the uncertainties. The luminosity implies a star formation rate of $\simeq 220\pm 60$ $M_{\odot}$/yr. Given its dust mass range of $10^{8.7}-10^{9.2}$ $M_{\odot}$, this places the source somewhere between having the characteristics of a high redshift SMG or lower redshift ULIRG and the bulk population of $z<0.5$ galaxies detected in H-ATLAS, according to Rowlands et al. (2014).

• HELMS J001626.0+042613 - This double image system is well described by an isolated power-law density profile and a relatively compact source. Both reconstruction methods suggest faint extended source structure but this only contributes a few per cent of the main source flux. The system has the lowest magnification factor in our sample of only $4.1\pm 0.3$.

The peak of the source SED in this system is well bounded by the ALMA and SPIRE photometry giving robust temperature estimates. In the dual temperature SED, the warm component makes a larger contribution to the total dust mass than the other five sources but this is not well constrained owing to uncertainties in the shorter wavelength PACS photometry. The de-magnified source luminosity is $\log(L_{\rm FIR}/L_{\odot})=12.7\pm 0.1$ which agrees with the value quoted by N16. The $z=2.51$ source has a high SFR of $980$ $M_{\odot}$/yr and its SFR to dust mass ratio is consistent with a typical SMG/ULIRG as indicated in Figure 5.

• HELMS J004723.6+015751 - This system is one of two in our sample which require external shear in the lens model, consistent with the location of a smaller external galaxy 10 arcsec to the south. The source shows a compact, relatively featureless morphology with the hint of an extended structure to the north west.

The SPIRE and ALMA photometry of the source on their own indicate that the peak of the SED lies in the vicinity of the shortest wavelength data point at 250 $\mu$m. This is borne out by the inclusion of PACS photometry. As a result, the fitted dual temperature SED implies a dominant mass of cold dust at 26 K. The intrinsic source luminosity of $\log(L_{\rm FIR}/L_{\odot})=12.2$ is in agreement with that measured by N16. The SFR of $230\pm 60$ $M_{\odot}$/yr for this $z=1.44$ source compared with its relatively low dust mass places it in the upper envelope of SFR to dust mass ratios spanned by SMGs and ULIRGs according to Rowlands et al. (2014).

• HELMS J001615.7+032435 - The relatively low image signal-to-noise ratio in this cusp-caustic configuration lens results in an undetected counter-image which increases the modelling uncertainty for this system. Nevertheless, the most probable lens model is one with a significant external shear. This is consistent with several smaller nearby galaxies, mainly to the north-east, with colours similar to the lens which is, in turn, consistent with the larger Einstein radius of a group-scale lens.

In light of this, we attempted a lens model that includes external convergence provided by a singular isothermal ellipsoid (SIE) mass model. The best fit model we found places the SIE to the north-east with the result that the required external shear is reduced by approximately 30 per cent and the normalisation of the primary lens, $\kappa_{0}$, is lowered by approximately 20 per cent. The magnification is also reduced by approximately 30 per cent. However, the model is less favoured by the Bayesian evidence and there is a tendency for it to produce a brighter counter image which would have been detected in the ALMA data. The location and normalisation of the external SIE is, as expected, degenerate with the normalisation and shear of the primary lens. Further observations of the lensing galaxies are required to better characterise the lens model.

The ALMA and SPIRE photometry of the source in this lens system prefers an optically thick single temperature SED. However, with the inclusion of PACS fluxes, a marginally improved fit is obtained with a second weak but quite hot dust component, although the improvement in the fit is not significant given the additional SED parameters. The source has a luminosity of $\log(L_{\rm FIR}/L_{\odot})=12.5$ which agrees with that of N16 who used an optically thin single component SED. The SFR to dust mass ratio of this source is extremely high, placing it nearly $3\,\sigma$ above the mean in the distribution of ratios measured in the SMG/ULIRG population.

5 Summary and Discussion

We have modelled ALMA imaging data of six strong galaxy-galaxy gravitational lens systems originally detected by the Herschel Space Observatory. For each lens system, we have carried out modelling of both the cleaned image data and the visibility data directly. We find only minor differences in the reconstructed source morphologies between the two methods. The expectation is that such differences will become more prominent as coverage of the uv-plane becomes more sparse, not least because this will generally lead to larger-scale image pixel covariances from beam sidelobes which are not included in the cleaned data. In Dye et al. (2015), modelling of the cleaned image was advocated on the basis that the uv-plane was very well sampled in that particular case and because image plane modelling is substantially more computationally efficient than uv-plane modelling generally. In the present work, the uv-plane is less well sampled in comparison and hence the decrease in efficiency by modelling the visibility data is less severe. Nevertheless, image plane modelling is still at least an order of magnitude quicker than uv-plane modelling and gives very similar results.

In our fitting of a smooth power-law mass density profile, we have found that the lenses are all close to isothermal and that the recovered model parameters are in broad agreement between both methods. However, one system with particularly poor signal to noise shows mildly significant discrepancies in the slope and normalisation of the power-law profile, although these two parameters are typically quite degenerate. A more exhaustive investigation into the origin, prevalence and strength of such discrepancies along with differences in the reconstructed source is left for future study.

We have used the lens magnification factors obtained from the modelling to demagnify the submm source photometry. Fitting rest-frame SEDs to this photometry, we have determined the dust temperature, dust mass, luminosity and inferred star formation rate of the lensed sources. Using both an optically thick single-temperature SED and an optically thin SED with two temperature components has allowed an estimate of the range of dust mass possible for each source. Taking the mid-point of this range in each case, we find that five of the six sources have a ratio of star formation rate to dust mass which is in excess of the mean ratio of the SMG/ULIRG population as determined by Rowlands et al. (2014).

The extent of this excess is shown in Figure 5 which plots the SFR obtained by scaling the far-IR luminosity using the relation given by Kennicutt & Evans (2012) against dust mass. The figure shows that two of the sources in our sample are at least as extreme as the H-ATLAS lensed source SDP.81 investigated by Dye et al. (2015). These lie in the upper envelope of the distribution of SFR-to-dust mass measured by Rowlands et al. (2014). Since our computed SFR is simply a scaled version of far-IR luminosity, the underlying fact is that these sources have a high luminosity for the quantity of gas available for star formation. This is often an indication that a component of the source’s luminosity comes from an active galactic nucleus but we are unable to comment further on this possibility without additional observations.

If we convert the rest-frame 850$\mu$m flux density of our sources to H2 gas mass (see table 4) using the empirical scaling relation given by Hughes et al. (2017), we find that the five sources located above the Rowlands et al. SFR-to-dust mass relationship also lie on or above the mean relationship between SFR and H2 gas mass determined by Scoville et al. (2016). If dust is indeed an accurate tracer of molecular gas as these scaling relationships suggest, then the implication is that these sources possess a higher star formation efficiency (SFE). Treating the range in dust mass for each source as a 1-sigma error and fitting a line parallel to the SMG/ULIRG relationship in Figure 5 to the mid-point of the dust mass range for all six sources, the increase in SFE is a factor of 5 relative to that implied by the SMG/ULIRG relationship of Rowlands et al. and a factor of 40 relative to $z<0.5$ H-ATLAS galaxies.

An alternative explanation to the SFR-to-dust mass offset being the result of an enhanced SFE could be that the gas-to-dust ratio in these sources is higher. Similarly, the results would be explained if the dust mass opacity coefficient were lower by the factors mentioned above. Both of these possibilities seem to disagree with measurements of gas mass from CO detections at low and high redshift (see, for example, Dunne & Eales, 2001; Magdis et al., 2012; Rowlands et al., 2014; Scoville et al., 2014, 2016; Grossi et al., 2016; Hughes et al., 2017). These studies indicate a tight correlation between CO line intensity and 850$\,\mu$m luminosity, thereby implying a constant H2 gas-to-dust mass ratio. However, a caveat is that this assumes a fixed value of the ratio of H2 surface gas mass density to CO line intensity, $\alpha_{\rm CO}$. Sandstrom et al. (2013) find a weak dependence of $\alpha_{\rm CO}$ on metallicity in local galaxies, such that lower metallicity tends to correspond to higher values of $\alpha_{\rm CO}$. If this holds in high redshift SMGs, whilst a lower metallicity would not affect the CO-to-dust ratio, the ratio of H2 gas-to-dust would be increased, leading to an enhanced SFR-to-dust mass ratio.

An additional point to note is that interpreting a higher SFR to gas mass ratio as a higher SFE when the total molecular gas mass is used assumes that star formation occurs throughout the full extent of molecular gas. Determinations of dense molecular gas mass traced by HCN emission show a correlation between far-IR luminosity and HCN line intensity that is much tighter than the correlation between HCN and CO line intensity (see for example, Gao & Solomon, 2004; Privon et al., 2015). SFR therefore appears to depend on dense molecular gas mass rather than total molecular gas mass traced by CO. In light of this, and assuming universal star formation physics, a more probable interpretation of the high SFR to gas mass ratios we find is that the sources in our sample have a significantly higher dense molecular gas mass fraction. This conclusion was also reached by Oteo et al. (2017) who carried out a similar analysis of two H-ATLAS lensed sources.

Papadopoulos & Geach (2012) provide evidence to suggest that high density molecular gas is more prevalent in galaxy mergers than quiescently forming systems and that its fraction can be used to determine the mode of star formation. Inspection of the reconstructed morphologies (Figure 1) of the two sources in our sample with extreme SFR to gas mass ratios (i.e., HELMS J001615.7+032435 and H-ATLAS J142935.3-002836) does indeed reveal signs of disturbed morphology, but no more so than others in the sample. Nevertheless, increasing the number of gravitational lens reconstructions of such systems with high magnification factors offers the ability to further investigate such hypotheses. This becomes especially true with the inclusion of source kinematics measured via molecular lines.

Acknowledgements

SD acknowledges support from the UK STFC Ernest Rutherford Fellowship scheme. LD acknowledges funding from the European Research Council Advanced Investigator grant COSMICISM and the ERC Consolidator grant CosmicDust. MN acknowledges financial support from the European Union’s Horizon 2020 research and innovation programme under the Marie Skłodowska-Curie grant agreement No 707601. MJM acknowledges the support of the National Science Centre, Poland through the POLONEZ grant 2015/19/P/ST9/04010. LM acknowledges support from the South African National Research Foundation through the South African Research Chairs Initiative. This project has received funding from the European Union’s Horizon 2020 research and innovation programme under the Marie Skłodowska-Curie grant agreement No. 665778. Some of the spectroscopic redshifts reported in this paper were obtained with the Southern African Large Telescope (SALT) under proposal 2015-2-MLT-006. This paper makes use of the following ALMA data: ADS/JAO.ALMA#2013.1.00358.S. ALMA is a partnership of ESO (representing its member states), NSF (USA) and NINS (Japan), together with NRC (Canada) and NSC and ASIAA (Taiwan) and KASI (Republic of Korea), in cooperation with the Republic of Chile. The Joint ALMA Observatory is operated by ESO, AUI/NRAO and NAOJ. This research has made use of the NASA/IPAC Infrared Science Archive, which is operated by the Jet Propulsion Laboratory, California Institute of Technology, under contract with the National Aeronautics and Space Administration.