Modelling high-resolution ALMA observations of strongly lensed dustystar forming galaxies detected by Herschel

We present modelling of ~0.1arcsec resolution Atacama Large Millimetre/sub-millimeter Array imaging of seven strong gravitationally lensed galaxies detected by the Herschel Space Observatory. Four of these systems are galaxy-galaxy scale strong lenses, with the remaining three being group-scale lenses. Through careful modelling of visibilities, we infer the mass profiles of the lensing galaxies and by determining the magnification factors, we investigate the intrinsic properties and morphologies of the lensed sub-millimetre sources. We find that these sub-millimetre sources all have ratios of star formation rate to dust mass that is consistent with or in excess of the mean ratio for high-redshift sub-millimetre galaxies and low redshift ultra-luminous infrared galaxies. The contribution to the infrared luminosity from possible AGN is not quantified and so could be biasing our star formation rates to higher values. The majority of our lens models have mass density slopes close to isothermal, but some systems show significant differences.


INTRODUCTION
Sub-millimetre (sub-mm) galaxies (SMGs) (see Casey et al. 2014, for a review) play host to some of the most intense star formation rates in the Universe. They are observed to be abundant at redshifts > 1 and to contribute approximately 20 per cent of the cosmic star formation rate density up to a redshift of ∼ 4 (Swinbank et al. 2013;Lapi et al. 2011). The rest-frame UV and optical radiation associated with young stars is highly obscured by the dust content of sub-mm galaxies ★ Herschel is an ESA space observatory with science instruments provided by European-led Principal Investigator consortia and with important participation from NASA † E-mail: jacob.maresca@nottingham.ac.uk (Dudzevičiūtė et al. 2020). Dust grains absorb this radiation, causing them to be heated. This energy is reprocessed and emitted as thermal continuum emission that we can observe in the sub-mm and mm regimes. Sub-mm galaxies inhabit a key role in our picture of galaxy evolution; massive SMGs at high-redshift evolve through the population of quiescent galaxies observed at lower redshift (Simpson et al. 2014;Toft et al. 2014), and through the mechanism of gas-poor minor mergers, helps to explain the build up of massive elliptical galaxies observed in the present day Universe (Oogi & Habe 2012;Guo & White 2008;Lapi et al. 2018). In the local Universe, Ultra Luminous Infrared Galaxies (ULIRGs), are often seen as analogues to the high-redshift sub-mm galaxies due to their strongly dust-obscured UV luminosity and high infrared luminosity, and whilst they are considerably less abundant, they have comparable bolometric lumi-nosities (Alaghband-Zadeh et al. 2012;Rowlands et al. 2014). The study of high-redshift sub-mm galaxies has benefited greatly from the advent of large interferometer arrays such as the Atacalma Large Millimeter/Sub-millimeter Array (ALMA), allowing observations to reach resolutions of < 0.1 arcsec, probing physical scales that were previously unreachable. Strong gravitational lensing provides a further boost in spatial resolution due to the magnification of the background source, which is typically within the range of 5-10 for SMGs (Spilker et al. 2016). In addition to this, there exists a strong lensing bias in the sub-mm regime, making it possible to find lensed sources in wide surveys with a simple cut in flux density above 100 mJy at 500 m (Blain 1996;Negrello et al. 2007Negrello et al. , 2010Perrotta et al. 2003). Using this technique, follow-up ALMA observations of strongly lensed sub-mm galaxies detected in wide area surveys, such as the Herschel Astrophysical Terahertz Large Area Survey (H-ATLAS; Eales et al. 2010), the Herschel Multi-tiered Extragalactic Survey (HerMES; Oliver et al. 2012) and the Herschel String 82 Survey (HerS; Viero et al. 2014), carried out using the Herschel Space Observatory (Pilbratt, G. L. et al. 2010), as well as the sub-mm surveys conducted by the Planck satellite (Cañameras et al. 2015) and the millimetre surveys of the South Pole Telescope (Carlstrom et al. 2011) have contributed to a rapidly increasing understanding of galaxy formation at its early stages Sun et al. 2021;Harrington et al. 2021;Cañameras et al. 2017a,b). Many earlier observations were focused on extremely luminous sources, but with the increased sensitivity in instruments such as ALMA, it is becoming possible to investigate more typical main-sequence star-forming galaxies.
High-resolution sub-mm follow-up observations allow for precision lensing analyses, greatly benefiting our characterisation of source properties sensitive to the lens model, such as the intrinsic luminosity, star formation rate and gas and dust masses. Studies of the molecular gas and dust in lensed sub-mm galaxies allow us to test models of star formation in the early Universe (Cava et al. 2018;Dessauges-Zavadsky et al. 2019). Recent studies have begun to reveal the compact nature of dust in sub-mm galaxies (Puglisi et al. 2019;Tadaki et al. 2020) and the disparity in size of these regions when compared with local ULIRGs. Similarities in size, number density and clustering properties with compact quiescent galaxies is suggestive of an evolutionary connection Dudzevičiūtė et al. 2020), and thus understanding of quenching in sub-mm galaxies is important to understand how they may become red and dead. The advent of observatories such as ALMA have revolutionised our understanding of high-redshift lensed sources when compared to what could be achieved with optical imaging. A classic example of this is the lensed system SDP.81, first discovered within the H-ATLAS sample and then observed with ALMA Rybak et al. 2015a;Rybak et al. 2015b;Negrello et al. 2010).
As our sample size of known strongly lensed sub-mm galaxies increases, so does our range of redshifts at which they have been observed (Wang et al. 2007;Riechers et al. 2013) thanks to a very negative K-correction such that high-redshift galaxies have approximately constant brightness in the sub-mm regime between redshifts ∼ 1 − 8. Higher redshift sources tend to be lensed by higher redshift lenses, due to the scaling of the lensing cross-section with redshift Turner et al. (1984). Increasing the redshift range then allows us to study the mass profiles of lenses at earlier epochs, when galaxy evolution is more rapid and less well understood (Dye et al. 2014;Negrello et al. 2014).
Reconstruction of the background lensed source from interferometer observations can be achieved with two distinct approaches. There are methods that model the visibilities directly in the uv-plane (Bussmann et al. 2012;Bussmann et al. 2013;Dye et al. 2018), and those that model the cleaned image plane data Inoue et al. 2016;Yang et al. 2019). The benefit of working in the image plane is that the task is significantly less computationally expensive, however, due to the incomplete coverage of the uv-plane, spatially correlated noise is introduced which can in principle bias the inferred lens model. Working directly with the visibility data avoids this problem at the cost of longer modelling times.
In this work, we have carried out lens modelling in the uv-plane of seven ALMA observations, four of which are galaxy-galaxy scale, and the remaining three are group-scale lenses. These systems were originally detected by Herschel within H-ATLAS and the extension to the HerMES field, HerMES Large Mode Survey (HeLMS) (Asboth et al. 2016;Nayyeri et al. 2016). We have investigated the intrinsic source properties, such as luminosity and SFR, by determining the magnification factors. Additionally, we investigate the morphologies of the reconstructed sources.
The layout of this paper is as follows: Section 2 describes the ALMA observations and other sources that were drawn upon for this work. Section 3 details the methodology of our lens modelling and Section 4 presents the results of our work. In Section 5 we compare our results to other similar studies. Finally, in Section 6 we summarise our findings and discuss their interpretation. Throughout this paper we assume a flat ΛCDM cosmology using the 2015 Planck results (Planck Collaboration et al. 2016), with Hubble parameter ℎ = 0.677 and matter density parameter Ω = 0.307.

