# SHARP – V. Modelling gravitationally lensed radio arcs imaged with global VLBI observations

C. Spingola<sup>1,★</sup> J. P. McKean,<sup>1,2</sup> M. W. Auger,<sup>3</sup> C. D. Fassnacht,<sup>4</sup> L. V. E. Koopmans,<sup>1</sup> D. J. Lagattuta<sup>1,5</sup> and S. Vegetti<sup>6</sup>

Accepted 2018 May 16. Received 2018 April 24; in original form 2017 December 14

## ABSTRACT

We present milliarcsecond (mas) angular resolution observations of the gravitationally lensed radio source MG J0751+2716 (at $z = 3.2$) obtained with global very long baseline interferometry (VLBI) at 1.65 GHz. The background object is highly resolved in the tangential and radial directions, showing evidence of both compact and extended structure across several gravitational arcs that are 200–600 mas in size. By identifying compact sub-components in the multiple images, we constrain the mass distribution of the foreground $z = 0.35$ gravitational lens using analytic models for the main deflector [power-law elliptical mass model; $\rho(r) \propto r^{-\gamma}$, where $\gamma = 2$ corresponds to isothermal] and for the members of the galaxy group. Moreover, our mass models with and without the group find an inner mass-density slope steeper than isothermal for the main lensing galaxy, with $\gamma_1 = 2.08 \pm 0.02$ and $\gamma_2 = 2.16 \pm 0.02$ at the $4.2\sigma$ level and $6.8\sigma$ level, respectively, at the Einstein radius ($b_1 = 0.4025 \pm 0.0008$ and $b_2 = 0.307 \pm 0.002$ arcsec, respectively). We find randomly distributed image position residuals of about 3 mas, which are much larger than the measurement errors ($40 \mu\text{as}$ on average). This suggests that at the mas level, the assumption of a smooth mass distribution fails, requiring additional structure in the model. However, given the environment of the lensing galaxy, it is not clear whether this extra mass is in the form of sub-haloes within the lens or along the line of sight, or from a more complex halo for the galaxy group.

**Key words:** gravitational lensing: strong-techniques: interferometric-galaxies: active-radio continuum: galaxies.

## 1 INTRODUCTION

In observational cosmology, gravitational lensing is the only method that allows one to directly probe the projected matter (including dark matter) density distribution of galaxies over an extended range of scales, independent of its dynamical state (see Treu 2010; for a recent review). Ever since the first discovery of multiple imaging (Walsh, Carswell & Weymann 1979), gravitational lensing has been used to investigate a broad range of astrophysical questions, from the structure of the large-scale matter distribution to the physical properties of the individual lensing galaxies (Koopmans et al. 2009; Giocoli et al. 2013; Sonnenfeld et al. 2015). In particular, gravitational lensing is a powerful technique to test models for global

halo profiles that are predicted from hierarchical galaxy formation simulations (e.g. Navarro, Frenk & White 1997). However, the constraints provided by observations of the multiple images alone are often not sufficient to determine a univocal lens mass model. It is for this reason that parametric models, motivated by the observed general properties of typical galaxies in the local Universe, are generally used to overcome this obstacle (Keeton 2001a; Wucknitz & Refsdal 2001).

The simplest macro model that can reproduce the relative positions and flux ratios of the multiple images is the singular isothermal ellipsoid (SIE) plus an external shear to account for neighbouring galaxies, which is often sufficient to describe the projected mass of elliptical lensing galaxies. For example, the majority of the lenses in the optical CfA-Arizona-Space Telescope Lens Survey (CASTLES) can be modelled by simple ellipsoidal mass distributions with an

<footnote n="1"><sup>1</sup>Kapteyn Astronomical Institute, University of Groningen, Postbus 800, NL-9700 AV Groningen, the Netherlands</footnote>

<footnote n="2"><sup>2</sup>ASTRON, Netherlands Institute for Radio Astronomy, Oude Hoogeveensedijk 4, NL-7991 PD Dwingeloo, the Netherlands</footnote>

<footnote n="3"><sup>3</sup>Institute of Astronomy, University of Cambridge, Madingley Rd, Cambridge CB3 0HA, United Kingdom</footnote>

<footnote n="4"><sup>4</sup>Department of Physics, University of California, Davis, CA 95616, USA</footnote>

<footnote n="5"><sup>5</sup>Univ Lyon, Univ Lyon1, Ens de Lyon, CNRS, Centre de Recherche Astrophysique de Lyon UMR5574, F-69230 Saint-Genis-Laval, France</footnote>

<footnote n="6"><sup>6</sup>Max Planck Institute for Astrophysics, Karl-Schwarzschild-Strasse 1, D-85740 Garching, Germany</footnote>

<footnote>* E-mail: spingola@astro.rug.nl</footnote>
external shear, and the number of gravitational lenses that require a deviation from this model are few (Muñoz et al. 1998; Falco et al. 1999). Among all of the known galaxy-scale gravitational lenses, those with extended images of the background lensed source provide the most constraints to test the mass model, and, therefore, they can be used to investigate whether an elliptical power-law density model is a more accurate description of the matter content of the lensing galaxy (e.g. Koopmans et al. 2009). Also, from comparing the constraints derived for data at progressively higher angular resolution, Lagattuta et al. (2012) find that the parameter space is highly constrained with better quality imaging.

In this respect, interferometric observations at radio wavelengths can currently provide the highest angular resolution imaging available, and surveys of radio-loud gravitational lens systems at high angular resolution have been carried out extensively (Burke 1990; Hewitt 1992; Patnaik 1993; Winn, Hewitt & Schechter 2001; Browne et al. 2003; Myers et al. 2003). However, most of the lensed sources discovered in this way have compact structure and the unresolved multiple images place only a few constraints on the lens mass model. This is because many of these objects were discovered through a systematic search of flat-spectrum radio sources, which are typically compact when observed with the very large array (VLA at 8.46 GHz; 170 mas beam size) and the Multi-Element Radio Linked Interferometry Network (MERLIN at 5 GHz; 50 mas beam size). In a few cases, the radio structure of the background source is resolved into Einstein rings and extended gravitational arcs (e.g. Hewitt et al. 1988; Langston et al. 1989; Biggs et al. 2001), as for example, the lenses MG J0414+0534 (MacLeod et al. 2013), MG B2016+112 (Koopmans et al. 2002; More et al. 2009), JVAS B1933+503 (Suyu et al. 2012), and CLASS B1555+375 (Hsueh et al. 2016). In these cases, the extra constraints provided by the extended arcs revealed that the macro model could not be explained by a simple SIE, but additional mass structures are necessary to reproduce the relative position and fluxes of the observed images.

Gravitational lensing is a powerful method to directly detect and quantify any deviation of a smooth mass model for the primary lens and it can be used with this aim in two different ways. One method uses the flux ratios of the multiple images of compact background sources to find evidence for peculiar magnifications, which can be due to a perturbation from small-scale structures in the lensing galaxy (Mao & Schneider 1998). For example, the study of seven radio-loud flux-ratio anomalous lenses by Dalal & Kochanek (2002) finds that the mass fraction of sub-structure required to reproduce the flux ratios within their sample was $f_{sub} = 2^{+5.0}_{-1.6}$ per cent (90 per cent confidence levels). The other method to detect perturbations to a smooth mass model consists of observing astrometric anomalies of the multiple images, namely observing lensed images in different positions from what is expected from a smooth mass distribution. Given the high angular resolution from very long baseline interferometry (VLBI) data at cm-wavelengths (2–10 mas beam size), deviations of a smooth macro model can be found via astrometric perturbations of the multiple images; the additional mass structure can perturb the deflection angle and, therefore, the relative positions of the images can be shifted. For example, simulations predict that a dark matter sub-halo of a mass $10^8 M_{\odot}$ in a Milky Way-sized galaxy can produce an astrometric perturbation of ~10 mas in the lensed images (Metcalf & Madau 2001). Even if the sub-halo has a lower mass (>10⁵ $M_{\odot}$), it is still possible to observe local independent distortions in the VLBI images of lensed radio jets (Metcalf 2002). These distortions, which appear as bends in the jets, can be detected by measuring the local curvature of extended multiple images and noting differences between points that correspond to

the same source position. Therefore, observations with mas (and sub-mas) angular resolution are key to testing the smoothness of the macro models.

