A Geometric Probe of Cosmology – II. Gravitational Lensing Time Delays and Quasar Reverberation Mapping Revisited

The time delay due to geometric and gravitational differences in the light paths between multiple images of strongly gravitationally lensed quasars is a direct measurement of cosmological distances, and therefore of cosmological parameters (Refsdal, 1964; Blandford & Narayan, 1992). Time delay cosmography is now an established independent method by which cosmological information can be determined, e.g. Tewes et al. (2013); Courbin et al. (2018); Birrer et al. (2019), and Suyu et al. (H0LiCOW; 2017) with a current 2.4 per cent precision measurement of $H_{0}$ (Wong et al., 2020). However, gravitational lens modelling depends implicitly upon the assumed underlying mass distribution. In particular, the “mass sheet degeneracy” involving the degeneracy of ill-constrained lens models leads to systematic uncertainties in estimations of $H_{0}$ (Falco et al., 1985; Saha, 2000; Kochanek, 2020; Birrer et al., 2020; Chen et al., 2020).

In a first paper (Ng & Lewis, 2020) we presented a geometric probe of cosmology independent of the lensing potential, by considering differential time delays over images, originating from spatially-separated photometric signals within a strongly lensed quasar. In the original work, we presented the predicted signal integrated across all wavelengths as a function of time, from considering a special case of a thin face-on ring geometry of the Broad Line Region (BLR) of the quasar. In this paper we re-examine the methodology and demonstrate some significant difficulties, most notably under-determination. We also consider analytically the picture in the line-of-sight velocity (i.e. spectral) and time delay space, referred to as a “velocity-delay map”; as well as the signals obtained at any spectral or temporal slice for axisymmetric BLR geometries with the inclined ring or disk, spherical shell, and double cone as examples.

The paper is organised as follows: we briefly review the definitions and the background in Section 2 and give an overview of the methodology and difficulties in Section 3. We look at the effect of the differential lensing on the velocity-delay maps of different BLR geometries in Section 4 and demonstrate the insensitivity of the required BLR parameter estimation to the differential lensing in Section 5.

An observer sees a photometric signal from an unlensed point source at line-of-sight physical distance $\xi_{z}$ and position $\bm{\beta}$ on the sky, at time $t=t_{i}$. When the source is lensed, the observer sees this signal in a given lensed image X at a position $\bm{\theta}_{X}$ on the sky and at time $t=t_{i}+\tau_{X}$, where $\tau_{X}$ is the time delay (e.g. Schneider et al., 1992; Blandford & Narayan, 1986, 1992):

Here $D=\frac{D_{d}D_{s}}{D_{ds}}\propto H_{0}^{-1}$ is the “lensing distance” or more commonly the “time-delay distance”, a ratio of angular diameter distances (subscript $d$, $s$, $ds$ denoting the angular diameter distance from the observer to the lens, to the source and from the lens to the source, respectively). The lensing distance is thus the factor containing all the cosmological information. We denote the speed of light by $c$, and the redshift of the lensing mass by $z_{d}$. The two-dimensional vector positions of image X in the lens plane, and the source in the source plane, are given by $\bm{\theta}_{X}$ and $\bm{\beta}$ respectively; scaled such that their magnitudes are the observed angular positions relative to the observer-lens axis. The dimensionless lensing potential is denoted by $\psi$.

Although the time delay $\tau$ is not an observable quantity, the time delays between images of strongly lensed quasars may be measured and is a conventional method employed to test cosmology. The observed time delay between two images A and B are given by

Since Equation (2) is dependent on the dimensionless lensing potential $\psi$ however, time delay measurements are limited by the assumptions and accuracy of the lens model. Simple gravitational lens systems are rare, and observational constraints on the lens model are limited due to the existence, typically, of only two or four images per system (Schneider & Sluse, 2013; Birrer et al., 2019).

In Paper I (Ng & Lewis, 2020), we instead consider lensed quasars as a finite rather than point source in combination with reverberation mapping (Blandford & McKee, 1982; Peterson, 1993; Cackett et al., 2021) to determine quasar structure. Broad emission lines are a feature of quasar spectra and correspond to gas clouds, known as the Broad Line Region (BLR), surrounding the central emitting accretion disk at some distance. Photons travelling outwards from the central source are absorbed and re-emitted by the BLR gas. The broad emission lines therefore respond to variations in the continuum luminosity of the central source with a time delay determined by the BLR geometry.