DATA
The seven ALMA observations modelled in this work are from the ALMA programme 2013.1.00358.S (PI: Stephen Eales) and are described in detail within Amvrosiadis et al. (2018). The observation targets for the original ALMA programme were selected from the H-ATLAS and HeLMS surveys for having the brightest 500 m flux densities and with spectroscopic redshifts > 1. Of the 16 sources observed during ALMA cycle 2, 14 display obvious lensing features. Six of those sources are modelled in the work by Dye et al. (2018), one of them we leave for future work, and the remaining seven are modelled here. Of these seven sources, one was identified in H-ATLAS whilst the remaning six are from HeLMS.
The spectral setup employed by ALMA was identical for each of the lensing systems observed. The band 7 continuum observations were comprised of four spectral windows, each with a width of 1875 MHz and centred on the frequencies 336.5 GHz, 338.5 GHz, 348.5 GHz and 350.5 GHz. The central frequency of 343.404 GHz corresponds to a wavelength of 873 m. Each spectral window is comprised of 128 frequency channels, resulting in a spectral resolution of 15.6 MHz. The ALMA configuration utilised forty two 12m antennae with an on-source integration time of approximately 125 s. Upon combining all four spectral windows, this achieves an angular resolution of 0.12 arcsec and RMS values of approximately 230 Jy/beam and 130 Jy/beam for the H-ATLAS and HeLMS sources, respectively. The synthesised beam shape for HeLMS J005159.4+062240 is the least elliptical of our sample, with a ratio of major to minor axis of ∼ 1.2. For the observation of HeLMS J235331.9+031718, the ratio of major to minor axis of the synthesised beam is the most elliptical of our sample at ∼ 1.7, with the remaining observations having ratios of ∼ 1.5. The maximum recoverable scales for our observations are between 1.32 arcsec and 1.46 arcsec.
In this work, we have used the visibility data provided by the ALMA science archive, and re-calibrated them using C Table 1. The list of the seven lensing systems modelled in this work, along with their lens galaxy redshifts, , and their background source redshifts . Where Appropriate, the redshift of both lensing galaxies have been provided, distinguished by the numbered subscript. The references from which the lens and source redshifts were obtained are as follows: a Bussmann et al. (2013), b Negrello et al. (2017), c Nayyeri et al. (2016), and d Okido et al. (2020).
The range of redshifts given for HeLMS J235331.9+031718 is based on the range of source redshifts in this paper, since there is no available redshift measurement for this source. A dash indicates missing redshift information for a lens whilst 'N/A' is used to indicate a second mass profile was not used in our modelling procedure.  -Mullin et al. 2007) and the scripts provided by the archive. Baselines flagged as bad by the ALMA data reduction pipeline were excluded from the analysis. The CASA task tclean was used to create images in order to measure the flux of the sources at 873 m. The images were constructed using a natural weighting scheme and were primary beam corrected. To well sample the minor axis of the primary beam, an image pixel scale of 0.02 arcsec and 0.03 arcsec was used for the H-ATLAS and HeLMS sources respectively. For the creation of these images, we used CASA version 6.1.2.7.