With these aims, the strong lensing at high angular resolution program (SHARP; Fassnacht et al., in preparation) has carried out high angular resolution observations at optical, near-infrared (NIR), and radio wavelengths of known gravitational lenses with extended source structure. The main goal is to detect and measure possible perturbations to the macro models associated with the main lensing galaxy halo to constrain models for galaxy formation, dark matter, and cosmology. Thus far, SHARP has focussed on using high angular resolution observations with the adaptive optics system on the W. M. Keck 10-m Telescope and the *Hubble Space Telescope (HST)* to test global mass models (Lagattuta, Auger & Fassnacht 2010; Lagattuta et al. 2012; Chen et al. 2016), infer the properties of luminous and dark dwarf galaxies (sub-structures; McKean et al. 2007; Vegetti et al. 2012), and investigate the mass perturbations due to galactic-scale discs (Hsueh et al. 2016, 2017). Here, we extend SHARP to higher angular resolution (by a factor of 30–60) and longer wavelengths by presenting imaging of the most extended gravitational arc on mas-scales known with VLBI at cm-wavelengths.

The first target is the strongly lensed radio-loud quasar MG J0751+2716, which is one of the most promising targets for studying the smoothness of the lensing mass distribution on mas-scales. The gravitational lens was discovered as part of the VLA follow-up of sources found from the MIT–Green Bank survey, and shows large gravitational arcs at high surface brightness when observed with MERLIN at 5 GHz (50 mas beam size; Lehar et al. 1997). Optical imaging shows that the lensing galaxy is a satellite of a bright cluster galaxy (BCG) and is part of a small group of galaxies. The redshift of this group of galaxies was found to be $z_{group} = 0.35$ (Tonry 1998; Momcheva et al. 2006) and the redshift of the background source is $z = 3.200$ (Tonry 1998; Alloin et al. 2007). The macro models proposed to date require a significant external shear, likely due to the group of galaxies in the field (Lehar et al. 1997; Momcheva et al. 2006; Alloin et al. 2007). The large extent of the arcs coupled with the bright emission of this source (~350 mJy at 1.7 GHz) make MG J0751+2716 an excellent lens system to study the level of deviations from a smooth macro model.

In this paper, we present new global VLBI observations of MG J0751+2716 at 1.65 GHz with the main aim of investigating the global mass model and determining a precise radial density profile of the mass distribution (accounting for the galaxy group). Our paper is organized as follows. In Section 2, we describe the observations and imaging results, and in Section 3, we illustrate the improved macro models that can be determined with mas-scale angular resolution imaging. Our discussion and summary of results are presented in Sections 4 and 5, respectively. Throughout this paper, we assume $H_0 = 67.8 \text{ km s}^{-1} \text{ Mpc}^{-1}$, $\Omega_M = 0.31$, and $\Omega_\Lambda = 0.69$ in a flat Universe (Planck Collaboration XIII 2016).

## 2 OBSERVATIONS

In this section, we present the global VLBI observations of MG J0751+2716 and optical imaging of the field galaxies that we use as additional components to the lensing mass model.

### 2.1 Very long baseline interferometry data

MG J0751+2716 was observed at a central frequency of 1.65 GHz with the global VLBI array on 2012 October 21 for a total time of
18.5 h (project GM070; PI: McKean). The observation comprised 24 antennas from the European VLBI Network (EVN) and the Very Long Baseline Array (VLBA), and included the large (>50 m) Lovell, Effelsberg, Robledo, and Green Bank telescopes. The scans on the target were about 3 min in duration, which were interleaved by scans of about 2 min on the phase-reference source J0746+273. Several observations of the bright calibrator sources 4C39.25 and DA193 were taken throughout the run for fringe finding during correlation and for the bandpass calibration at the data reduction stage. The data were recorded at 512 Mbits s⁻¹ and correlated at the Joint Institute for VLBI in Europe (JIVE) to produce eight spectral windows (IFs) with 8 MHz bandwidth and 32 channels each, through both circular polarizations (RR, LL). A visibility averaging time of 2 s was used. This time and channel resolution limited the effective field of view of the observations to about 16 and 10 arcsec, respectively, from the phase centre, which easily encompassed all of the expected structure of the target.

The data set was initially edited, calibrated, and reduced using the EVN pipeline and the Astronomical Image Processing Software (AIPS) to produce/apply standard calibration tables. However, during the fringe-fitting process, three antennas (Shanghai, Urumqi, and Svetloe) were lost because they had baselines with a signal-to-noise ratio of <5. After the initial calibration was completed, new models for the calibrators, and target were obtained, which were then used to rerun the fringe-fitting process; this additional step improved the corrections for the residual fringe rates and delays on all antennas, with the exception of Svetloe, which could not be recovered. The phase-referenced data set for MG J0751+2716 was then imaged and self-calibrated. A solution interval of 120–30s was used to perform several iterations of phase-only self-calibration. Finally, amplitude self-calibration was applied using at first a long solution interval (lasting the whole observation for each antenna) that was gradually reduced to 30 min to remove any residual uncertainties from the antenna gains. We note that the self-calibration process leads to a global shift in the absolute position of the lensed images by about 1 mas, that is, a fraction of the synthesized beam size. However, this does not affect our gravitational lens modelling, because we use the relative positions of the lensed images, which are not changed by the self-calibration process.

As the system is quite complex and extended, the imaging was performed using multiscale cleaning, which is more efficient at modelling extended structures (Cornwell 2008; Rich et al. 2008), within the WSCLEAN algorithm (Offringa et al. 2014). Our final image of MG J0751+2716 is presented in Fig. 1, which was obtained by using uniform weights for the visibilities; the off-source rms is 41 µJy beam⁻¹ and the peak surface brightness is 2.9 mJy beam⁻¹. The restored beam is 5.5 × 1.8 mas² at a position angle −9.8 deg east of north.

In Fig. 1, we see that the extended arcs that were previously detected with MERLIN by Lehar et al. (1997) are now resolved into several sub-components that are connected via diffuse jet structure. Components A and C are resolved into four sub-components, while components B and D are resolved into six sub-components. The pair of merging images (A4 and C4) are highly distorted in the tangential direction with a low-flux density emission. Components A(1→3) have a similar morphology of components D(1→3), distorted both in the radial and tangential direction. In addition, the counter image of the doubly imaged part (components B6 and D6) of the radio source is detected for the first time. We also detect a new source components (B7) at the 4σ level in the doubly imaged region. Never before have such extended gravitational arcs been observed at such a high angular resolution. This demonstrates the excellent UV

coverage and surface brightness sensitivity provided by the global VLBI array (Fig. 2), which is fundamental for a detailed study of the structure of extended arcs on mas-scales from objects like MG J0751+2716. For example, the global VLBI array sensitivity is 2.5 times better than an EVN only observation and 10 times better than a VLBA only observation. The total flux density of MG J0751+2716 was determined in the image plane by placing an aperture over the area that contains the arcs and the double component, and was found to be $S_{1.65GHz} = 350 \pm 35$ mJy.

We conservatively assume an uncertainty on the absolute flux density scale of ~10 per cent, based on the gain and system temperature variations during the observation.

## 2.2 Hubble Space Telescope data

As MG J0751+2716 is known to be gravitationally lensed by a group of galaxies, we use high resolution optical imaging from the HST to provide the relative positions, ellipticities, and position angles of the group galaxies as an input to our mass modelling. The archival optical observations of MG J0751+2716 were obtained as part of the CASTLES program (GO-7495; PI: Falco) using the Wide Field Planetary Camera 2 (WFPC2) through the F814W filter. These observations were processed in ASTRODRIZZLE within the IRAF package using standard procedures. The final drizzled image has a pixel scale of 0.045 arcsec pixel⁻¹ and is shown in Fig. 3. In order to estimate the relative position and magnitude of the galaxy group members in this HST image, we use the software SExTRACTOR (Bertin & Arnouts 1996).

## 3 LENS MODELLING

The high angular resolution of the data and the wealth of extended structure that has been detected from MG J0751+2716 can provide many constraints to the gravitational lens mass model. In principle, such an analysis should include all of the structure that is observed as a constraint. However, this requires producing a model for the lens and background source that fits the observed data, which in this case are the visibilities. Since producing such a model is computationally expensive, we instead start with a simpler case: we generate a model by matching the conjugate positions of compact sub-components seen in each multiply imaged arc. The more sophisticated grid-based modelling of the data will be presented in a future paper.

### 3.1 Parametric lens modelling