For an unlensed quasar, if the continuum emission (across all wavelengths) from the central source occurs at $t=t_{i}$ (i.e. a Dirac delta signal in time), then the BLR emission occurs at $t=t_{i}^{\prime}=t_{i}+\tau_{\textsc{rm}}$, where the reverberation mapping time delay is

Here $\hat{r}_{z}$ is the unit vector towards the observer and $\bm{r}=(r_{x},r_{y},r_{z})^{T}=(D_{s}\delta\bm{\beta},\delta\xi_{z})^{T}$ is the spatial position of a given BLR cloud relative to the central source. The cosmological redshift factor $(1+z_{s})$ corrects the time delay in the quasar rest frame for the time delay as measured in our observer frame. If we hold $\tau_{\textsc{rm}}$ fixed, then the “iso-delay surface” is a paraboloid of revolution around $\hat{r}_{z}$. Equation (3) therefore describes a family of nested paraboloids parametrised by $\tau_{\textsc{rm}}$. For the unlensed quasar, the time at which the flux from the BLR is seen is determined by the intersection of the BLR geometry with the iso-delay paraboloids.

Considering the effect of the lensing time delay, the time that we see the continuum signal, originating from the centre of the source at $\bm{\beta}$, in image X is given by

The time that we see a signal in image X originating from a BLR cloud at position $\bm{r}$, using the shorthand notation ${\tau}^{\prime}_{X}\equiv\,\tau(\bm{\theta}_{X}+\delta\bm{\theta}_{X},\bm{\beta}+\delta\bm{\beta},\xi_{z}+\delta\xi_{z})=\tau_{X}+\delta\tau_{X}$, is similarly given by

The differential lensing time delay within image X, may be found via a Taylor expansion to first order in $\delta\bm{\theta}_{X}$ and $\delta\bm{\beta}$ (the term corresponding to a line of sight displacement $\delta\xi_{z}$ may be omitted as it is many orders of magnitude less than that corresponding to a shift in the source plane; see Ng & Lewis (2020); Tie & Kochanek (2018)), using the lens equation $\bm{\nabla}_{\bm{\theta}}\psi=\bm{\alpha}_{X}=\bm{\theta}_{X}-\bm{\beta}$ where $\bm{\alpha}_{X}\equiv\bm{\alpha}(\bm{\theta}_{X})$ is the scaled deflection angle, as

The interval $\mathcal{T}_{X}(\bm{r})$ between the signal from a BLR cloud at $\bm{r}$ and the continuum from within a given image X is given by

where $\bm{e}_{X}\equiv\left(\frac{(1+z_{d})}{(1+z_{s})}\hat{\alpha}_{X,x},\frac{(1+z_{d})}{(1+z_{s})}\hat{\alpha}_{X,y},1\right)^{T}$ and the physical deflection angle is $\hat{\bm{\alpha}}_{X}\equiv\hat{\bm{\alpha}}(D_{d}\bm{\theta}_{X})=\tfrac{D_{d}}{D_{ds}}\bm{\alpha}(\bm{\theta}_{X})$. The magnitude of $\bm{e}_{X}$ is the eccentricity $e_{X}$ of the conic section which generates the quadric surface of revolution described by Equation (7). We regain the unlensed case, i.e. $\tau_{\textsc{rm}}$, simply by setting the deflection angle to $0$. The effect of differential lensing is therefore to deform each iso-delay paraboloid into a one of the sheets of a two-sheeted hyperboloid of revolution about the $\hat{e}_{X}$ axis, since $e_{X}>1$ when the quasar is lensed and $e_{X}=1$ when the quasar is unlensed. For the lensed quasar, the time at which the flux for the BLR is seen is determined by the intersection of the BLR geometry with a family of iso-delay hyperboloids.

The difference in the differential lensing time delays of an image pair is independent of $\bm{\beta}$ and $\psi(\bm{\theta}_{X})$ to first order in $\delta\bm{\beta}$ and $\delta\bm{\theta}_{X}$:

In general, we will write $\mathcal{T}_{BA}(\bm{r})\equiv\mathcal{T}_{B}(\bm{r})-\mathcal{T}_{A}(\bm{r})$ or $\delta\tau_{BA}(\bm{r})\equiv\delta\tau_{B}(\bm{r})-\delta\tau_{A}(\bm{r})$ explicitly as a difference to emphasise that the time delay difference involves two separate measurements which are both functions of position. Yonehara (1999); Yonehara et al. (2003) proposed and Goicoechea (2002) investigated constraining the locations of spatially separated flares within gravitationally lensed quasars via measurements of this time delay difference.

The time delay difference between corresponding photometric signals in images A and B is solely determined by cosmology through the lensing distance, the geometry of the lensing configuration and the spatial separation within the source in the plane of the sky. This unusual property inspired the key proposal of paper I: to extract cosmological information via the ratio of angular diameter distances

whilst alleviating the systematic uncertainties associated with lens modelling in using time delay measurements. The lens redshift and image positions can be measured directly. The time delay within an image $\mathcal{T}_{X}$ is found by subtracting the time $t_{X}$ we see a particular signal in the continuum flux, from the time $t^{\prime}_{X}$ at which we see the corresponding signal in the flux of a broad emission line originating from a responding BLR cloud at $\bm{r}$.

The technique as presented, centred on the usage of reverberation mapping to determine the position of the BLR cloud $\bm{r}$ responding at a given measured time delay in each image $A$ and $B$, is ill-motivated rather than one with limitations which could be directly mitigated by future improvements in observations. We summarise the interlinked difficulties as follows

1. 1.

   Reverberation mapping only gives general constraints on the geometric structure of the BLR. It does not uniquely map the location of individual points in the BLR, meaning that Equation (9) is underdetermined.

2. 2.

   The projection of the BLR cloud position on the image-image axis can not be determined from reverberation mapping, but requires e.g. spatially resolving each image.

3. 3.

   In the general case $\tau_{\textsc{rm}}\gg|\delta\tau_{AB}|$; the flux from each image is largely insensitive to differential lensing.

4. 4.

   The determination of the geometric parameters is not a priori independent of the source position.

A main challenge of exploiting Equation (9) lies in the task of pinpointing the position of the BLR cloud responding at a given measured time delay in each image. This position also enters separately in the dot product in the denominator of the same equation, i.e. requiring the projection of the BLR cloud position on the image $A$ - image $B$ axis on the sky, which assuming a spatially unresolved BLR is not directly measurable. For example, Goicoechea (2002) determined the position of a flare only along the image-image direction for a given image pair.

A key assumption of Yonehara et al. (2003) was that these flares are discrete, spatially separated and not physically correlated. If there is no physical correlation between discrete flares, then the shape of the flux arising from each of the flares is unique. Together with either lack of or minimal superposition of the flux, this allows each flare to be identified between images.

This assumption is violated when considering reverberation mapping of the BLR, which does not pinpoint the location from discrete flares. Rather, reverberation mapping utilises the response of an extended region to the same continuum signal to give general constraints on the BLR structure and kinematics in terms of geometric parameters. Determining the time at which we see the response from a particular position in a given image, as required by Equation (9), is in general an underdetermined problem as the response is a superposition of spectral flux from over the entire BLR.

The second consequence of the assumption that flares are not physically correlated is that the time between the discrete flares, denoted by $\Delta t_{\mathrm{burst}}$ in Yonehara (1999), in general can be arbitrarily small whilst their spatial separation is arbitrarily large. This allows the fulfilment of the condition $\Delta t_{\mathrm{burst}}\lesssim|\delta\tau_{AB}|$ required for the time delay difference to have an appreciable effect. Choosing to identify $\Delta t_{\mathrm{burst}}$ as a characteristic difference in $\tau_{\textsc{rm}}$ between points on the BLR, this condition was fulfilled in Paper I in the special case of a face-on thin ring BLR geometry which responds simultaneously, but no longer holds for a general configuration in which the difference in $\tau_{\textsc{rm}}$ is on the scale of the light crossing time over the BLR structure. In the general case $\tau_{\textsc{rm}}\gg|\delta\tau_{AB}|$ and so the dependence of the flux from each image on differential lensing is very small.