ID
For the calculation of intrinsic source properties, photometry from our ALMA data was used in combination with a number of other datasets. Sub-mm photometry, obtained by the Herschel Space Observatory, making use of both the Spectral and Photometric Imaging Receiver (SPIRE; Griffin et al. 2010)  The available lens and source redshifts for the seven systems modelled in this paper can be found in Table 1. The observed photometry can be found in Table 2.

The semi-linear inversion method in the uv-plane
The standard image plane approach of the semi-linear inversion (SLI) method makes use of a pixelised source plane. For a given lens model, the image of each pixel is formed and the linear superposition of these images that best fits the data determines the source surface brightness distribution. Analogously to the image plane version, when working with interferometer visibility data, a model set of visibilities is formed for an image of each source pixel. The linear combination of these model visibilities determines the source surface brightness distribution for a particular lens model.
We used the source inversion method implemented within PyAutoLens (Nightingale et al. 2021) which is based on the operator approach described within Powell et al. (2020). An interferometer visibility dataset is comprised of samples of complex visibilities. The surface brightness in the source plane is given by the vector , with each element corresponding to the surface brightness of a source plane pixel. The parameterised projected surface mass density of the lens model is given by the vector and the mapping of the source light to the image plane is described by the operator L( ). The sky brightness is therefore simply L( )s. The response of an interferometer is encoded into the operator D, which performs the Fourier transforms to convert the pixelised sky brightness distribution into a set of complex visibilities. The observed data , can therefore be described by the combination of these effects = DL( )s + n. (1) Assuming uncorrelated noise in the visibility data, the diagonal matrix C −1 represents the covariance. Using a similar method to Dye et al. (2018) to determine the 1 uncertainties on the visibilities, we used the CASA task statwt to empirically measure the visibility scatter computed over all baselines.
Combining this with the set of model visibilities DL( ) allows us to write the 2 statistic as With the addition of a prior on the source denoted by the operator R and with the regularisation strength given by , Powell et al. (2020) show that we can write the regularised least-squares equation as where the maximum a posteriori source inversion matrix, i.e the solution matrix for MP is given by the quantity in square brackets in equation 3. This linear system of equations is in principle straightforward to solve but becomes extremely memory intensive for large numbers of visibilities and/or large numbers of source plane pixels. For this reason, a Direct Fourier Transform (DFT) is replaced by a non-uniform Fast Fourier Transform (NUFFT), constructed out of a series of operators. In PyAutoLens, the Fourier transforms are performed using the NUFFT algorithm PyNUFFT (Lin 2018), which greatly increases computational efficiency over the DFT method. This substitution results in a modified version of equation 3, where the solution for MP is given by a series of operators evaluated by use of an iterative linear solver. The linear algebra package PyLops (Ravasi & Vasconcelos 2020), is used to achieve this (see Powell et al. 2020;Nightingale et al. 2021, for more details on this methodology and the specific implementation used here).