We use the publicly available code GRAVLENS (Keeton 2001a,b) to model the compact radio components. As primary constraints to the starting model, we use the relative positions of four source components that are quadruply imaged (from the four extended arcs) and two components that are doubly imaged. The relative positions of these components are listed in Table 1 and are shown in Figs 4 and 5. We also give the peak surface brightness of each component in Table 1, but these are not used as model constraints for all of the components except of the two doubly imaged components; since they are both compact and fairly isolated, we use their relative flux ratios as additional model constraints.

We first identify groups of lensed images that correspond to the same background source component. Next, we fit each of these image components with a single elliptical Gaussian model in the image plane using the task JMFIT within AIPS to determine their position. In this way, we take into account the extended nature of each source component and obtain a weighted position for the individual lensed
<figure><img src="image_1.png" /><figcaption>**Figure 1.** Global VLBI imaging of MG J0751+2716 at 1.65 GHz obtained by using uniform weights and multiscale cleaning in WSCLEAN. The off-source rms is 41μJy beam⁻¹ and the peak surface brightness is 2.9 mJy beam⁻¹. The restored beam is 5.5 × 1.8 mas² at a position angle −9.8 deg, and is shown within the white box in the bottom left-hand corner.</figcaption></figure>

images. We apply a different approach for the groups of image components 4 and 7, which are either highly distorted or not detected in more than one of the lensed images. We instead use these image components as a test for our best model, as opposed to using them as constraints. For source component 4, we could confidently identify the more compact image components B4 and D4, while A4 and C4 are highly distorted by the lens (e.g. see Fig. 1). Therefore, we use the model-predicted position for A4 and C4 to fit a Gaussian in the image plane and obtain the weighted average position for these two images, which we show in Figs 4 and 5. We follow the same method for finding the predicted position and flux density of image component D7. The relative model-predicted flux density of image component D7 is about 0.16 of image component B7, which is detected at the 4σ level. The flux density of image component B7 is 223 ± 55 μJy. Therefore, the model-predicted flux density for D7 is ~36 μJy, which would correspond to a non-detection in our image, and is consistent with the data.

The positional uncertainty on each component is calculated in the standard way by using the major and minor axes of the elliptical Gaussian determined with JMFIT, and their signal-to-noise ratio (calculated using the peak surface brightness given in Table 1 and the

rms of the uniform weighted image). We assume larger positional uncertainties for the relative declination of the components A4 and C4, because they show a significant distortion in the declination direction. With this assumption, they do not significantly effect the final χ² of the lens model.

Moreover, since we use the relative position of the images as constraints to the lens modelling, we can neglect the systematic errors in the estimate of the absolute positional accuracy. The main lensing galaxy (G3 following the nomenclature of Lehar et al. 1997; see also Fig. 3) is known to be an elliptical galaxy, as shown from the surface brightness profile at near-infrared and optical wavelengths (e.g. Lehar et al. 1997; Spingola et al. in preparation) and from the optical spectrum (Tonry 1998). Therefore, we model the lensing potential as a power-law ellipsoid density profile, which has been a remarkably good fit to early-type lensing galaxies that were observed with the HST and Keck adaptive-optics imaging (Lagattuta et al. 2012; Vegetti et al. 2014; Oldham et al. 2017). However, in the case of MG J0751+2716, there is known to be a significant perturbation to the lensing potential, which could be due to the group of galaxies associated with the main lensing galaxy (Lehar et al. 1997; Momcheva et al. 2006; Alloin et al. 2007). Therefore, we
<figure><img src="image_2.png" /><figcaption>**Figure 2.** The UV coverage of the global VLBI observations of MG J0751+2716 at 1.65 GHz.</figcaption></figure>

take this into account by considering two models; Model 1 includes the main lensing galaxy and an additional external shear component, while Model 2 includes the lensing galaxy, external shear, and five additional haloes representing all spectroscopically confirmed group-member galaxies (see Fig. 3 and Section 3.2).

As there is a known degeneracy between the ellipticity and the external shear, we perform the optimization by using the following method. We first find the optimal position angles and lens strength for G3 by using the values obtained by Lehar et al. (1997) as an initial guess. We then keep the ellipticity and shear fixed, and optimize for the position angle of each. For the next step, we optimize for all of the parameters by keeping the ellipticity and shear fixed. Finally, we leave all parameters free and optimize for everything. From the positions of the four quadruply imaged components and the two doubly imaged components measured from our global VLBI observations, we obtain 36 positional constraints that are also probing the lensing potential over a large region (e.g. see Fig. 4).

For Model 1, we have 12 variables to describe the six source positions, and the lens mass model has 8 variables [mass scale (b); position (x<sub>L</sub>, y<sub>L</sub>); ellipticity (e), position angle (θ); power-law slope (γ), external shear (Γ), and its position angle (Γ<sub>θ</sub>)]. For Model 2, we have the same number of model parameters, except that we fix the position of the galaxies relative to G3 based on their centroid positions from HST imaging (G3 position is free; x<sub>L</sub>, y<sub>L</sub>), and we fix the mass of the galaxies relative to G3 (which is free; b) using their relative optical magnitudes and scaling relations according to their Hubble type (Faber–Jackson or Tully–Fisher, e.g. McKean et al. 2005; More et al. 2008). We infer the best-fitting parameters and the uncertainty on the values (68 per cent confidence level) from a Markov-Chain Monte Carlo (MCMC) sampler implemented in GRAVLENS.

## 3.2 Results

The results for Models 1 and 2 are given in Table 2, and the observed and model-predicted image positions, with the critical and caustic curves, are shown in Figs 4 and 5, respectively. The residual image

positions between the data and models are shown in Fig. 6. The probability density distribution for each parameter is shown in Figs 7 and 8 for Models 1 and 2, respectively. We now briefly describe the results for the two models considered here.

### 3.2.1 Model 1 – A single lensing galaxy and external shear

We first test whether the image configuration can be explained with a simple ellipsoidal power-law density profile, with an external shear. This model is similar to previous single-lens models that assumed an isothermal density profile (Lehar et al. 1997), but has much tighter constraints on all of the parameters given the increased precision provided by the mas-scale resolution of our global VLBI observations.

The Einstein radius of G3 is found to be $b = 0.4025^{+0.0007}_{-0.0008}$ arcsec, with an ellipticity of $e = 0.159 \pm 0.001$ at a position angle of 36 $\pm$ 3 deg (east of north). The ellipticity and position angle of the surface brightness profile as measured in the HST F814W imaging using the software GALFIT (Peng et al. 2010) are 0.35 $\pm$ 0.04 and 16 $\pm$ 2 deg, respectively. This misalignment of $\sim$20 deg between the mass distribution and the light profile is not unexpected due to the significant external perturbation from the other galaxies in the group, which can affect the shape of the lensing potential (Keeton, Christlein & Zabludoff 2000; Kochanek 2002). This level of misalignment has also been found in other lensing groups (e.g. CLASS B2108+213; More et al. 2009; McKean et al. 2010) and in lensing systems with a substantial external shear (e.g. Gavazzi et al. 2012). Indeed, there is a significant amount of external shear needed by this model ($\Gamma = 0.084 \pm 0.005$), and its position angle ($\Gamma_{\theta} = 79 \pm 2$ deg) suggests that the BCG of the group, to the west of the main lensing galaxy, is the principal cause of this external perturbation (see Fig. 3). Such a high external shear is typical for galaxies lying in a group or cluster of galaxies (e.g. Keeton, Kochanek & Seljak 1997; Oguri, Keeton & Dalal 2005; Auger et al. 2007). Finally, the best-fitting model suggests a power-law density slope for the main lensing galaxy that is steeper than isothermal ($\gamma = 2.08 \pm 0.02$) at the $4.2\sigma$ level.

This model does not provide a satisfactory fit to the observed images, as shown in Table 1 and Fig. 6. The most difficult images to fit are A4 and C4, because they are highly distorted and, therefore, are not properly represented by a single elliptical Gaussian component; any small change in the position of the source-component four will have a significant change in the position of image-components A4 and C4 due to its position relative to the caustic. Also, these images have larger uncertainties on their positions with respect to the other images (Table 1), therefore they are not adequately constrained during the optimization and they do not dominate the $\chi^2$ minimization. Moreover, we also find for the other images a significant mismatch between their observed and predicted positions, with respect to the astrometric uncertainty, with mean offsets of the order of $\sim$200$\sigma$.

### 3.2.2 Model 2 – A lensing group of galaxies and external shear