To see this clearly, in the time delay within the image $\mathcal{T}_{X}$ given by Equation (7) the terms $\frac{D_{d}}{D_{ds}}$, $(1+z_{s})$ and $(1+z_{d})$ are all on the order of unity, whereas $\alpha_{X}$ is on the order of at a minimum $10^{-6}$ to a maximum of $10^{-3}$ radians considering galaxy scale to cluster scale lensing. That is, in general the time scale for the differential lensing (minutes to days, depending on the lens mass) is very small compared to the time scale for reverberation mapping (tens to hundreds of days): $\delta\tau_{X}\ll\tau_{\textsc{rm}}$ implying $\mathcal{T}_{X}\sim\tau_{\textsc{rm}}$.

Finally, the determination of the geometric parameters using reverberation mapping is also not a priori independent of the distance ratio $\frac{D_{d}}{D_{ds}}$ nor the source position $\bm{\beta}$. We show in Section 5 of the present paper however that geometric parameters may be inferred accurately without either further data (e.g. image of the Einstein ring) or lens modelling.

As an additional note, Tie & Kochanek (2018) showed in the context of delayed emission across an accretion disk, whilst the differential lensing time delay given by (6) is negligible, the terms equivalent to the reverberation mapping terms in this work no longer cancel between images when considering the effect of microlensing. Rather, differential magnification over the source is included as an extra weighting factor to these terms which is image and time dependent (as the source moves relative to the stars causing the microlensing). Microlensing thereby adds time delays on the scale of the light crossing time of the relevant quasar structure to the standard time delay (2), potentially causing systematic problems for time delay cosmography; although no evidence for this effect has been observed from current data (e.g. Wong et al., 2020). We do not explore this microlensing effect as it does not impact the above conclusions centring on the potential-free expression for the differential lensing time delay.

Reverberation mapping provides constraints on the overall geometry of the BLR since measurements of $T_{X}$ provides information on $\bm{r}$; and in addition, the observed wavelength within the broad emission line is proportional (assuming relativistic effects are negligible¹¹1A full analysis would consider the impact of relativistic effects which are significant for small BLR radii, such as gravitational redshifting, the relativistic Doppler effect and relativistic beaming, e.g. Corbin (1997); Goad et al. (2012). A consequence of the relativistic Doppler effect is that in general the observed wavelength is not simply proportional to line-of-sight velocity.) to the line-of-sight velocity of the BLR cloud. If the velocity of a BLR cloud may be related to its position via a model for the velocity field, then the line-of-sight velocity gives additional information on the position of the BLR cloud as well as the kinematics. A “cloud” refers to small individual entities in the emission line region in reverberation mapping literature; we interchangeably refer to either a cloud or a point particle of the BLR distribution.

We first review the unlensed case (Blandford & McKee, 1982; Peterson, 1993; Netzer, 1990, 2013). Let $\bm{w}$ be the velocity coordinates of a BLR cloud, whereas $v_{z}$ denotes the line-of-sight velocity variable. The broad emission line spectral flux responds to the continuum flux $C(t-\tau_{\textsc{rm}})$, as given in generality by

where the line-of-sight projected 1D velocity distribution $g(\bm{r},v_{z})$ is defined in terms of the full velocity distribution $f(\bm{r},\bm{w})$ as

The responsivity term $j(\bm{r})$ is in general dependent on the position (or if isotropic, simply the radius) of the cloud and contains the physics of the photoionisation of the BLR. The continuum flux received by a gas cloud at time $t^{\prime}$ and position $\bm{r}$, emitted by the central source at an earlier time $t^{\prime}-\tfrac{r}{c}$, is $j(\bm{r})C(t^{\prime}-\tfrac{r}{c})=\frac{\varepsilon(\bm{r})}{4\pi r^{2}}C(t^{\prime}-\tfrac{r}{c})$ and plays the role of the emissivity of the cloud where $\varepsilon(\bm{r})$ is the reprocessing coefficient of the cloud.

The general reverberation mapping problem involves finding the deconvolution of