Lens Model
We have used the elliptical power-law density profile, which is a generalised form of the singular isothermal ellipsoid that is commonly used to fit strong lens profiles (Keeton 2001). When it improves the Bayesian evidence, an external shear component is included to compensate for the influence of line of sight galaxies that may be outside our field of view. Where necessary, two elliptical power-law profiles have been used to model the group-scale lenses present in the sample. We find that in all cases this is sufficient to provide acceptable fits to the data and through the use of optical/near-IR imaging, we see no indication of more than two lenses being required to produce an accurate model. The surface mass density, , of this profile is given by where is the model Einstein radius in arc seconds, the powerlaw index and the elliptical radius defined as = √︁ 2 + 2 / 2 , where is the axis ratio (Suyu 2012). The lens orientation , is measured as the angle counter-clockwise from the east axis to the semi-major axis of the elliptical lens profile. The centre of the lens profile is given by the coordinates in the image-plane ( , ). The external shear field is parameterised by a shear strength and shear direction , measured counter-clockwise from east. The shear direction is defined to be perpendicular to the direction of the resulting image stretch. For single lens systems there are six lens model parameters (eight when including external shear) and 12 lens model parameters (14 when including external shear) in the case of groupscale lenses.

RESULTS
Fig 1 shows the model image, residual image and source reconstruction for each of the seven lenses modelled in this work. The lens model parameters from our fitting procedure are given in Table 3.
Different interferometer configurations are able to probe distinct scales and reach varying surface brightness limits. Fig 2 shows how for each system the inferred magnification is sensitive to this effect. Working down a list of source plane pixels, above a surface brightness threshold, ranked by flux density (i.e the product of their reconstructed surface brightness and area), the reconstructed image and average magnification was computed. The source flux fraction refers to the fraction of the total flux contained by a subset of these ranked pixels, e.g five pixels with the greatest flux might contain one third of the total flux of all the pixels in the source plane above a surface brightness threshold. Likewise, the image flux fraction refers to the fraction of the total flux in the image plane that these subsets of source plane pixels contribute when they are lensed by our bestfit model. This process was repeated 100 times using a randomised source plane pixelisation for each, to produce an averaged magnification profile. The total magnification is computed by calculating the flux in the source plane (for pixels above a surface brightness threshold) and the flux in the image plane due to the lensed image of the source plane pixels. Only the flux lying within the image plane mask contributes to the image plane flux and thus the total magnification.