For our second model, we explore a more realistic lens mass model for MG J0751+2716 in which we take into account the individual members of the group of galaxies explicitly. The mass distribution of the entire group is parametrized using an ellipsoidal power-law density profile for G3 (the main lensing galaxy), with five SIEs to represent each member of the group of galaxies (G1, G2, G4, G5, and G6), plus an external shear component. Their position and mass scales relative to G3 are fixed, as also is their ellipticity and position
<figure><img src="image_3.png" /><figcaption>**Figure 3.** (Left) HST WFPC2 F814W image of the MG J0751+2716 group of galaxies, where G3 is the main lensing galaxy. The field of view of this image is 15 × 15 arcsec². The nomenclature for the galaxies follows Lehar et al. (1997) and Alloin et al. (2007) for G1–G6, which are all spectroscopically confirmed group members (Momcheva et al. 2006). (Right) The dimensionless convergence map for Model 2 (see Section 3), showing the combined contribution of the individual group galaxies to the mass model. GA and GB are not included in Model 2 since they do not have spectroscopic information. The contours show regions of iso-convergence for $κ = 0.05, 0.1, 0.2$, and $0.4$.</figcaption></figure>

<table><caption>**Table 1.** The measured positions of the various components used for lens modelling, as indicated in Figs 4 and 5. The observed relative RA and Dec are determined by performing elliptical Gaussian fits using JMFIT within AIPS, and are measured with respect to component A₁ (phased referenced absolute position 07h51m41.487s, +27°16'31.621'). The position of component D₇ is the model-predicted position. The predicted positions from the lens models tested here are also relative to the observed position of component A₁. The peak surface brightness (mJy beam⁻¹) is given as a reference, and is not used as a constraint for the lens models. The reported error on the peak surface brightness is the nominal error of the Gaussian fit.</caption><thead><tr><th rowspan="2">ID</th><th colspan="2">Observed</th><th colspan="2">Model 1</th><th colspan="2">Model 2</th><th rowspan="2">I<br>(mJy beam<sup>-1</sup>)</th></tr><tr><th>α (mas)</th><th>δ (mas)</th><th>α (mas)</th><th>δ (mas)</th><th>α (mas)</th><th>δ (mas)</th></tr></thead><tbody><tr><td>A<sub>1</sub></td><td>0.00 ± 0.01</td><td>0.00 ± 0.02</td><td>-3.27</td><td>+2.86</td><td>-2.64</td><td>+2.56</td><td>4.91 ± 0.63</td></tr><tr><td>A<sub>2</sub></td><td>+6.87 ± 0.02</td><td>+35.60 ± 0.20</td><td>+2.05</td><td>+27.86</td><td>+1.28</td><td>+29.81</td><td>2.29 ± 0.67</td></tr><tr><td>A<sub>3</sub></td><td>+9.99 ± 0.02</td><td>+20.29 ± 0.10</td><td>+12.81</td><td>+21.78</td><td>+11.34</td><td>+24.63</td><td>2.27 ± 0.08</td></tr><tr><td>A<sub>4</sub></td><td>+47.50 ± 0.07</td><td>+21.50 ± 1.00</td><td>+52.89</td><td>+37.59</td><td>+64.14</td><td>+26.02</td><td>0.67 ± 0.02</td></tr><tr><td>B<sub>1</sub></td><td>-493.31 ± 0.01</td><td>-279.25 ± 0.03</td><td>-487.28</td><td>-277.09</td><td>-485.34</td><td>-279.18</td><td>2.78 ± 0.22</td></tr><tr><td>B<sub>2</sub></td><td>-414.61 ± 0.03</td><td>-298.87 ± 0.03</td><td>-410.23</td><td>-292.15</td><td>-413.28</td><td>-294.04</td><td>1.20 ± 0.02</td></tr><tr><td>B<sub>3</sub></td><td>-379.66 ± 0.04</td><td>-302.01 ± 0.20</td><td>-382.89</td><td>-298.69</td><td>-387.79</td><td>-300.46</td><td>1.20 ± 0.03</td></tr><tr><td>B<sub>4</sub></td><td>-261.65 ± 0.02</td><td>-304.10 ± 0.04</td><td>-262.34</td><td>-306.69</td><td>-262.03</td><td>-308.44</td><td>1.96 ± 0.18</td></tr><tr><td>B<sub>6</sub></td><td>-58.30 ± 0.30</td><td>-318.10 ± 0.40</td><td>-58.35</td><td>-320.85</td><td>-60.34</td><td>-321.27</td><td>0.11 ± 0.01</td></tr><tr><td>B<sub>7</sub></td><td>-164.37 ± 0.01</td><td>-300.44 ± 0.05</td><td>-163.19</td><td>-299.46</td><td>-162.94</td><td>-301.33</td><td>1.24 ± 0.06</td></tr><tr><td>C<sub>1</sub></td><td>-81.45 ± 0.01</td><td>+501.27 ± 0.02</td><td>-83.00</td><td>+501.61</td><td>-86.97</td><td>+499.99</td><td>3.25 ± 0.23</td></tr><tr><td>C<sub>2</sub></td><td>-29.85 ± 0.02</td><td>+451.82 ± 0.13</td><td>-25.68</td><td>+456.31</td><td>-26.95</td><td>+451.90</td><td>1.53 ± 0.05</td></tr><tr><td>C<sub>3</sub></td><td>-5.82 ± 0.03</td><td>+424.10 ± 0.01</td><td>-2.31</td><td>+426.19</td><td>-2.98</td><td>+421.22</td><td>1.44 ± 0.05</td></tr><tr><td>C<sub>4</sub></td><td>+59.50 ± 0.01</td><td>+278.50 ± 0.80</td><td>+84.93</td><td>+242.38</td><td>+85.21</td><td>+215.30</td><td>0.68 ± 0.05</td></tr><tr><td>D<sub>1</sub></td><td>-643.05 ± 0.02</td><td>+346.88 ± 0.03</td><td>-649.74</td><td>+346.88</td><td>-649.57</td><td>+350.76</td><td>2.48 ± 0.38</td></tr><tr><td>D<sub>2</sub></td><td>-639.02 ± 0.02</td><td>+328.94 ± 0.05</td><td>-643.21</td><td>+328.32</td><td>-644.16</td><td>+326.69</td><td>1.97 ± 0.12</td></tr><tr><td>D<sub>3</sub></td><td>-639.23 ± 0.02</td><td>+329.81 ± 0.05</td><td>-635.49</td><td>+327.41</td><td>-636.86</td><td>+326.23</td><td>2.04 ± 0.09</td></tr><tr><td>D<sub>4</sub></td><td>-607.80 ± 0.02</td><td>+321.78 ± 0.02</td><td>-606.77</td><td>+316.33</td><td>-607.30</td><td>+315.49</td><td>1.53 ± 0.12</td></tr><tr><td>D<sub>6</sub></td><td>-503.82 ± 0.03</td><td>+290.59 ± 0.04</td><td>-510.54</td><td>+291.12</td><td>-509.35</td><td>+283.97</td><td>0.52 ± 0.02</td></tr><tr><td>D<sub>7</sub></td><td>-574.03</td><td>+302.87</td><td>-574.03</td><td>+302.87</td><td>-574.03</td><td>+302.87</td><td>&lt;0.01</td></tr></tbody></table>

angle, based on the HST imaging (Table 2). A convergence map for this model is shown in Fig. 3, while the marginalized posterior PDF for the lens model parameters is shown in Fig. 8.

The immediate difference between Models 1 and 2, is a lowering of the Einstein radius of G3 to $b = 0.307 \pm 0.002$ arcsec and a lowering of the external shear to $\Gamma = 0.034 \pm 0.003$, as expected since we are now accounting for the external convergence of the system due to the group of galaxies. The change of the shear position angle to $\Gamma_{\theta} = -61.8 \pm 3.5$ deg can be attributed to the added

complexity of the galaxy environment. The ellipticity of G3 is found to increase to $e = 0.1605 \pm 0.0002$, and the misalignment between the optical surface brightness profile and the gravitational lensing mass profile is still of the order of 20 deg. It is not clear whether this misalignment is due to some additional mass structure that is not included in the mass model (see Section 4 for discussion), or whether this is evidence for some interaction in the group environment that has affected the G3 dark matter halo to a larger extent than the stellar component. Finally, the power law of the ellipsoidal
<figure><img src="image_4.png" /><figcaption>**Figure 4.** Model 1 for MG J0751+2716. The model consists of an ellipsoid power-law mass distribution for the main lensing galaxy G3 (red cross), with an external shear. The observed positions are the open circles and the predicted positions are represented by the crosses, with all positions given relative to component A1. Each colour corresponds to a different background source component (filled circles). The lens critical curve is shown by the thick line; the dashed line shows the source plane caustics. The grey lines are the $3\sigma$ contours of the extended emission detected from our global VLBI imaging, for reference. The mass model parameters and their uncertainties are given in Table 2.</figcaption></figure>