since all of the geometrical information is contained within the transfer function $\Psi(v_{z},t)$. When the transfer function is visualised using a heat map in $(v_{z},t)$-space, it is commonly referred to as a “velocity-delay map”. This deconvolution problem may be approached using maximum entropy fitting techniques, or regularised linear inversion; and then the outputted velocity-delay map may be compared qualitatively to the theoretical velocity-delay maps produced by simple models. Recovering the transfer function and thereby inferring the physical and kinematic distribution of the BLR using reverberation mapping has been, up until recent improvements of the quality of reverberation data sets, regarded as a technique in a developmental state (Shen et al., 2014; Mangham et al., 2019; Cackett et al., 2021). Rather than solving this ill-posed inverse problem, reverberation mapping has been often limited to finding only the mean radius for the BLR through measuring the mean time delay for emission lines using cross-correlation analyses, from which it is possible to constrain the mass of the central black hole. Forward modelling of the BLR also offers an alternative Bayesian approach to directly (without explicitly finding the transfer function) providing best-fit parameters of a flexible BLR geometry (Pancoast et al., 2011).

The transfer function, or the Green’s function for the system, is the line response to a Dirac-delta continuum pulse; replacing $C(t^{\prime}-\tfrac{r}{c})$ with $\delta(t^{\prime}-\tfrac{r}{c})$ gives

The second equality holds when there exists a velocity model $\bm{u}(\bm{r})$ whose corresponding line-of-sight velocity model is $u_{z}(\bm{r})$, such that we may write the velocity distribution as $f(\bm{r},\bm{w})=n(\bm{r})f(\bm{w}|\bm{r})=n(\bm{r})\delta(\bm{w}-\bm{u}(\bm{r}))$ where $n(\bm{r})$ is the number density of responding clouds. This gives the line-of-sight projected velocity distribution as $g(\bm{r},v_{z})=n(\bm{r})\delta(v_{z}-u_{z}(\bm{r}))$. We can interpret the general problem as one of finding a probability density function under a change of variables from positions $\bm{r}$ to $(v_{z},T_{X})$. As the dimensionality of the original set of random variables is not the same as the new set of random variables, the probability density functions are written using Dirac delta generalised functions.

Including the effect of differential lensing on the transfer function simply involves substituting the reverberation mapping time delay function $\tau_{\textsc{rm}}(\bm{r})$ with the time delay function $\mathcal{T}_{X}(\bm{r})$ and denoting the time variable with the subscripted $T_{X}$ to emphasise that it is measured with respect to the lensed continuum signal in each image $X$. It follows that the transfer function and the observed broad emission line spectral flux of each image are respectively given by

and

Although overall deductions about the BLR geometry and kinematics made be made from finding the transfer function or else by direct forward modelling, determining the time at which we see the response from a particular position in a given image, as is required by Equation (9), is in general an underdetermined problem. For the simple BLR models usually considered, e.g. where the (emissivity-weighted) number density distribution is modelled as uniform over the line, surface or volume to which the particles are confined and does not add an extra parameter, the full transfer function does not give information further than its support. For inhomogeneous BLR models which feature over- or under-densities of clouds in particular spatial regions however, the form of the transfer function may further constrain the number density distribution over that region. It can still be somewhat useful to find the theoretical transfer functions for idealised models for the visualisation of the method.

We therefore present the form of the transfer function in the case of a uniform thin inclined co-planar ring in Section 4.2 below, and leave the details as well as calculations for a disk and spherical shell to Appendix A. Throughout these calculations, we assume $j(\bm{r})$ to be a constant (e.g. Pancoast et al., 2011), i.e. that the radiation received and re-emitted by all BLR clouds is constant. A more physically accurate assumption of the dependence of $j(\bm{r})$ on, for example, radius does not impact the overall findings discussed in Section 3 of this paper. As light is simply redistributed in time due to the differential lensing, the overall effect on the spectral flux tends to be a small change in magnitude of the transfer function whilst spanning a correspondingly larger or smaller domain in time compared to reverberation mapping without lensing.