Intrinsic source properties
For each lens system we have determined the intrinsic properties of the background source. To achieve this, we have demagnified the sub-mm photometry (see Table 2) by the total source magnification factors given in Table 4, taking into account the uncertainties on our magnification values. With the source redshifts given in Table  1, we fitted the rest-frame photometry with two Spectral Energy Distributions (SEDs). Firstly, a single temperature optically thick SED of the form where is the flux density at frequency , 0 is the frequency at which the optical depth equals unity, is the dust emissivity index, is the dust temperature and ( , ) denotes the Planck function. Secondly, a dual temperature optically thin SED was fitted of the form where is the weighting of the cold and warm components (denoted with subscripts) and is the dust temperature of the two components (denoted by the subscripts). This allows for the computation of an upper and lower bound in the range of possible dust masses, which were determined using the method described in Dunne et al. (2011).
In this work we have used the 873 m ALMA flux density and computed the dust mass absorption coefficient by extrapolating the 850 m value of 850 = 0.077m 2 kg −1 (James et al. 2002) (see Dunne et al. (2000) for more details). The uncertainties in the dust mass absorption coefficient are known to be large, but a relative comparison between dust masses computed with this value is still valuable. The same value of 873 was used for all of our galaxies, which assumes that the physical properties of the dust, such as grain size and density, are constant for our sample of galaxies. The scaling   , coordinates ( , ) of the centroid of the mass profile with respect to the phasetracking centre of the observations ( and correspond to west and north respectively), the lens profile orientation measured as the angle counter clockwise from East to the semi major axis, , the lens profile axis ratio, , the density slope of the power-law, , the magnitude of the external shear field, , and the external shear field direction, , measured counter clockwise from east   Table 4. Intrinsic source properties. The columns are total magnification, , dust mass computed from the single temperature optically thick SED, ℎ , dust mass assuming a dual temperature optically thin SED, ℎ , temperature of the optically thick SED, ℎ , temperatures of the optically thin SED, ℎ , the optical depth at 100 m for the optically thick SED, 100 , the demagnified luminosity (computed as the integral of the optically thin SED from 3 to 1100 m), , the 2 gas mass computed using the scaling relation of Hughes et al. (2017), , and the SFR scaled from using the procedure given by Kennicutt & Evans (2012) with a Kroupa IMF. Dust masses are expressed as log( d / ), luminosity as log( IR / ), and gas masses as log( gas / ).
The source HeLMS J235331.9+031718 appears twice, displaying the intrinsic source properties for the upper and lower value of source redshift being considered (indicated in the ID column).  relations from Hughes et al. (2017) were used to calculate the H 2 gas mass. During the fitting of the optically thin SED, the temperature and normalisation of both dust components were allowed to vary. In the case of the optically thick SED, the temperature, normalisation and opacity at 100 m, 100 , were varied during the fit. Throughout, the emissivity index, , was fixed to 2.0 (see Smith et al. 2013). The fitted SEDs and the demagnified source photometry can be seen in Fig 3. The best fit parameters for these SEDs can be found in Table  4, along with the demagnified luminosity of the sources, computed by integrating the optically thin SED from 3-1100 m. Finally, the star formation rates of the sources are given, computed using the conversion from luminosity log(SFR) = log(L IR ) − 43.41 (7) given by Kennicutt & Evans (2012) which uses a Kroupa (Kroupa 2001) Initial Mass Function (IMF). This conversion assumes that all of the infrared emission originates from star forming regions, which may result in biased SFR values if an AGN is significantly contributing to the luminosity of the galaxy.