density profile for G3 has become even more superisothermal, with $\gamma = 2.16 \pm 0.02$, at the $6.8\sigma$ level.

Similar to Model 1, we find that there are significant deviations between the observed and model-predicted positions of the image-components (1–7 mas, that correspond to ~200σ offsets on average; see Table 1 and Fig. 6) and the positional offsets of Model 2 are 1 per cent larger than those of Model 1. We have carried out some additional tests of Model 2 to investigate what other mass structures could account for these differences, although in general, we find that the data are not sufficient to constrain these additional mass components, so they are not formally part of Model 2.

We first attempted to constrain a possible common halo for the group. Such a model was also tested by Alloin et al. (2007), who used a truncated pseudo-isothermal profile to account for the common group halo. Here, we use the more realistic case of including an NFW halo in addition to the six individual galaxies that make up the system. However, the position of the common halo and the mass-scale were poorly constrained by the data. Also, because we do not want to impose a bound state among these galaxies, we did not include this common halo as part of Model 2. Further observations at, for example, X-ray wavelengths may well constrain the position of a common halo, and reveal whether it is in a relaxed or disturbed dynamical state (e.g. Fassnacht et al. 2008).

The position angle of the external shear changes to $-61.2 \pm 2$ deg, which may suggest that any additional mass component causing

this residual shear should be in the direction of G2 and G6 (see Fig. 3). As a test, we included GA and GB (not spectroscopically confirmed as group members) in the model and found that there is a negligible change in the offsets between the observed and model predicted positions of the image-components, which is not surprising given the small contribution that they make to the total convergence. Indeed, by including these two galaxies the shear strength and its position angle do not change with respect to the Model 2 values. Therefore, we can exclude them as the additional mass component that may be responsible for the residual external shear.

Finally, we note that although the structure of the background source is relatively unchanged when the group is included, the position of the background source does change to the extent that it no longer sits within the Einstein radius of G3 (Fig. 5). Therefore, without the additional convergence provided by the group of galaxies, MG J0751+2716 would not be strongly gravitationally lensed.

# 4 DISCUSSION

In this section, we discuss our lens models for MG J0751+2716 obtained from the VLBI imaging presented here.
<figure><img src="image_5.png" /><figcaption>**Figure 5.** Model 2 for MG J0751+2716. The model consists of an ellipsoid power-law mass distribution for the main lensing galaxy G3 (blue cross) and five SIEs for the group galaxies, with an external shear. The observed positions are the open circles and the predicted positions are represented by the crosses, with all positions given relative to component A1. Each colour corresponds to a different background source component (filled circles). The lens critical curve is shown by the thick line; the dashed line shows the source plane caustics. The grey lines are the $3\sigma$ contours of the extended emission detected from our global VLBI imaging, for reference. The mass model parameters and their uncertainties are given in Table 2.</figcaption></figure>

<table><caption>**Table 2.** The minimum-$\chi^2$ parameters of the two parametric lens models for MG J0751+2716 presented in this paper: *b* is the lens strength (arcsec), $x_L$ and $y_L$ are the positions in RA and Decl relative to A1 (arcsec), *e* is the ellipticity, $\theta$ is the position angle of the ellipticity (east of north, degrees), $\Gamma$ is the external shear strength, and $\Gamma_{\theta}$ is the external shear position angle (east of north, degrees). The density slope of the ellipsoidal power-law mass distribution is given by $\gamma$, where for an isothermal profile $\gamma = 2$. For Model 2, the positions and ellipticities are fixed to the optical parameters for the group galaxies and their $\gamma$ is fixed to the isothermal case. Also, for Model 2, the Einstein radius of the group galaxies relative to G3 is fixed based on their relative optical fluxes and using the Faber–Jackson or Tully–Fisher relation, depending on Hubble type: G1 and G3 are early-type galaxies, the others are late-type galaxies. For G3 we report the best set of parameters recovered via the minimization with GRAVLENS (Best) and the average values with relative 95 per cent limits assessed by the MCMC chains implemented in GRAVLENS (Mean).</caption><thead><tr><th rowspan="2">Par.</th><th colspan="4">Model 1</th><th colspan="6">Model 2</th></tr><tr><th>Mean</th><th>G3<br>σ<sub>95</sub><sup>mean</sup><br>per cent</th><th>Best</th><th>Mean</th><th>G3<br>σ<sub>95</sub><sup>mean</sup><br>per cent</th><th>Best</th><th>G1</th><th>G2</th><th>G4</th><th>G5</th><th>G6</th></tr></thead><tbody><tr><td>b</td><td>0.40249</td><td>+0.00074<br>-0.00081</td><td>0.39810</td><td>0.3073</td><td>+0.0021<br>-0.0022</td><td>0.3136</td><td>≡</td><td>≡</td><td>≡</td><td>≡</td><td>≡</td></tr><tr><td>x<sub>L</sub></td><td>-0.3530</td><td>+0.0011<br>-0.0011</td><td>-0.3561</td><td>-0.35482</td><td>+0.00018<br>-0.00018</td><td>-0.35448</td><td>0.720<sup>+0.098</sup><sub>-0.096</sub></td><td>0.250<sup>+0.110</sup><sub>-0.120</sub></td><td>0.320<sup>+0.170</sup><sub>-0.160</sub></td><td>0.061<sup>+0.041</sup><sub>-0.041</sub></td><td>0.265<sup>+0.052</sup><sub>-0.050</sub></td></tr><tr><td>y<sub>L</sub></td><td>0.1594</td><td>+0.0010<br>-0.0010</td><td>0.1624</td><td>0.16240</td><td>+0.00019<br>-0.00019</td><td>0.16273</td><td>≡ +5.447</td><td>≡ +2.120</td><td>≡ +6.530</td><td>≡ -4.807</td><td>≡ +1.066</td></tr><tr><td>e</td><td>0.159</td><td>+0.014<br>-0.013</td><td>0.2269</td><td>0.1605</td><td>+0.0080<br>-0.0085</td><td>0.1896</td><td>≡ 0.34</td><td>≡ 0.40</td><td>≡ 0.40</td><td>≡ 0.60</td><td>≡ 0.40</td></tr><tr><td>θ</td><td>35.7</td><td>+3.4<br>-3.3</td><td>49.0</td><td>38.9</td><td>+1.9<br>-2.0</td><td>44.9</td><td>≡ +1.0</td><td>≡ -70.0</td><td>≡ +40.0</td><td>≡ -82.0</td><td>≡ -50.0</td></tr><tr><td>Γ</td><td>0.0837</td><td>+0.0049<br>-0.0053</td><td>0.06605</td><td>0.0343</td><td>+0.0026<br>-0.0026</td><td>0.03109</td><td></td><td></td><td></td><td></td><td></td></tr><tr><td>Γ<sub>θ</sub></td><td>79.2</td><td>+1.6<br>-1.5</td><td>78.99</td><td>-61.8</td><td>+3.4<br>-3.7</td><td>-59.8</td><td></td><td></td><td></td><td></td><td></td></tr><tr><td>γ</td><td>2.079</td><td>+0.019<br>-0.019</td><td>2.008</td><td>2.157</td><td>+0.023<br>-0.023</td><td>2.078</td><td>≡ 2.0</td><td>≡ 2.0</td><td>≡ 2.0</td><td>≡ 2.0</td><td>≡ 2.0</td></tr><tr><td>χ<sup>2</sup><sub>red</sub></td><td colspan="4">1.9</td><td colspan="6">1.4</td></tr></tbody></table>
<figure><img src="image_6.png" /><figcaption>**Figure 6.** The offset between the observed and the model-predicted positions in units of sigma (upper) and mas (lower) for Models 1 (left) and 2 (right). Each colour and symbol represents a different group of lensed images as indicated in the legend on the top-right of each panel. The error bars are shown in black and the two black dashed lines indicate the no offset position.</figcaption></figure>

## 4.1 Precision lens modelling with global VLBI observations

