A million-solar-mass object detected at a cosmological distance using gravitational imaging

Observation

The VLBI observation of JVAS B1938+666 used for this work was performed using a global VLBI array combining antennas from the Very Long Baseline Array (VLBA) and European VLBI Network (EVN) at 1.7 GHz for 14 h at a data recording rate of 512 Mbits s⁻¹ (ID: GM068; Principal Investigator: McKean). The observing strategy followed a standard phase referencing mode only for the VLBA part (phase calibrator J1933+654). Throughout the observations, scans on the fringe finder 3C454.3 were also included every ~4 h for the bandpass calibration. The correlation of the data from the 19 antennas was performed at the Joint Institute for VLBI-European Research Infrastructure Consortium (JIV-ERIC) and resulted in a dataset with eight intermediate frequencies, each with 8-MHz bandwidth and divided into 32 channels for two circular polarizations (RR and LL).

The data were processed within the Astronomical Image Processing System (aips)⁴³ package following a typical calibration procedure for phase-referenced observations. The resulting deconvolved image was obtained with a restoring beam of 7.4 mas × 4.7 mas at a position angle of 32.1 degrees east of north³³. The noise, assumed to be Gaussian and uncorrelated, was measured directly in the Fourier domain by subtracting time-adjacent visibilities to remove the source signal, then computing the root mean square (r.m.s.) within 30-minute time intervals. The triangle of highly sensitive baselines between Effelsberg, Jodrell Bank and Westerbork were flagged to prevent them from dominating the model. The baselines Green Bank–Hancock, Green Bank–Owens Valley, Green Bank–Pie Town and Los Alamos–Pie Town were also flagged for several 1- to 2-h time intervals owing to strong radio frequency interference.

The surface-brightness distribution of the background radio source consists of two radio lobes (hot spots), the brighter of which is lensed into an extended gravitational arc, whereas the fainter one is doubly imaged (Fig. 1). The physical properties of the source are derived and discussed in ref. ³³.

Bayesian inference

We use the visibility-plane gravitational lens modelling software PRONTO^{34,35,36,37,38} to jointly infer the pixellated source surface-brightness vector s and lens parameters η_H from the observed visibilities d, as well as to obtain the Bayesian log-evidence for each model parametrization H. At each likelihood evaluation (the first level of inference), we obtain the maximum a posteriori source s_MP for a given η_H and source regularization weight λ_s by solving

where

Here \(L({{\boldsymbol{\eta }}}_{H})\) is the lens operator, which maps light from the source plane to the lens plane, \({\lambda }_{{\boldsymbol{s}}}R_{{\boldsymbol{s}}}\) is the source prior covariance, which enforces the lens equation by penalizing strong surface-brightness gradients and \(C^{-1}\) is the noise covariance of the data. We solve equation (1) using a preconditioned conjugate gradient solver, where the Fourier operator \(D\) is implemented using a non-uniform fast Fourier transform. We refer the reader to ref. ³⁷ for further details on the method.

In the second level of inference, we sample the lens parameters η_H and source regularization weight λ_s. The posterior is

The likelihood (which is the evidence from the source-inversion step) is

This expression follows from the marginalization over all possible sources s when the noise and source prior are both Gaussian. As \({{{{\boldsymbol{R}}}}}_{{\mathbf{s}}}\) and \({{{{\boldsymbol{C}}}}}^{-1}\) are sparse, the terms containing them are straightforward to evaluate. Computing \(\log \det A\) is non-trivial; we approximate it using the preconditioner from the inference on s_MP as described in ref. ³⁷. The χ² term is

We optimize its evaluation using the fast-χ² technique from ref. ³⁸.

We compared different model parametrizations, H, using the log-evidence \(\log {{\mathcal{E}}}_{{\rm{H}}}\equiv \log P({\boldsymbol{d}}| {\rm{H}})\), which was obtained by marginalizing over all parameters. We used the nested sampling algorithm MULTINEST⁴⁴ to compute the log-evidence and to sample the posterior distributions of η_H and λ_s.

Our method takes both noise and priors to be Gaussian. From a practical standpoint, these assumptions make our analysis computationally tractable by allowing the marginalization over all possible configurations of the source and pixellated potential corrections (‘Gravitational imaging’) using a single linear solution. Although it is the case that non-Gaussian noise statistics can, in principle, arise from calibration errors or correlated instrumental effects, we do not expect this to be a major issue given the high signal-to-noise ratio of the observation. Nevertheless, we tested our analysis on several independent sets of calibration solutions, including on individual spectral windows, and found our results to be robust across all cases. Most reassuringly, we found that the gap in the bright arc induced by object \({\mathcal{V}}\) cannot be suppressed by the self-calibration process.