Object notes
H-ATLAS J083051.0+013225 -This lens system has an almost complete ≈ 1.5 arcsec Einstein ring with three major image components, along with a fainter central image. The central image is due to the rare line of sight configuration of the dual deflectors. Keck Adaptive Optics band imaging of this lensing system (see Calanog et al. 2014) shows two galaxies, but the lensing features are ambiguous due to the low signal to noise ratio. Upon superimposing the lensed emission from the background source as detected in our ALMA data, it is clear that both of these galaxies are interior to the Einstein ring. Long slit spectroscopy presented in Bussmann et al. (2013) provides evidence that these two galaxies are at different redshifts (0.626 and 1.002 respectively). Bussmann et al. (2013) presents the modelling of SMA data for this system, fitting two SIE models for the foreground deflectors and a Sérsic model for the background source. They infer a magnification factor = 6.9 ± 0.6, which is in agreement with the magnification we find with our best fit model of = 6.7 ± 0.5. Yang et al. (2019) presents ALMA band 4 data, along with the band 7 data presented in this paper, and finds a significantly higher magnification factor of = 10.5 +0.5 −0.6 , using a double SIE lens model and a dual disk model for the background source.
The star formation rate of this lensed source is extremely high  Table 4. The SED for HeLMS J235331.7+031718 is shown at a redshift of = 2. at 3700 ± 500 yr −1 , in reasonable agreement with the value reported by Zhang et al. (2018), who used the magnification factor from Bussmann et al. (2013) to determine the SFR. This compares to the lower values of SFR inferred by Yang et al. (2019) who find 600 ± 300 yr −1 and 900 ± 400 yr −1 for each of the components in their source model. Our value of star formation rate is reduced to 2400 yr −1 when using the magnification factor of 10.5 found by Yang et al. (2019). A possible explanation for the large discrepancy in SFR is differential magnification of the source. Yang et al. (2019) find evidence of differential magnification for the compact and extended components of their source model, which they have taken into account in their SED modelling. In addition to this, whilst the lens model that we find is similar to that of Yang et al. (2019), we have allowed the power-law slope to vary rather than using SIE profiles. The higher values (>2) of slope that we find will tend to decrease the overall magnification of the source plane.
Our reconstructed source for this system shows significant disturbance, with a main component oriented North-South and a fainter component oriented East-West.
HeLMS J005159.4+062240 -This system is a group-scale lens with a doubly-imaged background source. Optical imaging from the Sloan Digital Sky Survey (SDSS) Data Release 14 (Abolfathi et al. 2018) reveals two lensing galaxies within the lensed arcs. Tergolina (2020) investigated these objects, finding that one was a red and passive early-type galaxy (ETG) ( = 0.60246 ± 0.00004), whilst the other is the host of a quasar ( = 0.59945 ± 0.00009) (Okido et al. 2020). We were able to find a lens model constructed from two elliptical power-law density profiles and an external shear field that reconstructs the background source as a single component. Our reconstructed source appears elongated with a disturbed morphology.
With a magnification factor of = 10.2 ± 0.3, this system has the highest total magnification in our sample. Our measurement of the intrinsic source luminosity, log( IR / ) = 12.6 ± 0.1, is in reasonable agreement with the far-IR luminosity given in Nayyeri et al. (2016), once the lens magnification has been taken into account. The peak of the SED is bounded by the SPIRE photometry and the PACS flux densities indicate the presence of a warm dust component. As such, the two temperature SED provides a better fit to the data. With a star formation rate of ≈ 590 ± 50 yr −1 and a dust mass range of 10 8.2 − 10 8.9 this source's SFR to dust mass ratio is typical of high-redshift SMGs and low redshift ULIRGs according to the empirical relations in Rowlands et al. (2014).
HeLMS J234051.5-041938 -This system displays a nearlycomplete ≈ 0.75 arcsec Einstein ring with three image components. A singular power-law with external shear provides a good fit to the data, with the power-law index preferring relatively high values of around 2.3. The source has evidence of a disturbed morphology, displaying extended faint emission to the west and east.
The peak of the source SED is constrained by our ALMA and SPIRE photometry. The PACS flux densities do not indicate the presence of a significant warm dust component, with both the single and double temperature SEDs providing excellent fits to the data. This source's intrinsic luminosity of log( IR / ) = 13.1 ± 0.1 agrees well with the luminosity calculated by Nayyeri et al. (2016). A star formation rate of ≈ 1900 ± 300 yr −1 and a dust mass range of 10 8.3 − 10 9.1 places this source above the typical dust mass to SFR ratio for high-redshift SMGs and local ULIRGs.
HeLMS J232439.5-043935 -This system exhibits an unusual image configuration with three distinct image components, the two western most of which show faint extended features, whilst the eastern component is more compact. Two power-law density profiles with an external shear field provide an excellent fit to the data, reconstructing a source with faint extended emission to the north.
The peak in flux density for this source occurs close to the 250 m SPIRE measurement, with the constraints coming from the PACS measurements. The SED does not show clear evidence of a significant second temperature component, with both SED models providing good fits to the measured flux densities.
With an intrinsic luminosity of log( IR / ) = 13.0 ± 0.1, our estimate is consistent within our uncertainties of the value given by Nayyeri et al. (2016). A star formation rate of ≈ 1500 ± 300 yr −1 and a dust mass range of 10 8.5 − 10 9.1 means that this source has a higher than typical ratio of SFR to dust mass for high-redshift SMGs and local ULIRGs.
HeLMS J233255.4-031134 -A doubly-imaged source with some faint extended emission emanating from the southern image, that is well fit by a single power-law density profile and an external shear field. The reconstructed source is relatively compact and featureless.
The peak of the SED is constrained by the SPIRE photometry, with the PACS flux densities not showing any major warm dust component. The SED is well described by an optically thick model. The intrinsic source luminosity we have obtained of log( IR / ) = 13.1±0.1 is in agreement with Nayyeri et al. (2016). The SFR of ≈ 2000 ± 300 yr −1 for this source agrees with the value quoted by Zhang et al. (2018) within our stated uncertainties. A dust mass range of 10 8.5 − 10 9.1 means that this galaxy also lies above the mean SFR to dust mass ratio as given by Rowlands et al. (2014).
HeLMS J233255.6-053426 -This quadruple image system is well described by a power-law density profile embedded in an external shear field. The source exhibits a relatively featureless compact morphology. The peak of the source SED is well constrained by the ALMA and SPIRE photometry. The PACS photometry indicates the presence of a warmer dust component, as shown by the relatively high 100 m PACS flux density measurement, and the significantly better fit of the two temperature SED. Our measurement of the intrinsic source luminosity, log( IR / ) = 12.5 ± 0.1, is in agreement with the far-IR luminosity given in Nayyeri et al. (2016) using our magnification factor of 9.2 to de-magnify the quoted value. With a star formation rate of ≈ 500 ± 50 yr −1 and a dust mass range of 10 8.5 − 10 9.0 this source's SFR to dust mass ratio is consistent with typical high-redshift SMGs and low redshift ULIRGs as indicated by Fig 4. HeLMS J235331.9+031718 -This double image system has an extremely small image separation, with an Einstein radius of ≈ 0.1 arcsec. This system is described by a single power-law density profile and an external shear field. The inferred slope of the power-law density profile is relatively low at = 1.64 ± 0.04, contributing to the high magnification of this source, which also lies near a lensing caustic cusp. The source itself appears to be mostly compact with an extended feature to the south east, which is readily visible in the observed lensed image.
There is no redshift measurement for this source and so we have opted to use the range of redshifts ( ∼ 2 − 3.7) present in our sample to calculate redshift dependent quantities. There are also no PACS flux density measurements for this source, and so we rely on the SPIRE and ALMA measurements to constrain the SED. The peak of the SED appears to lie within the SPIRE wavelengths, and fitting the optically thick SED for the range of redshifts considered gives a temperature range of 37 − −59 K. Without the PACS measurements and with the extra free parameters of the optically thin SED, it is not possible to meaningfully infer the presence of a warmer dust component.