MG J0751+2716 is one of the few quadruply imaged radio-loud gravitationally lensed quasars that show extended arcs on VLBI-scales. Our deep imaging detects the extended arcs at high significance, showing the complex surface brightness structure of the background source in unprecedented detail (Fig. 1). Never before have such extended (200–600 mas) gravitational arcs been detected at an angular resolution of a few mas. This detection allowed us to confidently identify lensed emission corresponding to the same source component, providing a very large number of constraints to the mass model that also sampled a large radial and tangential extent (Table 1). Moreover, these observations detect the faint counter-image B6 for the first time, providing a new additional constraint to the radial mass distribution of the lens. This component could be detected because of the excellent µJy beam⁻¹ sensitivity of the data, thanks to the large data recording rate (0.5 Gbit s⁻¹) and number of antennas used for this observation (Fig. 2). The advent

of even larger recording rates (at 2–4 Gbit s⁻¹) and global VLBI arrays that contain over 25 antennas will routinely provide µJy beam⁻¹ surface-brightness sensitivities and excellent UV coverage in a single synthesis observation.

Using the constraints provided by the global VLBI imaging of MG J0751+2716, we were able to infer the lens parameters with a high precision. By explicitly including the group of galaxies in the macro model, the uncertainties on the parameters are reduced significantly. For example, the lensing galaxy position is recovered with a precision of 0.6 and 0.1 per cent, respectively, for Models 1 and 2 (see Table 2), which also corresponds to a factor of 10 improvement in precision with respect to previous modelling that used MERLIN observations (50 mas FWHM beam size; rms 89 µJy beam⁻¹; Lehar et al. 1997). Moreover, our models constrain the Einstein radius with ~0.6 per cent precision and the ellipticity at the order of 0.5 per cent (Table 2). The slope of the mass density distribution $\gamma$ is found to be steeper than isothermal at the 4.2σ
<figure><img src="image_7.png" /><figcaption>**Figure 7.** The marginalized posterior probability distribution function (PDF) for the lens model parameters of Model 1. The contours enclosing the 99, 95, and 68 percentiles indicate the distribution between two parameters of the lens model. The PDF of each parameter is shown at the top of each column. The meaning of the parameters, the maximum-likelihood model values for each parameter, and their uncertainties, are presented in Table 2.</figcaption></figure>

level for Model 1 and at the 6.8σ level for Model 2. The two main consequences of this precise lens modelling are explained in more detail below. However, even if the lens parameters are recovered at sub-per cent precision, they have significantly different values in Models 1 and 2. Therefore, we would like to highlight that the recovered parameter values are precise, but at least one, and possibly both, of the models are incorrect descriptions of the data, as we discuss further in Section 4.2. In other words, whilst the statistical uncertainties are quite small, the systematic uncertainties due to our model choices may be up to two orders of magnitude larger.

Our findings demonstrate that high-resolution and high-sensitivity observations are vital for testing complex mass models, as opposed to the standard assumption of a smooth power-law elliptical mass density distribution. However, such in-depth studies of the global mass distribution using mas-resolution observations have

been mainly performed on compact (lensed) radio cores, as a result of the selection criterium of most lensing surveys at radio wavelengths (e.g. JVAS/CLASS; Browne et al. 2003; Myers et al. 2003). For example, using the same parametric lens modelling method applied here, the precision on the mass model parameters for CLASS B0712+472 is of the order of ~10 per cent when using the positions of the four compact lensed images measured with the VLBA at 1.7 GHz (10 mas FWHM beam size) as constraints (Hsueh et al. 2017). Also, by using the position of the four images of CLASS B1555+375, as measured with MERLIN at 5 GHz (50 mas FWHM beam size), the precision of the mass model parameters is of the order of ~20 per cent (Hsueh et al. 2016).

As such, previous studies of the mass distributions of gravitational lenses have been limited by the number of systems that show extended structure, either due to the intrinsic source morphology
<figure><img src="image_8.png" /><figcaption>**Figure 8.** The marginalized posterior PDF for the lens model parameters of Model 2. The contours enclosing the 99, 95, and 68 percentiles indicate the distribution between two parameters of the lens model. The PDF of each parameter is shown at the top of each column. The meaning of the parameters, the maximum-likelihood model values for each parameter, and their uncertainties, are presented in Table 2.</figcaption></figure>

or low brightness of any extended emission. To improve on this, we have started a high-sensitivity VLBI campaign of a carefully selected sample of gravitationally lensed radio sources with radio-bright Einstein rings or potentially extended arcs. The most promising sources have been followed-up with global VLBI imaging at mas resolution and high sensitivity, and will be presented in forthcoming papers. Nevertheless, the next generation of interferometers (i.e. Square Kilometer Array; SKA) will allow the discovery of ~$10^5$ gravitational lenses with both compact and extended structure, increasing by several orders of magnitude the number of systems suitable for testing mass distributions on mas-scales (Koopmans, Browne & Jackson 2004; McKean et al. 2015). In addition, the next generation of optical/infrared telescopes will have a sensitivity and angular resolution that is comparable to global VLBI observations,

for example the European-Extremely Large Telescope (E-ELT), which will allow a multi-wavelength test of lensing macro models at radio and infrared wavelengths.

## 4.2 Evidence for additional mass structure

In Fig. 6 we show the image position residuals for Models 1 and 2. We find that the image position residuals do not correlate with any particular group of lensed images, but they are scattered almost uniformly in RA and Dec. A possible explanation for the offset between the observed and model-predicted positions is the presence of some additional mass structure within the lensing galaxy or along
the line of sight that has not yet been taken into account in the lens
models presented here.

Based on the methodology provided by Despali et al. (2017), we find that the combination of source and lens redshifts for MG J0751+2716 leads to a projected number density of low-mass ($M \sim 10^6 M_{\odot}$) line-of-sight haloes per arcsec$^2$ of the order of a few tens for a CDM scenario. This additional mass can be contained, for example, in sub-haloes associated with the G3 dark matter halo, which can change the deflection angle and, therefore, shift the position of the lensed images from what is expected from a smooth mass distribution (Wambsganss & Paczynski 1992; Metcalf & Madau 2001; Metcalf 2002; Inoue & Chiba 2003, 2005b, 2005a; Sluse et al. 2012). This method has been used to quantify the level of sub-structures in the intermediate mass regime ($\sim 10^8 M_{\odot}$) from adaptive optics imaging of extended gravitational arcs (Vegetti et al. 2012) and from spectro-imaging of the narrow-line region of lensed quasars (Nierenberg et al. 2014). However, only VLBI observations can directly resolve the small-scale astrometric shifts due to very low mass haloes ($\sim 10^6 M_{\odot}$; McKean et al. 2015) because of the excellent astrometric information provided by the data (Chen et al. 2007; Keeton & Moustakas 2009).

For example, VLBI observations of CLASS B0128+437 at mas
resolution revealed astrometric offsets of between 5 and 10 mas
(much larger than the intrinsic astrometric precision) that have been
ascribed to the presence of substructure within the main lensing
galaxy (Biggs et al. 2004). Also, multifrequency global VLBI obser-
vations of MG B2016+112 confirmed that the astrometric anomaly
observed for this system could be entirely attributed to a luminous
satellite associated with the lensing galaxy (Koopmans et al. 2002;
More et al. 2009). In the case of MG J0414+0534, the high reso-
lution imaging from VLBI was used to infer the position and mass
of a dark sub-halo (MacLeod et al. 2013). Moreover, the radio-loud
lensing systems CLASS B1933+503, CLASS B1555+375, and
CLASS B0712+472 have extended arcs, which also show hints of
a disturbed morphology or anomalous flux ratios at radio wave-
lengths (Jackson et al. 1998; Marlow et al. 1999; Norbury et al.
2001). Subsequently, optical imaging of these systems revealed
that the lensing galaxy is a late-type galaxy and by adding the disc
as the additional mass component the position of the lensed images
could be completely recovered (Suyu et al. 2012; Hsueh et al. 2016,
2017).

Nevertheless, given the complexity of the mass model for this
lensing system, we cannot draw stringent conclusions on the origin
of the observed astrometric anomaly. Only a Bayesian grid-based
analysis that takes into account the flux density distribution of the
entire lensed arcs and performs lens-potential corrections directly
in the visibility plane can test whether our parametric models are
too simplistic for this system, or if there is the need for extra-mass
in the model.

**4.3 Evidence in favour of the two-phase galaxy formation scenario**