Similarly, although there exist techniques for applying highly informative generative machine-learning models as non-Gaussian priors for ‘realistic’ sources and potentials (for example, ref. ⁴⁵), these are susceptible to systematic biases learned from the training data (which are themselves taken from numerical simulations). False-positive subhalo detections are a particularly relevant concern for this work. Instead, we opt for highly uninformative Gaussian priors, which are agnostic to the shapes of the source surface brightness and the linearized potential corrections, which thus maximizes the ability of the data to drive the inference.

Parametric lens model components

We now introduce the surface mass density profiles used for parametric modelling in this work. All quantities labelled κ are in units of the critical surface mass density for strong lensing,

where D_s, D_l and D_ls are the angular diameter distances from observer to source, from observer to lens and from lens to source, respectively. Assuming a Planck 2015 (ref. ⁴⁶) cosmology, \({\varSigma }_{{\rm{cr}}}=1.50\times 1{0}^{11}\,{M}_{\odot }\,{{\rm{arcsec}}}^{-2}\) for JVAS B1938+666 (z_s = 2.059, z_l = 0.881).

We model the smooth galaxy-scale mass distribution (the macromodel) as an elliptical power-law mass distribution (for example, ref. ⁴⁷), with a projected surface mass density of

which we define in terms of an effective Einstein radius R_E and the elliptical radius ξ² ≡ x²q + y²/q, where q is the axis ratio, and γ is the three-dimensional logarithmic slope, where γ = 2 is isothermal. The additional parameters x₀, y₀ and θ₀ set the position and orientation of the profile. We use FASTELL⁴⁸ to compute the corresponding deflection angles.

To account for non-ellipticity in the lens galaxy (for example, diskiness or boxiness), we include multipole perturbations of orders m = 3 and 4, parametrized as

which we write in circular polar coordinates here for readability. The coefficients a_m and b_m encode the amplitude and orientation of each multipole term and γ is fixed to that of the underlying elliptical power-law model (equation (7)). We use Gaussian priors for a_m and b_m with μ = 0 and σ = 0.01, a choice motivated by the results of numerical simulations of massive elliptical galaxies (for example, ref. ⁴⁹), which is also consistent with previous lens modelling results^23,38,50. In addition, we include an external shear term in the macromodel, with strength Γ and direction θ_Γ.

For parametric modelling of the low-mass perturbers \({\mathcal{A}}\) and \({\mathcal{V}}\) (Table 1), we use the following spherically symmetric mass profiles. Models PJ_free and PJ_tidal are truncated isothermal profiles (PJ; for example, refs. ^47,51) with a projected surface mass density of

where r_t is the truncation radius, x = r/r_t and m_tot is the total mass. We leave r_t as a free parameter for PJ_free, whereas for PJ_tidal we set r_t to be the tidal radius

Here R is the three-dimensional radius to the centre of the lens and M_lens (<R) is the mass of the main lens enclosed within R; see ref. ²⁰. Because we have no knowledge of the three-dimensional position of a low-mass perturber within the lens, we take R to be the projected distance to the perturber in the plane of the lens.

In the case of NFW_sub, we model perturber \({\mathcal{A}}\) as a spherical NFW profile^47,52, with a surface mass density of

where x = r/r_s, κ_s = ρ_sr_s/Σ_cr and

We found that using a mass-concentration-redshift relationship for CDM haloes (for example, ref. ⁵³) leads to virial concentrations that are much too low to produce the required lensing effect for a halo of a given virial mass (see also refs. ^20,23). Therefore, we leave the concentration of NFW_sub as a free parameter in this work. The redshift is fixed to that of the main lens.

GI

GI (refs. ^{17,24,34,39,54}) is a technique for recovering pixellated (non-parametric) corrections to the lensing convergence. GI is well suited for finding local overdensities that have not been accounted for in the parametric mass model, η_H, which remains fixed during the GI procedure. Our GI formulation augments the lens operator (for example, equation (1)) with an extra block column representing linearized degrees of freedom in the lensing potential:

where

is the matrix of source gradients and