Since the available data comes from observations of each unresolved lensed image, the most probable BLR geometric parameters need to be estimated together with lensing parameters. Consequently, the estimates of BLR parameters which are to then be used for estimation of the distance ratio $\frac{D_{d}}{D_{ds}}$ according to the proposed method are, through the physical deflection angles, implicitly dependent on the source position $\bm{\beta}$. However, we find that it is possible to recover the desired BLR geometric parameters even if the lensing parameters (i.e. the deflection angles) are not constrained, as may be expected from the prior discussions on the smallness of the effect of the lensing for all physically reasonable scales.

As an illustration of this component of the method, we demonstrate finding these parameters using Bayesian parameter estimation, assuming that the transfer function has already been recovered, of the least degenerate, idealistic scenario of contriving the BLR model as restricted to a thin inclined ring. The parameters for this model which we wish to input into the expression for the distance ratio, Equation (10), are the black hole mass $M$, the inclination angle $i$ and BLR radius $R$.

In Bayesian data analysis, the posterior distribution $\mathcal{P}(\bm{\Theta}_{M}|\bm{\mathcal{D}},\mathcal{I})$ is the probability distribution of a given set of parameters $\bm{\Theta}_{M}$ belonging to a model $M$ given some data $\bm{\mathcal{D}}$ and prior information $\mathcal{I}$. The posterior distribution may be determined through Bayes Rule:

where $\mathcal{P}(\bm{\mathcal{D}}|\bm{\Theta}_{M},\mathcal{I})$ is the likelihood, $\mathcal{P}(\bm{\Theta}_{M}|\mathcal{I})$ is the prior, and the evidence or marginal likelihood is $\mathcal{P}(\bm{\mathcal{D}}|\mathcal{I})=\int_{\Omega_{\bm{\Theta}_{M}}}\mathcal{P}(\bm{\mathcal{D}},\bm{\Theta}_{M}|\mathcal{I})\mathcal{P}(\bm{\Theta}_{M}|\mathcal{I})d\bm{\Theta}_{M}$. The integral is over the domain of the parameter space $\Omega_{\bm{\Theta}_{M}}$ and as we are only considering a single model, we drop the subscript $M$ for the rest of this work. The posterior may be estimated via numerical methods, usually involving an algorithm that approximates the posterior by generating and weighting a set of discrete samples. Since evaluating the multi-dimensional integral to compute the evidence is a key difficulty, common numerical methods such as Markov Chain Monte Carlo (MCMC) rely on drawing these samples from a distribution proportional to the posterior distribution. In this work we use the emcee (Foreman-Mackey et al., 2013) implementation of the MCMC algorithm. We assume non-informative priors, such that the posterior is directly proportional to the likelihood

and the maximum a posteriori estimate of $\bm{\Theta}$ is equivalent to finding the maximum likelihood estimate.

We define our likelihood function, first specifying our model. Let us assume that the original continuum function $C_{\text{original}}(t_{X}^{\prime})=\delta(t_{X}^{\prime}-t_{i})$ is a Dirac Delta function which is then convolved or blurred in the time domain by a Gaussian function such that:

where $\epsilon$ is the width or standard deviation of the blurred continuum, which allows for numerical sampling of the posterior. Since for a thin inclined ring, the transfer function $\Psi_{X}(v_{z},T_{X})$ is given by a sum of Dirac-delta signals (26), the line response with $t_{i}=0$ becomes

where

and

and so, fixing $\epsilon$, the parameters associated with this model are $\bm{\Theta}=\left(\log_{10}\left(\tfrac{M}{M_{\odot}}\right),i,\log_{10}\left(\tfrac{R}{\textrm{m}}\right),\{\hat{\alpha}_{X,x},\hat{\alpha}_{X,y}\}\right)$ where $\{\hat{\alpha}_{X,x},\hat{\alpha}_{X,y}\}$ denotes the set of all the physical deflection angles over all the available images (through which there is an implicit dependence on $\frac{D_{d}}{D_{ds}}$ and $\bm{\beta}$). For the purpose of parameter estimation, the model is considered a function of the parameters $\bm{\Theta}$ instead of the variables $(v_{z},t_{X}^{\prime})$, which are evaluated at values fixed by each data point and then binned: $L_{x}(v_{z},t_{X}^{\prime})\to L_{X,jk}(\bm{\Theta})$, where $j$ and $k$ respectively denote the time and line-of-sight velocity bins of the data.