Both Model 1 and Model 2 find that the mass density profile for G3 is steeper than isothermal with $\gamma_1 = 2.08 \pm 0.02$ (4.2$\sigma$ level) and $\gamma_2 = 2.16 \pm 0.02$ (6.8$\sigma$ level), at the Einstein radius $b_1 = 0.4025 \pm 0.0008$ and $b_2 = 0.307 \pm 0.002$ arcsec, respectively. These density slopes are consistent within $1\sigma$ of the distribution of slopes from the Sloan Lens ACS Survey (SLACS) sample of lenses (Auger et al. 2010), which on average is $\gamma = 2.078 \pm 0.027$, therefore mildly steeper than isothermal (Koopmans et al. 2009; Auger et al. 2010; Barnabè et al. 2011; Sonnenfeld et al. 2013).

Therefore, the MG J0751+2716 lensing galaxy has a similar den-
sity profile when compared to other early-type galaxies that act as
strong gravitational lenses at a similar redshift.

Nevertheless, this density profile slope can be considered as ev-
idence for the so-called *two-phase* galaxy formation scenario for
G3, which is a low-mass early-type galaxy (with velocity disper-
sion $\sigma_v = 101$ km s$^{-1}$; Alloin et al. 2007) in a very rich environment
(Fig. 3, but also see Tonry 1998; Momcheva et al. 2006; Alloin et al.
2007; Momcheva et al. 2015; Wilson et al. 2016). It has been shown
that early-type satellite gravitational lenses within groups or clus-
ters of galaxies (as in the case of G3) are better modelled with
mass density profiles that are steeper than isothermal (Rusin et al.
2002; Dobke, King & Fellhauer 2007; Auger et al. 2008). Also, the
SLACS lenses show a superisothermal mass density profile when
associated with a perturbing companion galaxy (in this case it is
likely G1), suggesting that this steepening can be attributed to their
possible interaction (Auger 2008). Moreover, from observations and
simulations it was found that low-mass (and compact) early-type
galaxies have a mass slope $\gamma$ steeper than isothermal, while high-
mass early type galaxies have a shallower $\gamma$ (Barnabè et al. 2011;
Sonnenfeld et al. 2012; Remus et al. 2013; Dutton & Treu 2014;
Tortora et al. 2014).

The theoretical scenario for explaining $\gamma > 2$ in low-mass non-
isolated galaxies (called *two-phase* scenario) regards the interaction
with their companion galaxies, in particular, their merger history
(Guo & White 2008; Oser et al. 2010; Johansson, Naab & Ostriker
2012; Remus et al. 2013; Dubois et al. 2013). In this framework,
at the early stages of galaxy formation gas-rich mergers lead to an
enhanced *in situ* star formation; the dissipative process due to *in situ*
star formation cause a mass density profile that is superisothermal
with an increased baryonic matter content in the central regions
of the galaxy. Then, after $z \sim 2$, the merging events lead to a
reordering of the early-type galaxy into an isothermal state, because
the dissipative processes are not dominant anymore, and the galaxy
growth is principally driven by gas-poor mergers. The data in hand
for MG J0751+2716 are consistent with this model, given the dense
environment of the system and a density profile that is steeper than
the isothermal case.

**5 CONCLUSIONS**

We have presented sensitive global VLBI observations at 1.65 GHz of the radio-loud quasar MG J0751+2716, which is gravitationally lensed by a foreground group of galaxies to produce gravitational arcs that are extended by 200–600 mas; these data represent the highest angular resolution imaging of extended gravitational arcs from a gravitational lens. Our observations demonstrate that mas resolution observations of gravitationally lensed radio sources can provide a large number of constraints to the lensing mass model, which can be used to search for any deviation from a globally smooth mass distribution. By using the positions of four quadruply imaged components and two doubly imaged components, we investigate a simple single-lens model and a more realistic mass model for the group of galaxies that is associated with the main lensing galaxy. We find that from these constraints, we are able to infer the lens model parameters with a precision of less than a per cent, even though our models are not accurate enough to fit the positions of the observed images to the measurement error level. Furthermore, both models suggest an inner density slope for the main lensing galaxy that is steeper than isothermal. This is consistent with studies of other low-mass early-type satellite galaxies in dense environments, and is in agreement with the *two-phase* galaxy formation scenario.
This is important, because more than 50 per cent of galaxies are found to lie in groups, at least locally, and to date, there is not a complete picture of the total projected mass distribution of galaxy groups.

Due to the excellent sensitivity and high angular resolution of the VLBI imaging, we find there is a discrepancy between the observed and predicted positions of the lensed images, with an average position rms of the order of 3 mas for the simple parametric models tested here. At this stage, it is not clear if these deviations are due to some additional mass structure in the form of a population of low mass sub-haloes that are either part of the lensing group or along the line of sight, or if the complexity of the group environment is not being fully taken into account by the parametric models. In a future paper, we will present modelling with a grid-based source surface brightness distribution, which is fitted directly with the visibility data, that will allow the complete set of extended gravitational arcs to be used as constraints. In addition, with the improved source model, we will be able to test non-parametric lens models using grid-based corrections to the gravitational potential, based on the methodology of Vegetti & Koopmans (2009), which will shed light on the cause of the astrometric anomaly seen in the compact lensed components of MG J0751+2716.

## ACKNOWLEDGEMENTS

LVEK is supported through an NWO-VICI grant (project number 639.043.308). CDF acknowledges support from the U.S. National Science Foundation, via grant AST-1715611. This project has received funding from the European Research Council under the European Union's Horizon 2020 research and innovative programme (grant agreement No 758853). DJL acknowledges support from the European Research Council (ERC) starting grant 336736CALENDAS. The National Radio Astronomy Observatory is a facility of the National Science Foundation operated under cooperative agreement by Associated Universities, Inc. The European VLBI Network is a joint facility of European, Chinese, South African and other radio astronomy institutes funded by their national research councils. Scientific results from data presented in this publication are derived from the following EVN project code: GM070. Based on observations made with the NASA/ESA Hubble Space Telescope, obtained from the Data Archive at the Space Telescope Science Institute, which is operated by the Association of Universities for Research in Astronomy, Inc., under NASA contract NAS 5-26555. These observations are associated with program 7495.

## REFERENCES

* Alloin D., Kneib J.-P., Guilloteau S., Bremer M., 2007, A&A, 470, 53

* Auger M. W., 2008, MNRAS, 383, L40

* Auger M. W., Fassnacht C. D., Abrahamse A. L., Lubin L. M., Squires G. K., 2007, AJ, 134, 668

* Auger M. W., Fassnacht C. D., Wong K. C., Thompson D., Matthews K., Soifer B. T., 2008, ApJ, 673, 778

* Auger M. W., Treu T., Bolton A. S., Gavazzi R., Koopmans L. V. E., Marshall P. J., Moustakas L. A., Burles S., 2010, ApJ, 724, 511

* Barnabè M., Czoske O., Koopmans L. V. E., Treu T., Bolton A. S., 2011, MNRAS, 415, 2215

* Bertin E., Arnouts S., 1996, A&AS, 117, 393

* Biggs A. D., Browne I. W. A., Muxlow T. W. B., Wilkinson P. N., 2001, MNRAS, 322, 821

* Biggs A. D., Browne I. W. A., Jackson N. J., York T., Norbury M. A., McKean J. P., Phillips P. M., 2004, MNRAS, 350, 949

* Browne I. W. A. et al., 2003, MNRAS, 341, 13

* Burke B. F., 1990, in Mellier Y., Fort B., Soucail G., eds, Lecture Notes in Physics. Vol. 360, Gravitational Lensing. Springer Verlag, Berlin, p 127,

* Chen J., Rozo E., Dalal N., Taylor J. E., 2007, ApJ, 659, 52

* Chen G. C.-F. et al., 2016, MNRAS, 462, 3457

* Cornwell T. J., 2008, IEEE J. Sel. Top. Signal Process., 2, 793

* Dalal N., Kochanek C. S., 2002, ApJ, 572, 25

* Despali G., Vegetti S., White S. D. M., Giocoli C., van den Bosch F. C., 2017, MNRAS, 475, 5424

* Dobke B. M., King L. J., Fellhauer M., 2007, MNRAS, 377, 1503

* Dubois Y., Gavazzi R., Peirani S., Silk J., 2013, MNRAS, 433, 3297

* Dutton A. A., Treu T., 2014, MNRAS, 438, 3594

* Falco E. E. et al., 1999, preprint (arXiv Astrophysics)

