Determination of Lens Mass Density Profile from Strongly-Lensed Gravitational-Wave Signals

As the interferometers detecting gravitational waves are upgraded, improving their sensitivity, the probability of observing strong lensing increases. Once a detection is made, it will be critical to gain as much information as possible about the lensing object from these observations. In this work, we present a methodology to rapidly perform model selection between differing mass density profiles for strongly lensed gravitational wave signals, using the results of the fast strong lensing analysis pipeline GOLUM. We demonstrate the validity of this methodology using some illustrative examples adopting the idealised singular isothermal sphere and point mass lens models. We take several simulated lensed signals, analyse them with GOLUM and subject them to our methodology to recover both the model and its parameters. To demonstrate the methodology's stability, we show how the result varies with the number of samples used for a subset of these injections. In addition to the analysis of simulations, we also apply our methodology to the gravitational wave event pair GW191230--LGW200104, two events with similar frequency evolutions and sky locations, which was analysed in detail as a potential lensing candidate but ultimately discarded when considering the full population and the uncertain nature of the second event. We find a preference for the singular isothermal sphere model over the point mass, though our posteriors are much wider than for the lensed injections, in line with the expectations for a non-lensed event. The methodology developed in this work is made available as part of the Gravelamps package of software.


INTRODUCTION
In the wake of a rapidly increasing number of gravitational wave (GW) detections, with 90 confirmed detections as of the latest catalog of events, GWTC-3 (The LIGO Scientific Collaboration et al. 2021), and with the coming years set to deliver significant improvements in the sensitivity of the existing network of ground-based GW detectors, currently comprising Advanced LIGO (Aasi et al. 2015), Advanced Virgo (Acernese et al. 2015), and KAGRA (Akutsu et al. 2021;Aso et al. 2013;Somiya 2012), the opportunities are manifest for the nascent field of GW astronomy to pursue a wide variety of scientific investigations.
One example of a particularly active area of GW research is that of GW lensing (Hannuksela et al. 2019;Abbott et al. 2021;The LIGO Scientific Collaboration et al. 2023;Janquart et al. 2023b).Analogously to the deflection of light by a strong gravitational field (Einstein 1936)-which formed one of the first major tests of Einstein's general relativity (GR) (Dyson et al. 1920)-GW signals may be lensed by intervening objects such as stars, black holes, galaxies and galaxy clusters (Oha-nian 1974;Wang et al. 1996;Takahashi & Nakamura 2003;Smith et al. 2017;Mishra et al. 2021).
As of the end of third observing run of the current ground-based detector network, searches for lensed GW signals have not yet yielded a confirmed detection (Abbott et al. 2021;The LIGO Scientific Collaboration et al. 2023;Janquart et al. 2023b).However, with the increased detection rate of GW events expected during the next observing run, there is a signficant probability that one or more of these detections may be identified as strongly lensed (Smith et al. 2017(Smith et al. , 2018;;Ng et al. 2018;Li et al. 2018;Wierda et al. 2021;Xu et al. 2022;Smith et al. 2023).In this scenario the GW signal would be lensed by an intervening galaxy (Dai et al. 2017;Ng et al. 2018;Robertson et al. 2020) or galaxy cluster (Smith et al. 2017(Smith et al. , 2018;;Robertson et al. 2020), forming multiple distinct "images".These images would be potentially magnified, time shifted, and overall phase shifted copies of the original signal that may be detected as repeated events (Dai & Venumadhav 2017;Janquart et al. 2021;Liu et al. 2021;Lo & Magana Hernandez 2021;Janquart et al. 2023a).
Detection of such a lensed event would not only be of enormous scientific interest in itself, but would also immediately leverage additional scientfic benefit from the information carried by lensed GW signals from compact binary coalescence events-much as the detection of GW signals from such events in the first place has opened the door to completely new ways in which to investigate these phenomena.Examples of the kinds of investigations that could be enabled by lensed GW signals include, but are not limited to, improving the localisation of detected GW events (Hannuksela et al. 2020;Wempe et al. 2022), precision cosmography (Sereno et al. 2011;Liao et al. 2017), or tests of the speed of gravity (Fan et al. 2017), or general relativity itself (Ezquiaga & Zumalacárregui 2020;Goyal et al. 2021).
Moreover, investigations of lensed GW systems may also reveal import information about the nature of the lens.For example, it has been demonstrated that the mass density profile of the lens will influence the properties of lensed GW signals (Takahashi & Nakamura 2003;Cao et al. 2014), so that observation of the latter may allow the former to be constrained.Investigations have therefore been carried out into how effectively GW observations can allow the characterisation of the signal and thus constraint of the lens model parameters in both the microlensing (Mishra et al. 2021)-where lenses would be inidividual compact objects such as stars or black holes-and strong lensing regimes (Herrera-Martín et al. 2019;Hannuksela et al. 2020;Wempe et al. 2022;Tambalo et al. 2022) regimes for specific models of lens.Previous work has also taken the additional step of performing model selection for microlensed GW signals (Wright & Hendry 2022).However, it is important to note that within that work, the presented methodology was applicable to only signle images and this proved to be insufficient to distinguish between lensing models when observing the individual images of strongly-lensed GW signals.
In this work, we present a new methodology for rapidly performing model selection on the results of model-agnostic joint parameter estimation analyses of strongly lensed multiplets.We demonstrate the effectiveness of this methodology by using it to analyse 65 simulated pairs of lensing signals and, in addition, investigate the pair GW191230 and LGW200104 which were considered by the LVK lensing search to be the pair of events with the highest likelihood of lensing-albeit we stress that the non-lensed hypothesis was still preferred and that LGW200104 is disfavoured as a real GW event by its extremely low p astro (Janquart et al. 2023b).
The work is structured as follows.In Section 2, we introduce the basics of strong lensing including the the-oretical basis, some potential models of lens, as well as the means by which strong lensing may be searched for model-agnostically in GW data.Section 3 then outlines the methodology by which we may perform model selection on the results of this model agnostic analysis.Section 4 demonstrates the testing that has been done to verify the validity of the methodology both on the set of injections and on the trigger data outlined above.Finally, we present the conclusions and future outlook from this work in Section 5.

STRONG LENSING OF GRAVITATIONAL WAVES
In general, lensing of GWs is an amplification process.A standard, non-lensed GW signal-which is described by the phase amplitude ϕ obs (w, η) where w is the dimensionless frequency and η is the displacement of the source from the optical axis-experiences an amplification to become the observed lensed signal, ϕ L obs (w, η).The ratio between these wave amplitudes is termed the amplification factor (F ) and is such that the relation between the unlensed (h(f )) and lensed (h L (f )) GW strains is: (1) This amplification factor varies depending upon the mass density profile of the lensing object.However, it can be calculated for any profile using a single general expression (Schneider et al. 1992;Takahashi & Nakamura 2003): where y is a dimensionless form of the displacement from the source (η as described above) and the function T (x, y) yields the dimensionless time delay.
In the case of strong lensing, the mass of the lensing object is very high (which yields w ≫ 1) and the geometric optics approximation is valid.This results in only the stationary points of the time delay function contributing to the integral, reducing it to a summation over these points (Nakamura & Deguchi 1999).These stationary points correspond to the individual images produced in the strong lensing process and the amplification of the j th such image is given by (Takahashi & Nakamura 2003;Dai & Venumadhav 2017;Ezquiaga et al. 2021) As can be seen, this yields three separate observable effects on the lensed signals: the magnification (µ j ) of the signal's amplitude, the time delay (t j ) between the signals due to the difference in path length between the images, and the Morse phase (n j ).This final property is an overall phase shift of the waveform that may take one of three distinct values, n = 0, 1/2, 1 depending on whether the image corresponds to a minimum, saddle point, or maximum of the time delay function respectively.

Lensing Models
As has been mentioned above, the general functional form of the amplification factor is given by Eq. ( 2), but its final expression will vary for each different model for the mass density profile of the lens.There are many such models that have been considered to describe possible lenses.Here, we briefly describe the two models that have been used in this work as examples, to demonstrate the validity of the methodology.These models were chosen due to their ability to be solved analytically which allows for direct testing of the method as well as the fact that they produce two and only two images which again allows a complete treatment of these two models as testing cases.In follow-up work we will expand our analysis to consider more general mass density profiles that provide large multiplet image sets, so that the effect of detecting a subset of images may be taken into consideration.

Point Mass
The simplest model of the lensing object is that of a point mass.Such a lens produces two images, regardless of the lens-source configuration.One image will be a minimum of the time delay function, and one will be at a saddle point.Consequently, following on from Eq. (3), the amplification factor is (Takahashi & Nakamura 2003) The image magnifications for the point mass case are given in terms of the dimensionless source position, y, by µ ± = 1/2 ± (y 2 + 4)/(2y y 2 + 4) and the time delay between the two images is given by ∆T = y( y 4 + 2)/2 + ln ( y 2 + 4 + y)/( y 2 + 4 − y) .

Singular Isothermal Sphere
As a step up in complexity from the point mass profile, the Singular Isothermal Sphere (SIS) is widely used to describe the dark matter halos of galaxies, due to its ability to describe the flat rotation curves observed for these systems (Binney & Tremaine 1987).The SIS profile models the galaxy as an extended luminous matter object embedded within a larger dark matter halo.
However, the weakness of this profile is that it suffers from a non-physical central singularity.
For an SIS lens, the configuration of source and lens may alter the number of images produced, with in one case a single image produced at higher source positions, and in the other case two produced at lower.Consequently, the amplification factor-again following Eq.( 3)-is piecewise and given by (Takahashi & Nakamura 2003) where in this case the magnifications are given by µ ± = ±1 + 1/y and the time delay is given by ∆T = 2y.

Strong Lensing Identification
There are a number of means by which one may identify a strong lensing multiplet.These range from examining the overlap between the GW source posteriors derived from the individual images, as outlined in Haris et al. (2018), which is speedy but does not provide confirmation of lensing status-and more importantly in our context does not provide detailed information on the observed lensing parameters-to full joint parameter estimation of the two image candidates using the method described in Liu et al. (2021); Lo & Magana Hernandez (2021); Janquart et al. (2023a), which is complete but computationally expensive.In this work, we focus on the methodology presented in Janquart et al. (2021Janquart et al. ( , 2023a) ) implemented in the python package GOLUM (Janquart et al. 2022a) as a middle ground in which events are provided with sufficiently accurate estimates on the observable lensing parameters, as well as providing confirmation of lensing status.Whilst we refer the reader to Janquart et al. (2021Janquart et al. ( , 2023a) for a complete explanation of how this identification is performed, we briefly outline the methodology here-for simplicity considering the two-image case only, although the method is valid for any number of images.
A pair of lensed images, h j L (t j ; θ, Λ j ), are described in terms of the binary parameters, θ and the individual image parameters, Λ j for the j th image.When examined under the lensing hypothesis, H L , under which the binary parameters should be identical, the joint evidence neglecting selection effects is given by (Liu et al. 2021;Lo & Magana Hernandez 2021) where the term p (θ, Λ 1 , Λ 2 ) is the prior, and the terms p (d j |θ, Λ j ) are the individual likelihoods (Veitch & Vecchio 2010).This may be compared with the joint evidence in the unlensed hypothesis, H U , which is the product of the individual likelihoods: The two are compared in the "coherence ratio": Calculation of this joint evidence and coherence ratio in full is the basis for complete joint parameter estimation such as that outlined in Liu et al. (2021); Lo & Magana Hernandez (2021).However, the process of calculating the coherence ratio may be sped up by instead considering the conditional evidence alongside the individual evidences.This allows acceleration of computation by means of importance sampling and a look-up table-the full details of which are described in Janquart et al. (2021Janquart et al. ( , 2023a)).To outline how this is done, the joint evidence is rewritten in terms of the conditional evidence as The conditional evidence p (d 2 |d 1 , H L ) may be evaluated as a "marginalised" likelihood of the form where the binary parameters, θ, have been replaced with the effective parameters, Θ which absorb the lensing magnification into the observed luminosity distance (D obs L = D L / √ µ, where D L is the source luminosity distance) and the time delay into the observed coalescence time (t obs c = t c +t, where t c is the coalescence time without lensing and t is the time delay), and the individual image parameters Λ j have been replaced with the relative image parameters Φ, i.e. the parameters relative to the first of the images.
The first term of the integrand is simply the likelihood of the second event averaged over the posterior samples of the first event.Whilst this term would already be faster to compute than the full joint evidence from Eq. ( 6), it is further sped up as the reuse of the posterior samples from the first event allows the construction of a lookup table allowing speedy computation of Eq. ( 9).
This culminates in the full calculation of the coherence ratio as We note that this coherence ratio is not a full Bayes factor as it does not include the selection effects outlined in Lo & Magana Hernandez (2021).Nevertheless, their inclusion can be done straightforwardly as a postprocessing step to the GOLUM analysis.

MODEL SELECTION OF LENSED GRAVITATIONAL WAVE SIGNALS
Strong lensing identification methodologies as outlined above have created a framework to answer the question of: "for a given pair of images, are these images lensed?" but have generally not answered the question of "by what?".In large part this is because the search methodologies are model-agnostic as to not misidentify an event pair because the method assumes an incorrect model.Some initial steps towards answering the question of the lensing object have been made, however, such as in Janquart et al. (2022b) where the authorial suggestion was to reweight the detection statistic using model information from various catalogs built from differing models such as those found in Haris et al. (2018), More & More (2022), or Wierda et al. (2021).Similarly, the methodology presented in Lo & Magana Hernandez (2021) describes the inclusion of the model for the calculation of selection effects to achieve a more complete Bayes factor in the context of "lensed vs unlensed".Whilst the first method can deal with more realistic models, it requires building an extended catalog when one wants to explore a realistic lens model, assuming particular source and lens populations and is therefore more prone to systematics than the direct application of a lens model.On the other hand, for the second method, it requires an analytical expression for the magnification probability distribution, which reduces the model it can consider.We note that the latter method could also make use of the results coming from a catalog as a model for the magnification probability distribution but it would then suffer the same caveats as the method presented in Janquart et al. (2022b).In both cases, the efforts to use these methods on candidate lensed signals is detailed in Janquart et al. (2023b).
An additional complication for more detailed analysis of strong lensing signals is the need for the detection and identification of multiple signals.For instance, analysis of individual images using a specific lens model as outlined in Wright & Hendry (2022) in the geometric optics case does not yield successful differentation between lens models due to the complete degeneracies between the effects of strong lensing and other binary source parameters when only a single image is considered.
A means of performing model selection based on the output of these detection pipelines, including the multiple images and that does not require resampling or the construction of extended catalogs which requires assumptions on population models, would be beneficial to future lensing searches, and it is such a methodology that we present here.
The goal of this methodology is to find the evidence for a given lens model, here termed H mod , which consists of a set of lens parameters, Ψ.By direct application of Bayes' theorem, the evidence for that model is given by: To simplify the following, we here define the evidence for a chosen model as Z = p (d 1 , d 2 |H mod ) and the model likelihood as L (H mod ) = p (d 1 , d 2 |Θ, Φ, Ψ, H mod ).We consider calculation of the inverse evidence which, as will be made clear, is easier to solve and once solved may be trivially inverted back to the evidence.This is given from the above as The posterior forming the numerator of Eq. ( 12) may be expanded further using the chain rule to yield The lattermost term of Eq. ( 13) is in fact, insensitive to the lens model since the apparent lensing parameters and the effective BBH parameters may be fully determined from the data.As such, p (Θ, Φ|d 1 , d 2 , H mod ) ≡ p (Θ, Φ|d 1 , d 2 ) and thus, Eq. ( 12) becomes This may be solved by sampling p (Θ, Φ|d 1 , d 2 ), computing the remaining terms for that sample and averaging the ratio over all samples, i.e.
The quantity over which we sample, p (Θ, Φ|d 1 , d 2 ), is the output of a model-independent joint parameter estimation pipeline (whether using a full joint parameter estimation or a GOLUM style approach).In addition, the prior and likelihood are known and easily calculable, with the joint likelihood being calculated using the implementation from Janquart et al. (2023a) with the produced lensed waveform.Finally, the numerator-the conditional probability of the lens parameters given the model-may be computed from the relations between the observables parameters and the model parameters as outlined in Sections 2.1.1 and 2.1.2.
As a result, Eq. ( 15) offers a means by which to compute the evidence for a given lens model directly from the output of the pipelines that would be used to claim a lensing detection.It also bypasses the need for sampling an extended parameter space meaning that the computation would be extremely rapid upon the completion of the model-agnostic search.
We implement Eq. ( 15) into the Gravelamps package (Wright et al. 2022) for the models outlined above in which Ψ reduces to the redshifted lens mass, M Lz , and the dimensionless source position, y.Whilst these models are relatively simple, they allow the testing of the methodolgy in a controlled environment in which the observable to lens parameter conversion may be analytically calculated and therefore may be more closely monitored when testing the methodology.However, the implementation within Gravelamps will allow for extending the method to more complex realistic models leveraging its already extant capabilities to rapidly compute lens amplification factors (Wright & Hendry 2022).

INVESTIGATION OF METHODOLOGICAL PERFORAMNCE
In order to investigate the performance of the methodology in real-world scenarios, it has been subjected to a number of tests.The primary source of investigation was the performance of the methodology on a constructed set of 75 lensed GW signal pairs with known parameters.This allowed investigation of the ability of the method to identify the lensing model as well as the recovery of the model-specific lensing parameters.The stability of the evidence calculation was also investigated by looking at the results of varying the number of samples considered.Finally, whilst a real lensing event remains unavailable, examination of the methodology on real-world data was performed by considering the event pair GW191230-LGW200104: the event pair that was determined by the lensing searches in Janquart et al. (2023b) to be the highest significance of the ultimately discarded candidates.

Injection Set Investigation
The main test of methodological performance was to investigate the results of subjecting a series of simulated lensing events to the method.The testing set consisted of 75 pairs of lensed GW signals.To consider a realistic testing set, the mass, spin, and redshift priors for these were chosen to reflect the inferred population from the GWTC-3 data, i.e. in-line with Abbott et al. (2022) 1 .The other parameters were selected from the usual priors on the parameters, see Janquart et al. (2021) for example.The events were chosen to be lensed by SIS lenses with a uniform prior on the redshifted lens mass between 10 6 M ⊙ and 10 9 M ⊙ and a power-law prior on the source position with α = 1 over the range of source positions resulting in multiple images, avoiding the direct hit and caustic cases, i.e. y = 0 and y = 1 respectively.
Once the sample source and lens parameters were drawn, the lensed pair signals were constructed and injected into a detector network consisting of the two LIGO and the Virgo detectors with noise representative of the expectations for O4.These were then subject to joint parameter estimation using the GOLUM pipeline (Janquart et al. 2021(Janquart et al. , 2023a) ) operating with the Nessai (Williams et al. 2021) nested sampler and the following priors on the model-agnostic lensing observables: • Relative Magnification: Uniform between 0.01 and 50 • Time Delay: Uniform over the range of the injected time delay ± 0.2 seconds • Morse Index : Uniform on the discrete values of 0, 0.5, and 1 The signal injection and recovery were done using the IMRPhenomXPHM waveform (Pratten et al. 2021).The resulting samples were then subjected to the methodology and the evidences for each of the models compared to arrive at the resultant Bayes factor.
Figure 1 shows the distribution of the recovered Bayes factors comparing the SIS and point mass lensing models in the form of both the raw histogram of the data as well as the inferred distribution constructed using a kernel-density estimator (KDE).All of the events considered demonstrate a positive log Bayes factor, indicating a consistent preference for the true SIS model.The minimum log Bayes factor of the considered events was 1.44, the maximum was 7.98, and the average log Bayes factor for the true SIS model across all of the events was 3.37.This is sufficiently large to demonstrate a consistent identification and preference for the correct model in these injections.
In addition to the model identification, the method produces candidate posteriors for the lens parameters in 1 To be precise, we took the maximum likelihood parameters of the population models given in (LIGO Scientific Collaboration and Virgo Collaboration and KAGRA Collaboration 2023).each of the models.These two may be interrogated to assess the ability of the method to constrain and recover these model specific parameters.Figure 2 illustrates an example of the recovered posteriors for the parameters of a single event from the testing set for the SIS model, with the true value indicated.In that case, the parameters demonstrate relatively tight constraints with the true value inside of the constructed posterior.
As expected, the tightness and accuracy of the constraint of the model-specific lensing parameters are correlated with the tightness and accuracy of the constraint of the model-agnostic lensing observables.This is demonstrated in Figure 3 which contains a subset of the total runs.In this figure, it can be seen that cases with relatively tight constraint of the relative magnification yield relatively tight constrain of the source position.Similarly in the case of the fifth pairing in the figure, the failure to constrain the relative magnification accurately in the model agnostic search results in a failure to constrain the source position accurately.

Stability of Methodological Performance
Whilst the methodology is relatively inexpensive to perform over all the samples from a given GOLUM investigation, preliminary results may be sought more quickly by using fewer samples.Similarly, the number of samples will differ between investigation of one event and the next, therefore one metric that needs to be investigated for the method is the number of samples that will yield a stable result.Figure 4 shows the evolution of the Bayes factor calculation from a number of events as progressively more samples are included.This was done by taking the full sample set and calculating the Bayes factor from every n th sample scaled by the maximal value calculated to allow multiple events to be overlaid.As can be seen from the figure, the estimated value of the Bayes factor shows significant fluctuation when a small number of samples are used and an increasing narrowness as a greater number of samples are included, demonstrating that the solution is stabilising, with an approximate preliminary value being representative at the inclusion of every 100 th sample.
In addition to the overall fluctuation of the Bayes factor, a number of events demonstrated that at low sampling rates, the estimated value of the Bayes factor can appear to cluster around a value that is different from the value to which the estimate later appears to converge to.An example of this is shown in Figure 5. Similarly, the figure shows an approximate point at which a preliminary result-from the inclusion of a limited number of samples-is robustly representative of the true Bayes factor.Our results again would indicate that including every 100 th sample is sufficiently representative of the Bayes factor in this way.
In addition to the overall fluctuation of the Bayes factor, a number of events demonstrated that at low sampling rates, the estimated value of the Bayes factor would appear to cluster around a deviated value from the later settling value.An example of this is shown in Figure 5. Similarly, this would place an approximate point at which a prelimninary result would be represen-tative of the final value at the inclusion of every 100 th sample.

Example Deployment on GW191230-LGW200104
As discussed in the introduction, as of the end of the third observing run, the LVK lensing searches which include joint parameter estimation analyses that would be appropriate for this method have yet to yield any confirmed detections (Hannuksela et al. 2019;Abbott et al. 2021;The LIGO Scientific Collaboration et al. 2023).Consequently, there is no real lensing data on which to perform a test case analysis.To represent such a test case, we deploy the method on the event pair GW191230-LGW200104 which was identified to have the highest Bayes factor of the ultimately discarded candidates from searches thus far (The LIGO Scientific Collaboration et al. 2023;Janquart et al. 2023b).Similarly to Janquart et al. (2023b), we stress that we do not claim that this pair is a genuine lensing event.However, we will treat it as though it were for the purposes of this test deployment.Similarly, due to the relatively high SNR of the trigger, LGW200104 it is treated as a real GW event despite the very low p astro of 1% which would suggest this is unlikely in actuality (Janquart et al. 2023b).
The pair has been re-analysed using the GOLUM pipeline and the Nessai nested sampler in order to yield an equivalent model-agnostic search result to those of the injection set.Performing the model selection method on the results yielded a result that favours the SIS model with a log 10 Bayes factor of 3.65.This is consistent with the finding of Janquart et al. (2023b) that the consideration of the SIS model did improve preference for the pair as compared to the raw model-agnostic  .The resulting value of the Bayes factors found by using every n th sample of the full set, scaled by the maximal value from the calculations to allow multiple events to be shown.As can be seen, a preliminary estimate of the Bayes factor, from applying the methodology using every 100 th sample, would be robustly representative of the "final" estimate using all of the samples. .These are fairly broad as might be expected when including a trigger that is likely not a genuine lensing event or indeed a genuine GW event.Additionally, we note that there is distinct railing in the reconstructed source position posterior against the di-rect hit-i.e.y = 0-case which may too be an indication of the probably unlensed nature of the candidate pair.However, should this have been a genuine lensed event pair the median values of these reconstructed parameter posteriors for the SIS model would suggest a 5.08 +0.94 −2.09 ×10 10 M ⊙ lens at a source position of 0.23 +0.15 −0.14 .

CONCLUSION
With the increasing sensitivity of the interferometers, the detection of gravitational lensing of a GW event becomes increasingly likely.When such a detection occurs, the nature of the lens, as well as the source, will be studied.One of the fundamental questions about the lens will be about how the mass is distributed.We present here a method by which to rapidly determine the mass density profile using the output of the model-agnostic strong lensing joint parameter estimation pipelines without the need to build extended catalogs of lenses, under the assumption of source and lens populations.This method has been implemented within the python package Gravelamps (Wright et al. 2022) to augment its pre-existing capabilities to determine mass density profiles for microlensing events.
We have demonstrated the efficacy of the method using a pair of relatively simple lens models-the point mass and the SIS lenses-though it is not restricted to these models.In order to test the methodology's performance it was confronted with an injected set of 75 lensed BBH pairs following the observed distribution of events during the third observing run of the LVK detectors.In each case, the correct lensing model was identified with the lens model-specific parameter recovery falling in line with the recovery of the observable model-agnostic lensing parameters.
To test the stability of the results calculated, we observed the behaviour of the Bayes factors as more and more samples are included.Preliminary results from the methodology would be representative of the final solution using all of the samples from the inclusion of every 100 th sample.
Finally, whilst there remains no real lensing data on which to deploy the methodology, we test deployment on actual data by examining the GW191230-LGW200104 pair.It had the highest Bayes factor of the ultimately discarded pairs.Our analysis indicated that if the pair were genuinely lensed, which we stress is unlikely, then the preferred lens of the two models examined would be a 5.08 +0.94  −2.09 × 10 10 M ⊙ SIS lens at a dimensionless displacement from the optical axis of 0.23 +0.15  −0.14 .This is in line with the increased support for the pair when considering the SIS model in we note that the SIE model was the most preferred in that case.With the demonstration of viability for the method, future work to be undertaken would be to investigate more realistic models for the mass density profile of the lens, such as the aforementioned SIE, or the Navarro-Frenk-White model.For these more sophisticated models, additional work would need to be done to rapidly go from the lensing observables to the model-specific parameters without the relatively simple analytical relationships that exist for the point mass and SIS models considered in this work.

Figure 1 .
Figure 1.Distribution of the log Bayes factors comparing the SIS and point mass lensing models for the injection set.Shown in black is the raw histogram of these with the red overlay showing the inferred distribution from a kerneldensity estimator.All events analysed show a positive Bayes factor which indicates that the SIS was correctly preferred in all cases.

Figure 2 .
Figure 2. Constructed candidate posteriors on the redshifted lens mass (left) and source position (right).The raw histogram of the samples is presented in black, the inferred posterior distribution using a KDE in blue, and the vertical red dotted line illustrates the true source position value.In this instance, the posteriors indicate a tight constraint with the true values of the parameters inside of the posterior and close to its mode-a successful recovery.

Figure 3 .
Figure3.Comparisons between the posterior on the relative magnification as determined by the GOLUM pipeline and the constructed posterior on the source position for the SIS model using the method presented in this work.The true value of these parameters is indicated by the red vertical dotted line.As can be seen, the constraint on the model parameter is dependent upon both the width and accuracy of the lensing observable constraint.
Figure4.The resulting value of the Bayes factors found by using every n th sample of the full set, scaled by the maximal value from the calculations to allow multiple events to be shown.As can be seen, a preliminary estimate of the Bayes factor, from applying the methodology using every 100 th sample, would be robustly representative of the "final" estimate using all of the samples.

Figure 5 .
Figure5.The results of the stability analysis for an event which demonstrates an initial deviation from the later stabilisation value.Similarly to the main set, by approximately the inclusion of every 100 th sample, the preliminary result would be representative of the final result

Figure 6 .
Figure 6.Top: Posteriors on the relative lensing parameters for the GW191230-LGW200104 pair as determined by the GOLUM investigation.Bottom: The reconstructed posteriors for the source position and redshifted lens mass as calculated by the methodology presented in this work.

ACKNOWLEDGEMENTS1
The authors thank Chris Van Den Broeck for useful discussions and insights on related topics.MW acknowledges the support of the Science and Technologies Facilities Council (STFC) of the United Kingdom.JJ is supported by the research programme of the Netherlands Organisation for Scientific Research (NWO).MH acknowledges additional support from the Science and Technologies Facilities Council (Ref.ST/L000946/1).
based upon work supported by NSF's LIGO Laboratory which is a major facility fully funded by the National Science Foundation.
The authors are grateful for computational resources provided by the LIGO Laboratory and the Leonard E Parker Center for Gravitation, Cosmology, and Astrophysics at the University of Wisconsin-Milwaukee and supported by National Science Foundation Grants PHY-0757058 and PHY-0823459, and PHY-1626190 and PHY-1700765 respectively.