Rather than work with the analytic expression (38) directly, which due to the ${V^{+}(R,v_{z})}^{-1}$ term is singular when $v_{z}^{2}=u^{\prime 2}(r)$, we instead obtain $L_{X,jk}(\bm{\Theta})$ numerically, utilising $L_{x}(v_{z},t_{X}^{\prime})\propto\frac{\partial^{2}N}{\partial u_{z}\partial\mathcal{T}_{X}}(v_{z},t_{X}^{\prime})$. For each BLR cloud we may draw a position sample from the thin ring number density distribution function. To each sample we then assign a time and a velocity from the time delay and velocity models (24a) and (24b), before performing a 2D Gaussian blur. The resultant time delay and line-of-sight velocity values were binned with bin widths of roughly 1d and 50 km s^-1 respectively.

We generate mock data by fixing the model parameters at a chosen set of true values $\bm{\Theta}_{\text{true}}$ and adding measurement noise. As photon statistics follow a Poisson distribution, the scaled Poisson distribution corresponding to the flux data has a mean at the model value $L_{X,jk}(\bm{\Theta_{\text{true}}})$ and a variance given by $\sigma_{jk}^{2}=aL_{X,jk}(\bm{\Theta_{\text{true}}})+\sigma^{2}_{\textrm{abs}}$, where $a$ is a constant and there is an additional background contribution $\sigma^{2}_{\textrm{abs}}$ to the noise. The parameters for the noise are chosen to correspond to a peak signal-to-noise ratio of approximately 10 with $a\approx 0.05\,\text{max}(L_{X,jk}(\bm{\Theta_{\text{true}}}))$. As the photon count is large, the distribution is approximated as Gaussian

and the probability density of an individual datum $\mathcal{D}_{X,jk}$ is

The likelihood function for the line response from a single image $X$, i.e. the joint probability density function $\mathcal{P}(\bm{\mathcal{D}_{X}}|\bm{\Theta},\mathcal{I})$, is given (assuming independent data) by the product of the probability densities of the individual measurements $\mathcal{D}_{X,jk}$:

The log-likelihood for the line response for a single image $X$ is thus

To find the posterior taking into account the data from all images (assuming independence), we can either iterate over the data from each subsequent image; or perform the analysis in one step, where the full likelihood function is the product of the individual ones

where $\{\bm{\mathcal{D}}_{X}\}$ is the set of data from all images.

We ran emcee until the Markov chain was longer than 100 times the integrated autocorrelation time estimate for the chain for each parameter, with the estimated autocorrelation time changing by less than 1%. Due to the highly restrictive assumptions of this BLR geometric model, the geometric parameter values are very well determined; whereas the deflection angles are unconstrained with the posterior distributions equivalent to the uniform priors. This result for realistic lensing scales could be surmised from the data by eye; although even if the scale of lensing is increased to unrealistically large magnitudes, the geometric parameter values remain very well determined for this model. The corner plot of the one and two dimensional projections of the posterior probability distributions in the case of galaxy-scale lensing is shown in Figure 6.

In summary, (Ng & Lewis, 2020) proposed a method of determining $\frac{D_{d}}{D_{ds}}$ without lens modelling, via reverberation mapping and measuring times delays of differentially lensed, unresolved quasars. This method suffers from the following:

1. 1.

   Reverberation mapping does not uniquely map the location from discrete flares as is required; rather it only gives general constraints on the geometric structure of the BLR. This is a problem of underdetermination.

2. 2.

   The projection of the BLR cloud position on the image-image axis is needed, which requires further data (i.e. spatially resolving each image).

3. 3.

   The time scale for differential lensing is on the order of minutes to days compared to tens to hundreds of days for the reverberation mapping time delay; $\tau_{\textsc{rm}}\gg|\delta\tau_{AB}|$ for a general BLR geometry.

4. 4.

   Finally, the determination of the geometric parameters is not a priori independent of $\bm{\beta}$. However, the effect of differential lensing is sufficiently small ($\tau_{\textsc{rm}}\gg|\delta\tau_{AB}|$), implying the geometric parameters may be inferred accurately without additional data or lens modelling.