* Fassnacht C. D., Kocevski D. D., Auger M. W., Lubin L. M., Neureuther J. L., Jeltema T. E., Mulchaey J. S., McKean J. P., 2008, ApJ, 681, 1017

* Gavazzi R., Treu T., Marshall P. J., Brault F., Ruff A., 2012, ApJ, 761, 170

* Giocoli C., Marulli F., Baldi M., Moscardini L., Metcalf R. B., 2013, MNRAS, 434, 2982

* Guo Q., White S. D. M., 2008, MNRAS, 384, 2

* Hewitt J. N., 1992, in Kayser R., Schramm T., Nieser L., eds, Lecture Notes in Physics. Vol. 406, Gravitational Lenses. Springer Verlag, Berlin, p. 44,

* Hewitt J. N., Turner E. L., Schneider D. P., Burke B. F., Langston G. I., 1988, Nature, 333, 537

* Hsueh J.-W., Fassnacht C. D., Vegetti S., McKean J. P., Spingola C., Auger M. W., Koopmans L. V. E., Lagattuta D. J., 2016, MNRAS, 463, L51

* Hsueh J.-W. et al., 2017, MNRAS, 469, 3713

* Inoue K. T., Chiba M., 2003, ApJ, 591, L83

* Inoue K. T., Chiba M., 2005a, ApJ, 633, 23

* Inoue K. T., Chiba M., 2005b, ApJ, 634, 77

* Jackson N. et al., 1998, MNRAS, 296, 483

* Johansson P. H., Naab T., Ostriker J. P., 2012, ApJ, 754, 115

* Keeton C. R., 2001a, preprint (arXiv Astrophysics)

* Keeton C. R., 2001b, preprint (arXiv Astrophysics)

* Keeton C. R., Moustakas L. A., 2009, ApJ, 699, 1720

* Keeton C. R., Kochanek C. S., Seljak U., 1997, ApJ, 482, 604

* Keeton C. R., Christlein D., Zabludoff A. I., 2000, ApJ, 545, 129

* Kochanek C. S., 2002, in Natarajan P., eds, The Shapes of Galaxies and their Dark Halos. p. 62

* Koopmans L. V. E., Garrett M. A., Blandford R. D., Lawrence C. R., Patnaik A. R., Porcas R. W., 2002, MNRAS, 334, 39

* Koopmans L. V. E., Browne I. W. A., Jackson N. J., 2004, New A Rev., 48, 1085

* Koopmans L. V. E. et al., 2009, ApJ, 703, L51

* Lagattuta D. J., Auger M. W., Fassnacht C. D., 2010, ApJ, 716, L185

* Lagattuta D. J., Vegetti S., Fassnacht C. D., Auger M. W., Koopmans L. V. E., McKean J. P., 2012, MNRAS, 424, 2800

* Langston G. I. et al., 1989, AJ, 97, 1283

* Lehar J. et al., 1997, AJ, 114, 48

* MacLeod C. L., Jones R., Agol E., Kochanek C. S., 2013, ApJ, 773, 35

* Mao S., Schneider P., 1998, MNRAS, 295, 587

* Marlow D. R. et al., 1999, AJ, 118, 654

* McKean J. P. et al., 2005, MNRAS, 356, 1009

* McKean J. P. et al., 2007, MNRAS, 378, 109

* McKean J. P. et al., 2010, MNRAS, 404, 749

* McKean J. et al., 2015, Proceedings of Advancing Astrophysics with the Square Kilometre Array (AASKA14), p. 84

* Metcalf R. B., 2002, ApJ, 580, 696

* Metcalf R. B., Madau P., 2001, ApJ, 563, 9

* Momcheva I., Williams K., Keeton C., Zabludoff A., 2006, ApJ, 641, 169

* Momcheva I. G., Williams K. A., Cool R. J., Keeton C. R., Zabludoff A. I., 2015, ApJS, 219, 29

* More A., McKean J. P., Muxlow T. W. B., Porcas R. W., Fassnacht C. D., Koopmans L. V. E., 2008, MNRAS, 384, 1701

* More A., McKean J. P., More S., Porcas R. W., Koopmans L. V. E., Garrett M. A., 2009, MNRAS, 394, 174

* Muñoz J. A., Falco E. E., Kochanek C. S., Lehár J., McLeod B. A., Impey C. D., Rix H.-W., Peng C. Y., 1998, Ap&SS, 263, 51

* Myers S. T. et al., 2003, MNRAS, 341, 1
Navarro J. F., Frenk C. S., White S. D. M., 1997, ApJ, 490, 493

Nierenberg A. M., Treu T., Wright S. A., Fassnacht C. D., Auger M. W., 2014, MNRAS, 442, 2434

Norbury M. A., Jackson N. J., Browne I. W. A., Wilkinson P. N., Marlow D. R., 2001, in Brainerd T. G., Kochanek C. S., eds, ASP Conf. Ser.Vol. 237, Gravitational Lensing: Recent Progress and Future Go, Astron. Soc. Pac., San Francisco, p. 97

Offringa A. R. et al., 2014, MNRAS, 444, 606

Oguri M., Keeton C. R., Dalal N., 2005, MNRAS, 364, 1451

Oldham L. et al., 2017, MNRAS, 465, 3185

Oser L., Ostriker J. P., Naab T., Johansson P. H., Burkert A., 2010, ApJ, 725, 2312

Patnaik A. R., 1993, in Surdej J., Fraipont-Caro D., Gosset E., Refsdal S., Remy M., eds, Liege International Astrophysical Colloquia. Vol. 31, p. 311

Peng C. Y., Ho L. C., Impey C. D., Rix H.-W., 2010, AJ, 139, 2097

Planck Collaboration XIII, 2016, A&A, 594, A13

Remus R.-S., Burkert A., Dolag K., Johansson P. H., Naab T., Oser L., Thomas J., 2013, ApJ, 766, 71

Rich J. W., de Blok W. J. G., Cornwell T. J., Brinks E., Walter F., Bagetakos I., Kennicutt Jr. R. C., 2008, AJ, 136, 2897

Rusin D., Norbury M., Biggs A. D., Marlow D. R., Jackson N. J., Browne I. W. A., Wilkinson P. N., Myers S. T., 2002, MNRAS, 330, 205

Sluse D., Chantry V., Magain P., Courbin F., Meylan G., 2012, A&A, 538, A99

Sonnenfeld A., Treu T., Gavazzi R., Marshall P. J., Auger M. W., Suyu S. H., Koopmans L. V. E., Bolton A. S., 2012, ApJ, 752, 163

Sonnenfeld A., Treu T., Gavazzi R., Suyu S. H., Marshall P. J., Auger M. W., Nipoti C., 2013, ApJ, 777, 98

Sonnenfeld A., Treu T., Marshall P. J., Suyu S. H., Gavazzi R., Auger M. W., Nipoti C., 2015, ApJ, 800, 94

Suyu S. H. et al., 2012, ApJ, 750, 10

Tonry J. L., 1998, AJ, 115, 1

Tortora C., La Barbera F., Napolitano N. R., Romanowsky A. J., Ferreras I., de Carvalho R. R., 2014, MNRAS, 445, 115

Treu T., 2010, ARA&A, 48, 87

Vegetti S., Koopmans L. V. E., 2009, MNRAS, 392, 945

Vegetti S., Lagattuta D. J., McKean J. P., Auger M. W., Fassnacht C. D., Koopmans L. V. E., 2012, Nature, 481, 341

Vegetti S., Koopmans L. V. E., Auger M. W., Treu T., Bolton A. S., 2014, MNRAS, 442, 2017

Walsh D., Carswell R. F., Weymann R. J., 1979, Nature, 279, 381

Wambsganss J., Paczynski B., 1992, ApJ, 397, L1

Wilson M. L., Zabludoff A. I., Ammons S. M., Momcheva I. G., Williams K. A., Keeton C. R., 2016, ApJ, 833, 194

Winn J. N., Hewitt J. N., Schechter P. L., 2001, in Brainerd T. G., Kochanek C. S., eds, ASP Conf. Ser.Vol. 237, Gravitational Lensing: Recent Progress and Future Go, Astron. Soc. Pac., San Francisco, p. 61

Wucknitz O., Refsdal S., 2001, in Brainerd T. G., Kochanek C. S., eds, ASP Conf. Ser.Vol. 237, Gravitational Lensing: Recent Progress and Future Go, Astron. Soc. Pac., San Francisco, p. 157

This paper has been typeset from a TEX/LATEX file prepared by the author.