DISCUSSION
Combining our sample with that of Dye et al. (2018), who carried out a similar analysis on a set of 6 strongly lensed sub-mm galaxies observed as part of the same ALMA programme, we can start to make more significant comparisons between our results and those found by other surveys. The HERschel ULIRG Survey (HERUS; Clements et al. 2017) is a sample of 43 local ULRIGS, selected at 60 m by the Infrared Astronomical Satellite (IRAS). The ALMA LESS survey is a set of cycle 0 and cycle 1 ALMA observations of the sub-mm sources detected in the LABOCA ECDFS Submm Survey (LESS). Comparing the relationship between the far-IR luminosity and the effective dust temperature (see Fig. 5) from the optically thick model between our sample, the ULIRG population from HERUS and the sub-mm sources in ALESS, we can see that our sources tend to possess both more extreme luminosities and higher Figure 5. Temperature distribution as a function of far-IR luminosity. The red crosses show the luminosities taken from this work (excluding HeLMS J235331.9+031718) and Dye et al. (2018). The spectroscopic redshift range of this combined dataset is approximately ∼ 1 − 4. The blue pluses show the data taken from the HERUS Survey (Clements et al. 2017) with a spectroscopic redshift range of approximately ∼ 0.02 − 0.26. The grey dots come from the ALESS survey (Swinbank et al. 2013) with a photometric redshift range of ∼ 0.5 − 6.5. dust temperatures. The median luminosity of the of the ALESS sample is L IR = (3.0 ± 0.3) × 10 12 L , compared with the median luminosity of the HERUS sample of L IR = (1.7 ± 0.3) × 10 12 L , and finally, that of our sample; L IR = (5.1 ± 2.5) × 10 12 L . Given our small sample size and large scatter on this measurement, overall, our luminosities are consistent with both the HERUS and ALESS samples. The median effective temperature of the ALESS sample is T = 31 ± 1K, which is consistent at the 2-level with that of the HERUS sample T = 38 ± 3K. The median dust temperature of our sample is significantly higher at = 59 ± 3K. The high dust temperatures can be interpreted as a product of the extreme SFRs present in our sample, though it is important to bear in mind that the significance of these results are likely explained by the combination selection effects due to our observing of the brightest galaxies detected by Herschel and the sample of ALESS galaxies being selected at 870 . Fig 4 shows the degree to which our sources exceed this ratio by plotting the SFR determined using the scaling relation of the far-IR luminosity given by Kennicutt & Evans (2012) against the dust mass. Given that our SFRs are derived from scaled far-IR luminosity, we can conclude that our sources have higher than expected luminosities for the amount of gas available for star formation. An obvious interpretation of this fact would be that a large fraction of the luminosity is due to an active galactic nucleus, but without additional imaging we cannot confirm this. It is also important to note that this excess could be at least in part explained by selection bias, as we have chosen the brightest Herschel sources to follow-up on.
Converting the rest frame 850 m flux density of our sources to H 2 gas mass (see Table 4), using the scaling relation given by Hughes et al. (2017), we find that our sources all lie on or above the mean relationship between SFR and H 2 gas mass as determined by Scoville et al. (2016). An interpretation of this is that our sources possess a higher star formation efficiency (SFE), provided that dust is an accurate tracer for molecular gas. Fitting a line parallel to the SMG/ULIRG relation shown in Fig. 4, treating the range of dust masses as the 1-error, we find an increase in SFE by a factor of 6 when compared with the value implied by the SMG/ULIRG relation from Rowlands et al. (2014). This factor increases to 50 when compared to the < 0.5 H-ATLAS galaxies.
An unknown fraction of our sources could have their SFR over estimated due to significant contamination in the IR by strongly obscured AGN. In order to quantify this effect, additional observational evidence would need to be considered, e.g estimating the stellar mass from broadband SED fitting, which would require careful lens light subtraction or additional follow-up observations in the X-ray hard band to reveal the obscured nucleus.