is the matrix of derivatives of the deflection angles α with respect to the pixellated lensing potential ψ. Operator \(L_{{\rm{GI}}}\) acts on the concatenated vector of the pixellated source and linearized potential corrections

Similarly, we define the block regularization matrix

For this work, we apply three different forms for the potential regularization operator \({R}_{{\boldsymbol{\psi }}}\), checking them against one another to ensure the robustness of our results. These penalize total mass, surface density gradients or surface density curvature. In each case, \(R_{{\boldsymbol{\psi }}}\) is assembled using standard finite-difference operators for the gradient and Laplacian on a Cartesian grid.

It can be shown that the proper linearization over ψ yields a linear system analogous to equation (1):

Formally, equation (18) represents a single Newton iteration in the space of s and ψ. After iterating until convergence, it can be shown that the likelihood with respect to the regularization weights, \(\log P({\mathbf{d}}| {\lambda }_{{\mathbf{s}}},{\lambda }_{{\boldsymbol{\psi }}})\) is obtained by substituting s → r, \({\lambda }_{{\mathbf{s}}}{R}_{{\mathbf{s}}}\to {R}_{{\rm{GI}}}\) and \({L}\to {R}_{{\rm{GI}}}\) in equation (4). A single likelihood evaluation in our GI procedure thus consists of the following.

1. (1)

   Initialize r = 0.

2. (2)

   Solve for Δr_MP.

3. (3)

   Update r = r + Δr_MP.

4. (4)

   Repeat steps 2 and 3 until convergence. We find that 10 iterations are generally sufficient for convergence; in practice, we use 20.

5. (5)

   Evaluate \(\log P({\mathbf{d}}| {\lambda }_{{\mathbf{s}}},{\lambda }_{{\boldsymbol{\psi }}})\) using equation (4).

Note that our GI implementation differs from that of ref. ⁵⁴, who used a single iteration per likelihood evaluation. Our approach of iterating to convergence allows the GI algorithm to more easily recover compact features in the lensing potential.

We use a standard simplex optimization algorithm to maximize the likelihood with respect to the regularization weights λ_s and λ_ψ, where each likelihood evaluation requires 20 subiterations as described above. We used starting values for λ_s and λ_ψ of 10⁹ and 10¹³, respectively, with initial logarithmic step sizes of 0.5 dex. After optimizing for λ_s and λ_ψ, we convert ψ into corrections to the lensing convergence (surface mass density) by means of \({\kappa }_{{\rm{GI}}}=-\frac{1}{2}{\nabla }^{2}{\boldsymbol{\psi }}\).

Our use of Gaussian priors for the linearized potential corrections was motivated by the need for an uninformative and computationally fast regularization term. A natural consequence of this choice is that ψ and hence κ_GI, prefer to resemble Gaussian random fields in regions where the data are uninformative. To aid in our interpretation of the GI results, we apply a 3σ_GI threshold to the convergence maps, treating any convergence over the threshold as a real feature. We compute σ_GI from the r.m.s. of the convergence corrections within the mask used for modelling the lensed images. The value of σ_GI must be empirically estimated in this way, as the operator-based iterative solver framework (‘Bayesian inference’) precludes direct manipulation of the posterior covariance matrix.

Expected number counts of detectable subhaloes

We express the differential mass function for WDM subhaloes in terms of the ‘half-mode mass’ M_hm, as

with α = −1.9, α₂ = 1.1, β = 1.0 and γ = −0.5 (refs. ^55,56). The mass m is the total mass of a PJ subhalo. The expected number of subhaloes in the mass range \({M}_{\min }=1{0}^{6}\,{{M}}_{\odot }\) to \({M}_{\max }=1{0}^{7}\,{{M}}_{\odot }\) is then

where R_E is the Einstein radius of the lens and M_lens(<2R_E) is the total projected mass of the lens within twice the Einstein radius. The variable f_sub is the fraction of dark matter (normalized with respect to the CDM mass function) contained in subhaloes; we use f_sub = 0.012 (ref. ⁵⁷). The variable A_sens is the area around the lensed arcs that we expect to be sensitive to the presence of subhaloes, which we obtain by thresholding the deconvolved image at 5σ above the noise for a total area of \({A}_{{\rm{sens}}}=1.06\times 1{0}^{-2}\,{{\rm{arcsec}}}^{2}\). This area roughly corresponds to one primary beam width (~5 mas) along each bright arc. Defining A_sens in such a tight region around the brightest parts of the lensed arcs is an intentionally conservative choice made to ensure that all subhaloes ≳10⁶ M_⊙ within this region have been detected. Our use of a constant sensitivity within A_sens is a rough but conservative approximation to more sophisticated sensitivity mapping techniques (for example, refs. ^28,56).