CONCLUSIONS
We have modelled seven ALMA observations of strongly lensed submillimetre galaxies. Four of these systems are galaxy-galaxy scale strong lenses, which are well described by a single power-law mass profile, whilst the remaining three are group-scale lenses, which have been successfully fitted with two power-law mass profiles. Where we found it improved the fit, an external shear term was also included in the lens model. In this work we have opted to model the visibility data directly rather than to work with CLEANed image plane data. Whilst the uv-plane method is more computationally expensive, it does not suffer from the image pixel covariances introduced due to incomplete sampling of the uv-plane. However, Dye et al. (2018) showed that with sufficient sampling of the uv space, both the image-plane and direct visibility modelling approaches produce very similar results. We have fitted smooth power-law density profiles (and in some cases two) for each of the lensing systems, and found that most of the lenses are close to isothermal. This result is expected in massive ETGs due to the combination of an inner Sérsic profile representing the baryonic component and an outer NFW profile representing the dark matter component (Lapi et al. 2012). In some instances, there are significant deviations from an isothermal power-law slope, which may be due to degeneracies between parameters in our model or reflect the true nature of the lens, but we leave a more thorough explanation of this to further work.
By obtaining the total magnification factors from our models, we have demagnified the available sub-millimetre source photometry. Fitting rest frame SEDs to this data allowed us to determine the dust temperature, dust mass, intrinsic luminosity and star formation rates of our lensed sources. To estimate a range of possible dust masses for these sources, we fitted the photometry with both single temperature optically thick and dual temperature optically thin SEDs. Using the midpoint of this range to calculate the SFR to dust mass ratio, we find that all seven of our sources lie above the mean ratio for the SMG/ULIRG population as described in Rowlands et al. (2014).
Galaxy morphology is strongly correlated with star formation history (Larson et al. 1980;Strateva et al. 2001;Lee et al. 2013). At ∼ 3 massive galaxies are mostly star forming disks, with the SFR peaking at 1 < < 2 (Madau & Dickinson 2014), and then dropping with the rapid growth of the fraction of massive quiescent galaxies. This period of galaxy evolution is extremely dramatic, and many different mechanisms have been proposed to explain the build up of of ETGs we see today. Galaxy mergers are one such mechanism as they are effective at disturbing the morphology and building a central bulge in a galaxy (Hopkins et al. 2006;Snyder et al. 2011). They have also been shown to trigger starbursts and AGN, which can lead to strong supernovae and/or AGN winds contributing to the quenching of a galaxy (Bekki et al. 2005). Our combined sample of strongly lensed sub-mm galaxies contains a majority of sources that display a disturbed morphology. 9 of the 13 galaxies in our sample are visually classified as being disturbed, with at least two of them having evidence of being mergers (H-ATLAS J142935.3-002836, H-ATLAS J083051.0+01322). Other observations of sub-mm galaxies have found similarly high fractions of disturbed morphologies, such as Chapman et al. (2004), who used high-resolution optical and radio imaging of 12 sub-mm galaxies to study their spatially extended star formation activity. It has been suggested that high density molecular gas is more commonly found in galaxy mergers than quiescent systems and that this can be used to predict the star formation mode of a galaxy (Papadopoulos & Geach 2012). We are not able to conclusively say what fraction of our sample's disturbed morphologies are a result of mergers, but the source with the most extreme ratio of SFR and gas mass (H-ATLAS J083051.0+01322) does display a significantly disturbed morphology and is identified as being a merger by Yang et al. (2019).