Parametric modelling: macromodel and \({\mathcal{A}}\) only

The results of the fully parametric modelling are summarized in Table 1. We first consider those fully parametric models consisting of the macromodel and the previously detected perturber \({\mathcal{A}}\) (‘Parametric lens model components’). Because \({\mathcal{A}}\) lies ~50 mas away from the nearest lensed radio emission, its redshift cannot be robustly constrained by our VLBI observation alone. For this work, we therefore assume that \({\mathcal{A}}\) lies at the redshift of the main lens.

The best parametric models for \({\mathcal{A}}\), in terms of Bayesian log-evidence \(\Delta \log {\mathcal{E}}\), are truncated isothermal subhalo models (PJ), with either free or fixed (tidal) truncation radii being equally preferred (\({\mathcal{A}}=\) PJ_free or \({\mathcal{A}}=\) PJ_tidal; \(\Delta \log {\mathcal{E}}=16\)). We additionally tested an NFW profile with free virial concentration (\({\mathcal{A}}=\)NFW_sub; \(\Delta \log {\mathcal{E}}=13\)) as well as a completely smooth macromodel (\({\mathcal{A}}=\)None; \(\Delta \log {\mathcal{E}}\equiv 0\)). For simplicity, we chose the best model with fewer free parameters, \({\mathcal{A}}=\) PJ_tidal, as our fiducial model for the GI procedure. We found that the choice of profile \({\mathcal{A}}\) does not affect the inferred properties of \({\mathcal{V}}\).

GI detection of \({\mathcal{V}}\)

The results of the GI modelling, applied to the parametric model with \({\mathcal{A}}=\) PJ_tidal, are shown in Fig. 1 and Supplementary Fig. 1. For three different regularization types (penalizing the convergence, gradient of convergence or curvature of convergence) and two different grid resolutions (N_GI = 512 and N_GI = 1,024, corresponding to pixel sizes of Δx_GI = 3.5 mas and Δx_GI = 1.8 mas, respectively), we consistently recover a compact density correction well above the 3σ_GI threshold used to identify bona fide features in the surface mass density that were unaccounted for during the initial parametric modelling. We label this feature \({\mathcal{V}}\), as it was discovered using a VLBI observation of JVAS B1938+666.

Defining \({m}_{80,{\mathcal{V}}}\) to be the cylindrical mass contained within a projected radius of 80 pc on the lens plane (‘Parametric modelling: macromodel, \({\mathcal{A}}\) and \({\mathcal{V}}\)’), we found \(8.3\times 1{0}^{5}\,{{M}}_{\odot }\le {m}_{80,{\mathcal{V}}}\le 1.8\times 1{0}^{6}\,{{M}}_{\odot }\) for the six GI runs. The scatter in GI masses is due to the large number of degrees of freedom in the pixellated convergence map, combined with the GI technique’s lack of prior knowledge on the approximate sphericity of gravitationally bound astrophysical objects. GI is an important technique for identifying and estimating the masses of positive density corrections to a parametric lens model; however, it is necessary to verify detection \({\mathcal{V}}\) with independent parametric modelling.

Parametric modelling: macromodel, \({\mathcal{A}}\) and \({\mathcal{V}}\)

We repeat the fully parametric modelling procedure, this time also including a parametric profile for perturber \({\mathcal{V}}\). The results are summarized in Table 1 and Supplementary Fig. 3. We assume, for this work, that perturber \({\mathcal{V}}\) lies at the lens redshift of z = 0.881, using a truncated isothermal profile (PJ; ‘Parametric lens model components’). We find that, regardless of the choice of profile for detection \({\mathcal{A}}\), the inclusion of a parametric model for \({\mathcal{V}}\) increases the log-evidence by at least \(\Delta \log {\mathcal{E}} > 332\) (>25σ) in all cases. Models \({\mathcal{A}}=\) PJ_free, \({\mathcal{V}}=\) PJ_free and \({\mathcal{A}}=\) PJ_tidal, \({\mathcal{V}}=\) PJ_free are equally preferred. As before, we chose \({\mathcal{A}}=\) PJ_tidal, \({\mathcal{V}}=\) PJ_free with fewer free parameters for a more detailed discussion, and treat it as our ‘best model’, with a significance of \(\Delta \log {\mathcal{E}}=364\) over the smooth macromodel and \(\Delta \log {\mathcal{E}}=348\) (26σ) over \({\mathcal{A}}=\) PJ_tidal, \({\mathcal{V}}=\)None. For the rest of the discussion, it is implied that \({\mathcal{A}}=\) PJ_tidal.

In the \({\mathcal{V}}=\) PJ_free parametrization, \({\mathcal{V}}\) has a total mass of \({m}_{{\mathcal{V}}}=(2.82\pm 0.26)\times 1{0}^{6}\,{M}_{\odot }\) and a truncation radius \({r}_{\rm{t},{\mathcal{V}}}=149\pm 18\,{\rm{pc}}\). To express the inferred mass in a way that is independent of \({r}_{\rm{t},{\mathcal{V}}}\), we found 80 pc to be the radius at which the enclosed (projected) mass \({m}_{80,{\mathcal{V}}}\) is uncorrelated with \({r}_{\rm{t},{\mathcal{V}}}\). For \({\mathcal{V}}=\) PJ_free, \({m}_{80,{\mathcal{V}}}=(1.13\pm 0.04)\times 1{0}^{6}\,{M}_{\odot }\), which is consistent with the gravitational imaging results (Fig. 3). For \({\mathcal{V}}=\) PJ_tidal, \({m}_{80,{\mathcal{V}}}=(1.07\pm 0.04)\) agrees with \({\mathcal{V}}=\) PJ_free, which indicates that the enclosed mass at 80 pc is indeed particularly well constrained by the data. The position of \({\mathcal{V}}\) is constrained to 194 μas precision in the y direction and 86 μas in the x direction. For all models, \({m}_{80,{\mathcal{V}}}\), \({x}_{{\mathcal{V}}}\) and \({y}_{{\mathcal{V}}}\) were consistent within the 1σ uncertainties regardless of the choice of profile \({\mathcal{A}}\).

The tidal truncation radius of \({\mathcal{V}}=\) PJ_tidal (equation (10)) is \({r}_{\rm{t},{\mathcal{V}}}=53\pm 1\,{\rm{pc}}\), which is a factor of three smaller than that inferred using PJ_free. Relative to \({\mathcal{V}}=\) PJ_free, \({\mathcal{V}}=\) PJ_tidal is disfavoured by \(\Delta \log {\mathcal{E}}=-16\), which suggests that this smaller value of \({r}_{\rm{t},{\mathcal{V}}}\) is too compact to produce the required lensing effect. This is likely to be a result of the use of the projected two-dimensional separation as a proxy for the three-dimensional distance between the lens galaxy and the perturber; any offset of \({\mathcal{V}}\) along the line of sight would increase the three-dimensional distance and hence r_t, alleviating this issue.

Unlike the Keck AO observation of JVAS B1938+666, in which perturber \({\mathcal{A}}\) lies on top of the infrared arc, the projected lens-plane distance between \({\mathcal{A}}\) and the nearest lensed radio emission is 400 pc; therefore, we expect only the cylindrical mass out to a projected radius of 400 pc, which we define as \({m}_{400,{\mathcal{A}}}\), to be well constrained by the VLBI data. For \({\mathcal{A}}=\) PJ_tidal, we find \({m}_{400,{\mathcal{A}}}=(5.0\pm 0.8)\times 1{0}^{7}\,{M}_{\odot }\) and \({r}_{\rm{t},{\mathcal{A}}}=243\pm 20\,{\rm{pc}}\). For \({\mathcal{A}}=\) PJ_free, \({m}_{400,{\mathcal{A}}}=(5.2\pm 0.7)\times 1{0}^{7}\,{M}_{\odot }\) and \({r}_{\rm{t},{\mathcal{A}}}=15\pm 20\,{\rm{pc}}\). Hence, \({m}_{400,{\mathcal{A}}}\) is well constrained and consistent between the two, whereas \({r}_{\rm{t},{\mathcal{A}}}\) is unconstrained and simply reflects the tidal radius relationship and the log-uniform prior (respectively) used for these profiles. For comparison, Despali et al. (ref. ²⁰ and subsequent private communication) find a mass of \({m}_{400,{\mathcal{A}}}=(7.7\pm 0.1)\times 1{0}^{7}\,{{M}}_{\odot }\) for their best subhalo model of the Keck AO observation.