Gravelamps: Gravitational Wave Lensing Mass Profile Model Selection

We present the package Gravelamps which is designed to analyse lensed gravitational wave signals in order to constrain the mass density profile of the lensing object. Gravelamps does this via parameter estimation using the framework of bilby, which enables estimation of both the lens and the source parameters. The package can be used to study both microlensing and macrolensing cases -- where the lensing mass distribution is described by a point mass and extended mass density profile respectively -- and allows the user to easily and freely switch between the full wave optics and approximate geometric optics description. The performance of Gravelamps is demonstrated via simulated analysis of both mircolensing and macrolensing events, illustrating its capability for both parameter estimation and model selection in the wave optics and hybrid environments. To further demonstrate the utility of the package, the real gravitational-wave event GW170809 was analysed using Gravelamps; this event was found to yield no strong evidence supporting the lensing hypothesis, consistent with previously published results.


INTRODUCTION
For most of the time that we have spent observing the universe, we have been able to do so using only the electromagnetic (EM) spectrum -firstly via the visible window and then in the 20th century expanding to cover the full range from gamma rays to radio. However, with the dawn of gravitational-wave (GW) astronomy and the first detections of signals from the mergers of compact binary systems made by LIGO and VIRGO, we have now opened an entirely new gravitational-wave window on the universe. The detection of gravitational waves was the culmination of a century of research, from the first theoretical description of the phenomenon (Einstein 1916) in 1916 up to the observation of GW150914 (Abbott et al. 2016), and onwards -with now approximately 100 confirmed detections (Abbott et al. 2019) (The LIGO Scientific Collaboration et al. 2021) , and the prospect of many hundreds more as the extended 2nd-generation ground-based detector network reaches its design sensitivity later this decade.
In addition to the development of the ground-based detectors themselves, in recent years a great deal of work has gone into creating computational tools capable of generating predicted GW waveforms. These are compared to interferometer data, in order to detect sources and infer their parameters -as discussed in more detail in Abbott et al. (2020). Two extensive toolkits of particular importance are the LIGO Algorithm Library, or LAL (LIGO Scientific Collaboration 2018), and Bilby (Ashton et al. 2019): the former allows the straightforward generation of many different types of waveform, as just one of its many functionalities, whilst the latter presents an easy interface with which to perform the aforementioned parameter estimation using a variety of nested sampling (Skilling 2004) (Skilling 2006)

algorithms.
A recent focus of attention for the nascent field of GW astronomy has been the gravitational lensing of gravitational waves (The LIGO Scientific Collaboration & The Virgo Collaboration 2021) (Broadhurst et al. 2019) ) (Janquart et al. 2021) (Seo et al. 2021). Just as light passing by a massive compact object is deflected by the warping of space-time, so too gravitational wave signals are affected by the distribution of intervening matter as they propagate towards us. Electromagnetic observations of gravitational lensing have played an important role in the development of general relativity, with the deflection of light being one of its four major tests (Dyson et al. 1923). In a similar way, it is hoped that future observations of the lensing of gravitational waves may prove to be not just an important further test of general relativity and other theories of gravity, but also a powerful diagnostic probe of the nature of dark matter (Mishra et al. 2021) (Urrutia & Vaskonen 2021).
Whilst dark matter may not be probed directly using light, EM observations of gravitational lensing have been used to reconstruct the distribution of dark matter in the lensing object (Schneider 1996) (Amorisco et al. 2021) and a number of authors have recognised that GW gravitational lensing could also provide a valuable adjunct to these already extant investigations . Already, authors have investigated how the mass density profile will affect the properties of lensed GW signals (Cao et al. 2014) (Takahashi & Nakamura 2003), and sought to demonstrate that these properties may be used to characterise the lensed events -thus constraining the lens model parameters (Herrera-Martín et al. 2019) (Mishra et al. 2021). By considering different lens models, which correspond to differing kinds of lensing objects such as isolated point masses, dark matter halos, etc, analysis of the lens model therefore also provides insight into the nature of the lensing object.
Extending the above work, presented here is the package Gravelamps, designed to analyse lensed gravitationalwave data using arbitrary lens models. Built upon the LAL and Bilby frameworks, Gravelamps performs parameter estimation on GW events that have already been confidently identified as lensed, based on the application of such pipelines as e.g. Golum (Janquart et al. 2021) or hanabi (Lo & Hernandez 2021). Specifically, Gravelamps can analyse a range of different mass density profile models, yielding estimates for the joint posterior distribution of the lens and source parameters for each model. The code also allows comparison between different models through calculation of their relative Bayes Factor, which quantifies how probable one model is compared to another, thus providing a means to quantitiatively evaluate which mass profile for the lens is most likely.
In the following sections, the reader will be guided through the theory behind Gravelamps and its usage, performance, design sensitivities, as well as results from it and how it may be extended. Section 2 will cover the theory of gravitational lensing in general, Section 3 will cover the specific profiles that are presently integrated into Gravelamps, Section 4 will develop upon this to describe the calculation of the amplification factor, Section 5 will then go into more detail about the code itself -giving more detailed explanations of how it was constructed and how it can be used. Section 6 will then present some results from usage of Gravelamps in a variety of possible scenarios. Finally, Section 7 will conclude and give some of the aforementioned detail on possible future extensions, as well as integration of the code into existing analysis pipelines.

GRAVITATIONAL WAVE LENSING
For an overview of gravitational lensing theory, the reader is referred to works such as (Mollerach & Roulet 2002) or (Schneider et al. 1992). Here we give a very brief introduction to the gravitational lensing of GWs. In this case the relationship between a lensed signal, h L (f ), and an unlensed signal, h (f ), takes the form (Takahashi & Nakamura 2003): where the quantity F (f ) is called the amplification factor and is the ratio between the wave amplitudes of the observed lensed φ L obs (w, η) and unlensed, φ obs (w, η) signals: where w is the dimensionless frequency, and η as defined above is the position of the source. The simplicity of the relation between the lensed and unlensed signal is owed to the fact that with the exception of the frequency, the source parameters defining the unlensed signal and the lens parameters defining the amplification factor are uncorrelated. What does determine the amplification factor, however, is the mass density profile of the lens object. This will be explored in more detail in the next section.

LENSING PROFILES AND THE AMPLIFICATION FACTOR
As was noted in the introduction, lensing objects can be a number of differing things: from isolated point masses to extended objects like dark matter halos. One of the ways in which these differences manifest is in the distribution of the mass of the object given by ρ(x). These density profiles result in different lensing signatures (Takahashi & Nakamura 2003), which in principle means that the lensing signature can provide information on the density profile. Thus, in principle, it is possible to determine whether a model profile is consistent with observed lensed data thus meaning that such an analysis can place important constraints on the nature of the lensing object. This determination of the model profile is that Gravelamps seeks to make.
Whilst the density profile is different for each lens model and this ultimately yields different forms of the amplification factor, the underlying mathematics is the same. The astrophysical distances involved allow the use of the so called thin-lens approximation to be applied. This means that instead of requiring the full three dimensional mass density profile, ρ(x), the surface mass density, Σ(ξ) may be considered instead. This is the projection of the three-dimensional density onto a two-dimensional plane perpendicular to the line-of-sight to the source and at the distance of the centre of mass of the lensing object, i.e (Mollerach & Roulet 2002) Σ(ξ) = ρ (ξ, z) dz. ( Should the extended mass distribution be spherically symmetric, some further simplification may be made. In this case, the surface mass density depends only upon the modulus of the impact parameter, ξ = |ξ|. Using an arbitrary normalisation length, ξ 0 , one may construct dimensionless quantities from the impact parameter and the source position (Schneider et al. 1992): with D OL and D OS being the angular distances between the observer and the lens and source respectively. Whilst the value of the amplification factor differs between the lensing profiles, it can be caculated from a general expression using the above dimensionless quantities (Takahashi & Nakamura 2003) (Schneider et al. 1992): where T (x, y) is the dimensionless time delay. The simplest form of this general expression is for the case of axially-symmetric lensing object, where it is given by (Takahashi & Nakamura 2003): where ψ(x) is the lensing potential, φ m (y) is the phase for the minimum time delay -chosen such that the minimum time delay is zero -and w is the dimensionless frequency defined by (Takahashi & Nakamura 2003): In the case where the dimensionless frequency is very high (i.e. w >> 1), the geometric optics approximation may be made. In this case the stationary points of the dimensionless time delay are the only contributors to the amplification factor integral, yielding a simpler expression for the amplification factor (Nakamura & Deguchi 1999): where µ j denotes the magnification of the j t h image, T j = T (x j , y), and n j = 0, 1/2, 1 when x j is a minimum, saddle point, or maximum of T (x, y) respectively.

Point Mass
The simplest possible lens case is that of a point mass. It is applicable in the cases of compact objects such as individual stars or black holes. It can also be used for extended objects, by means of the Birkhoff Theorem, in the case where the lens size is much smaller than the Einstein radius. It is one of the only cases in which the full wave optics calculation of the amplification factor, i.e. Equation 6, can be carried out analytically. Here the normalisation constant ξ 0 introduced above corresponds to the Einstein radius and is given by: This leads to the dimensionless frequency being given by w = 4M lz ω where M lz is the redshifted lens mass. This ultimately leads to the amplification factor in the full wave optics case being given by (Takahashi & Nakamura 2003): where 1 F 1 is the confluent hypergeometric function of the first kind, and φ m (y) is given by (x m − y) 2 /2 − ln x m where x m = (y + y 2 + 4)/2. In the geometric optics approximation case, the amplification factor is given by (Takahashi & Nakamura 2003) where the magnifications are given by µ ± = 1/2 ± (y 2 + 2)/(2y y 2 + 4) and the time delay between the two images is given by ∆T = y y 2 + 4/2 + ln(( y 2 + 4 + y)/( y 2 + 4 − y).

Singular Isothermal Sphere
The Singular Isothermal Sphere, or SIS, profile is the simplest and most widely used profile that has been designed to model the behaviour most commonly observed for galaxies -i.e. a flat rotation curve. This behaviour as a result of modelling the galaxy as a large, extended object containing luminous matter embedded within a dark matter halo. Whilst the key strength of the SIS profile is the ability to replicate this flat rotation curve (Gavazzi et al. 2007), it does suffer from a weakness in the form of a central singularity which is non-physical (Burkert 1995) (Shi et al. 2021). The full SIS density profile is given by (Binney & Tremaine 1987): where σ is the one-dimensional velocity distribution of stars around the galaxy being modelled. The normalisation chosen for the SIS profile, similarly to the point mass case, is the Einstein Radius, which in the SIS case with velocity disperison, v, is given by (Takahashi & Nakamura 2003): The redshifted lens mass in the SIS profile is given by (Takahashi & Nakamura 2003): which results in a consistent definition of the dimensionless frequency, w as with the point mass.
The amplification factor may be calculated from Equation 6 by identifying that in the SIS case ψ(x) = x and φ m (y) = y + 1/2. Expanding the second term in the exponential into an infinite sum yields: This equation is solvable. Using the identity that e z 1 F 1 (a, b; −z) = 1 F 1 (a, b; z) (Abramowitz & Stegun 1965), this ultimately yields (Matsunaga & Yamamoto 2006): In the geometric optics case, the source position determines the number of images that will be formed. If y < 1 double images occur, where if y ≥ 1 only a single image is formed meaning that the amplification factor is split between these cases, being given by (Takahashi & Nakamura 2003): where in this case the magnficiations are given by µ ± = ±1 + 1/y and the time delays are given by ∆T = 2y.

Navarro, Frenk, White
In the standard Concordance Model of Cosmology Aghanim et al. (2020) Carroll (2001 the dark matter is assumed to be cold and collisionless. Consistent with these properties, the Navarro, Frenk, and White, or NFW, profile is widely used to model the density profile of cold dark matter (CDM) halos, as described by a singular 'universal' scaling function. Whilst allowing for more complexity than the Singular Isothermal Sphere profile, the NFW profile suffers from the same central singularity, or cusp. The density profile is given by (Navarro et al. 1997): where ρ s is the central density and r s is the characteristic scale of the profile. This characteristic scale is also chosen as the normalisation constant ξ 0 for this profile.
Owing to the increasing complexity of the profile, the form of the lensing potential is more complicated being given by (Bartelmann 1996): Here, κ s is the dimensionless surface density which is given by 16πρ s (D OL D LS /D OS )r s . The time delay constant also becomes more complicated as the lens equation no longer has an analytical solution meaning it must be calculated numerically, meaning the amplification factor too must be calculated numerically from Equation 6 in the wave optics case, or Equation 8 in the geometric optics case. In the latter case, the number of images formed is dependent upon the value of y relative to a critical value -with three images being formed where |y| < y cr , and only one in the case where |y| > y cr . In the case of the source position being precisely at the critical value, the magnification becomes infinite.

GRAVELAMPS OVERVIEW
Gravelamps is an analysis pipeline designed to study arbitrary lensed gravitational-wave data and constrain properties of both the source and mass density profile of the lens . The pipeline estimates the Bayesian evidence for the user specified lens model. In so doing it also allows given that chosen model inference of both the source and lens parameters. By performing this analysis for multiple lens models, comparison may be made between these evidences to determine the most likely model for the data. The application of Gravelamps to some example scenarios, to illustrate its usefulness for identifying and characterising lensed gravitational-wave signals, will be presented in the following sections.

Languages and Toolkits Used
For the higher level parts -in terms of interacting with the user and performing the parameter estimation, Gravelamps is coded in python. This allows for flexibility and ease of reading as well as access to a number of useful modules. In use in Gravelamps are the extensively used numpy (Harris et al. 2020), scipy (Virtanen et al. 2020), and astropy (The Astropy Collaboration et al. 2018) modules. Most importantly, the use of python gives access to the bilby package (Ashton et al. 2019) thus ensuring the use of both reliable parameter estimation methods and native and easily extensible support for LAL waveforms. This creates a more cohesive experience for a user wishing to modify the program beyond the level of the options already available to them. The use of bilby also allows the pipeline to be scaled straightforwardly for usage on clusterised machines, through bilby pipe (Romero-Shaw et al. 2020). In particular, there is little difference between a run designed to be implemented on a local machine and one designed for a clusterised machine, e.g. to be submitted through the scheduler HTCondor, as is used widely in the LSC clusters.
For the lensing data generation, as can be seen in the previous section, the amplification factor calculation is not a simple task as it requires the calculation of the integral of a highly oscillatory function. Due to the intensive nature of the calculations, these were coded in C++ to benefit from the additional speed of a compiled language -as well as granting access to the C arbitrary precision library arb (Johansson 2017). In addition to being able to perform the necessary calculations with appropriate speed and accuracy, this approach had the benefit of allowing users to decide for themselves how much they wished to trade said speed and accuracy, as well as allowing the full wave optics to be calculated to as high a dimensionless frequency value as possible -well past the point at which the geometric optics approximation can successfully take over -without causing disruption to the amplification factor that is calculated.
In addition to those toolkits used directly in Gravelamps, there are other pieces of software which are critical to the functioning of the toolkits themselves -with the bilby package additionally depending upon such packages as matplotlib (Hunter 2007) and corner (Foreman-Mackey 2016) for plotting, and gwpy (Macleod et al. 2021) and lalsuite (LIGO Scientific Collaboration 2018) for gravitational wave data analysis. Further, through bilby, Gravelamps allows the use of a great many nested samplers: chosen for the work presented here was dynesty (Speagle 2020).

Design Intentions
As with almost all software packages, there were certain philosophies which underpinned how the software was constructed, in addition to the obvious requirements of functionality and speed. In the case of Gravelamps, the following were the most important considerations: 1. Openness: With such complicated integrals being necessary to perform the calculations to generate the amplification factor, it is tempting to consider platforms such as Mathematica which are designed to tackle extremely complex mathematics easily and with speed. However, a problem with such platforms is their proprietary nature.
Gravelamps is designed to be open for use in as many applications as possible without any licensing conflicts; as such the toolkits used were chosen to reflect this, as well as the license chosen for Gravelamps itself -with the source code for all of these being freely available (Wright 2021).

Simplicity:
In the spirit of the openness with which Gravelamps was designed, it was also designed to be as understandable to the user as possible. In particular, the code was designed with ease of readability in mind, and so that interactions with the software by the user should be as simple as possible. As such, there is little difference between a run on a local machine and that on a cluster -requiring only manipulation of the configuration INI file.
3. Extensibility: Whilst at this time, only the three density profiles discussed in detail in this paper are fully coded, Gravelamps as a platform is designed to be easily extensible -for example by including other density profiles or by considering multiple lenses.

Structure of an Analysis Run
Having been designed to be relatively simple, Gravelamps does not require that the user in any way modify the backend, in order to carry out an analysis run; instead the user may select all necessary settings through means of a simple INI file -contained within which there are options governing: • The output information • Lensing settings -such as the start, end and, number of points for the amplification factor grid if it being used as well as those settings pertaining to arithmetic precision • HTCondor settings for clusterised runs • bilby pipe settings for clusterised runs • Common analysis options such as the waveform model adopted (including whether to use direct or interpolated calculation of the amplification factor), trigger time chosen, etc • Whether or not to perform an unlensed analysis, and if so, the associated settings to be used • Event/Injection parameters

• Sampler settings
Here the user is given as much control over the run as possible without the requirement of delving into the code and modifying it directly for each analysis -although the latter is, of course, possible too. With a completed INI, the user need only call one of the preconstructed analysis functions, gravelamps local inference or gravelamps pipe inference, depending on whether they wish the run performed locally, or they wish submit files for the HTCondor Scheduler prepared.
Should users be implementing a clusterised version of the code, they will receive a DAG file that will call each of the necessary submit files in order. If running locally, the steps will be carried out directly. From the point at which the analysis run begins, the procedure is as follows: 1. Over the specified ranges of the dimensionless frequency and source position, generate a grid of values of the amplification factor for the chosen lens model 2. Generate an interpolating function to evaluate the amplification factor for any given dimensionless frequency and source position within the specified ranges 3. Generate injection data/Fetch event data -i.e. prepare the data that will be analysed 4. Run parameter estimation on the data using the chosen lensed waveform model In the case where the user is operating solely in the geometric optics regime, to assist with computational efficiency they may also choose to have the amplification factor calulated directly in the waveform, bypassing the need for steps 1 and 2.
The procedure outlined above results in acquiring both estimates of the lens and source parameters for the given model, and also an estimate of the evidence for that model. This can then be compared with the evidences for other models to give a quantitative indication of the preferred model. This latter step may include, if the user wishes, an estimate of the evidence for an unlensed model. In this unlensed case, however, Gravelamps should not be considered as providing a true Bayes factor for the relative probability of the lensed and unlensed hypotheses since no explicit evaluation of any selection effects relevant to the lensed-to-unlensed posterior odds ratio is carried out, based on e.g. a numerical injection study that models the probability of 'false positives' -i.e. pairs of candidate gravitational-wave events with parameters that could be consistent with them being strongly-lensed multiple images of a single source (The LIGO Scientific Collaboration & The Virgo Collaboration 2021). In its current form, Gravelamps conducts lens and source model parameter estimation on only a single such image. In the future we plan to explore its extension to the framework appropriate for analysis of multiple images, but in the remainder of this work we consider only the single image case.

Package Structure
Following the standard python package layout of submodules with specific purposes in mind, Gravelamps is split into gravelamps.inference and gravelamps.lensing covering the performance of the inference runs and the lensed waveform generation respectively.

gravelamps.inference
Contained within gravelamps.inference are the individual scripts that form the programs by which the user interacts with Gravelamps, i.e. gravelamps local inference and gravelamps pipe inference. In addition to the programs themselves, contained with gravelamps.inference are all of the functions which are concerned with handling the user's configurations as well as the generation of the HTCondor submit files for the clusterised code. These are contained within the helpers and file generators parts of the submodule.

gravelamps.lensing
gravelamps.lensing contains within it those parts of the source code pertaining to the construction of lensed waveforms.
It is itself split between a C++ and a python component. It also in a similar manner to the gravelamps local inference or gravelamps pipe infernece contains a pair of programsgravelamps generate lens local and gravelamps generate lens pipe that will generate lensing data without then constructing a lensed waveform and performing analysis on it.
The C++ parts of this submodule contains the source codes for the functions that calculate the amplification factor values, as well as the program that generates amplification factor data from which an interpolator may be constucted. These programs are compiled upon the installation of Gravelamps. Each of these programs functions in a similar maner: taking in as inputs dimensionless frequency and source position information contained within unique filestogether with user defined settings, such as the dimensionless frequency value at which to begin using the geometric optics approximation, or the arithmetic precision of the wave optics calculations. Eaah program then returns two files containing two unique files containing the amplification factors' real and imaginary components, which can then be used to generate a complex interpolator for the amplification factor function.
The python parts of the submodule concern both the utility functions that pertain to the construction of the dimensionless frequency and source position data used by the above. Most importantly, it also contains those functions that construct the lensed waveforms over both the wave and geometric optics regimes. This is done by firstly either using the data generated by the C++ programs to generate an interpolating function for the amplification factor, or should the user wish to work exclusively in geometric optics they may also instruct gravelamps.lensing to instead, by means of the ctypes module contained with python, directly use the calculation function from the C++ part. Once the amplification factor calculation function is established, it will then generate a base LAL waveform and using the amplification factor lens this waveform. Finally, within the python section, are those functions which perform physical calculations such as unit conversion, etc, whereas those calculations in gravelamps.inference are more statistical in nature.
Contained within the python parts of gravelamps.lensing are the utility functions that pertain to the construction of the lensed waveform directly, such as generating the dimensionless frequency and source position data over which interpolation will take place, as well as taking the base unlensed LAL waveforms and turning them into the corresponding lensed versions. In addition, gravelamps.lensing contains those functions which perform physical calculations such as unit conversion, etc, whereas those calculations in gravelamps.inference are more statistical in nature.
The ability to implement both the wave optics and geometric optics calculations is particularly useful, in view of the complexity of the calculations required for the full wave optics analysis. Whilst these can be run for higher values of the dimensionless frequency, the calculation requires ever increasing precision and greater numerical limits on the integrations and summations; this leads to an increasing slow down due to the necessity of performing more calculations to reach these higher thresholds. However, Gravelamps has been written in a sufficiently flexible state to leave this choice to the user rather than making it for them.

EXAMPLE USES OF GRAVELAMPS
One of the major focuses of Gravelamps is to be as versatile in terms of how it can be used to study lensing of gravitational wave signals. As such, the lensing data can be generated by themselves, i.e. simulated data can be analysed, as well as real data. The following will cover examples of each of these use cases to demonstrate how an analysis might be made.

Generation of Amplification Factor Data
As was discussed in the previous sections, the generation of a lensed waveform from an unlensed one is simply a matter of multiplying the strains of the unlensed waveform by the amplification factor function corresponding to the lensing model of interest. With the base, unlensed, waveforms being implemented by LAL, Gravelamps generates the amplification factor to then lensed these waveforms. As was also previously noted, there are two functions that can caclulate the amplification factor: the full wave optics approach or the geometric optics approximation. Both of these have benefits and drawbacks, and both are implemented within Gravelamps.  Figure 1 shows the results of using Gravelamps lens generation codes for both the wave and geometric optics cases, up to a dimensionless frequency of w = 100 over a range of source position values, for each of the Point Mass, SIS, and NFW lens models to illustrate the differences and applicability of both calculation styles.
In the case of full wave optics, the calculations are obviously more complete. However, they are significantly more computationally intensive, requiring the use of arbitrary precision calculations as well as the use of numerical integration in the most complex case. Even the SIS -the simplest of the extended object models requires the approximation of an infinite sum to a finite degree. This makes the generation of the wave optics amplificaiton factor slower -increasingly so with increasing dimensionless frequency where in order to retain numerical stability of the result, greater limits need to be placed on the integrations as well as increasing precision.
In the case of the geomtric optics approximation, the calculations are much simpler -resulting in sufficient speed that they may be used directly in the waveform generation as opposed to the interpolator approach used in the wave or mixed approaches. However, strictly this approximation is only applicable at higher dimensionless frequency values as at lower dimensionless frequency values they do not replicate the behaviours of the full wave optics function -noting that in the SIS and NFW cases the damped oscillatory nature of the higher source position value curves are replaced with a signle value that approximates the average of the oscillation. The increase in speed however, does warrant the use of geometric optics in higher dimensionless frequency regimes.
Gravelamps was designed to be a versatile system however, so the user is given as much control over these calculations as possible. The user from the INI file used to generate the amplification factor data is able to specify firstly whether to use the direct or interpolation methods. The former may be used for geometric optics alone, where the latter may be used for wave optics alone, or a mix of wave and geometric optics, (or indeed also for geometric optics alone). In the case of using the interpolator method, the INI file used to generate the amplification factor data is able to specify: • The minimum and maximum dimensionless frequency and source position • The number of dimensionless frequency and source position to be generated for the interpolator • The precision used in the arbitrary precision calcualtions (for wave optics only) • The upper limit of the summation and integration used in the SIS and NFW models respectively These choices are given to the user to allow them to decide for themselves how much speed/accuracy trading that they wish to adopt. This also allows for a hybrid approach, in particular preventing the user from being forced to use e.g. geometric optics in all cases; if the user has the computational resources and time available they are able to continue using the wave optics approach to as high a dimensionless frequency as they wish. It should be noted, however, that the higher dimensionless frequency wave optics calculations become ever more computationally expensive in order to retain numerical stability. This increase becomes particularly noticable between w = 100 and w = 1000. It is therefore recommended that should the user need to use wave optics above these values that they do so only when they are planning to use a single lens generation for multiple parameter estimation runs where the justification of such computational expense may be split between the runs.
The flexibility of Gravelamps extends to other features. The user is able to specify any terminal-callable process in the INI and also what arguments are needed, allowing the user the ability to generate new models easilyt, even without adding them to Gravelamps itself directly (although users are encouraged to incorporate any such new models into the codebase should they wish to do so). If the user chooses to use non-implemented models, they are able to use the gravelamps generate lens programs to generate amplification factor data. Placing the locations of the resultant files into the INI's optional input then allows these data to be accessed in inference runs.

Analysis of a Simulated Lensed GW150914-like Event
In this section we use Gravelamps to simulate a lensed waveform, with source properties modelled on those of GW150914, as an extension of the waveform generation capabilities of bilby -themselves wrappers for LAL -to give a wide variety of base, unlensed waveforms to which the generated amplification factor data can be applied. Selected here as an illustrative example is the IMRPhenomXP waveform model (Pratten et al. 2021).
To make Gravelamps more feature complete, it is also capable of simulating an event as one waveform, and then analysing it as another -allowing investigation of the case where the wrong lensing profile is mistakenly applied. Moreover, as already highlighted, Gravelamps is capable of generating predicted lensed waveforms under both geometric optics and wave optics scenarios.
The first analysis presented here simulates a GW150914-like event which is microlensed by a 50M mass point lens, and is correctly identified as a point-lensed event. The second analysis examines a similar scenario with a 1000M lens to show a more likely detection scenario. The final analysis explores the limitations of geometric optics alone, by showing that a 4.4 × 10 7 M Singular Isothermal Sphere, replicating a galaxy-scale 'macrolensing'-type strong lensing event, analysed purely using geometric optics is unable to discriminate conclusively between SIS and NFW lens models. In all cases, the lens is placed half-way between the observer and the source, with the source position, y, being set to 1.
The true source properties used in each of these illustrative cases are summarised in Table 1 below. Gravelamps has attempted to recover all of these parameters with the exception of the geocenter time. Discussed and presented below are a subset of these fuller results considering the most obviously related parameters -the masses, distance, and lens parameters (source frame lens mass and source position). Full configurations for the analyses presented here may be found within the Gravelamps repository Wright (2021) and were performed using the version located at Wright (2022) with the exception of those runs for Figure 4 which were performed with an updated version due to the identification of a bug specific to those runs during the review process. Figure 2 shows a subset of the full resulting parameter estimates inferred from the Gravelamps microlensing simulation when analysed with the correct model -i.e. an isolated point lens. As can be seen, each of the parameters has been recovered well, being tightly constrained to ranges consistent with their true values. This follows into the log Bayes Factor estimates -when compared with noise, analysis as a point lens, yields a value of 11063.213 ± 0.332. When analysis this simualated event as an SIS lensing event or a NFW lensing event with k s = 2 the Bayes Factors compared with noise are 11061.194 ± 0.334 and 11054.524 ± 0.345 respectively. This leads to an overall favouring of the Point Lens by a log Bayes Factor of 2.019 over the SIS case and 8.689 over the the NFW case. In the more likely scenario of a lensing mass on the scale of 1000M -specifically considered here was the case of M lens = 1000M , the results are even more conclusive and are presented in Figure 3. In this case, the log Bayes Factor comparing the signal to the noise case was 9208.585 ± 0.356 for the true point lens case and 9184.879 ± 0.352 when analysing the simulation as an SIS event -yielding a favouring of the point lens scenario with a log Bayes Factor of 23.706. Figure 4a shows the resulting parameter estimates from the Gravelamps SIS macrolensing simulation using purely geometric optics. As can be seen, the parameter estimation performance has decreased, particualrly amongst the lensing parameters. The broadening of the lens mass posterior is particularly noticeable. This is an expected result, as the inference of ther lens mass is dependent upon the change of the amplification factor as a function of the dimensionless frequency (in which the lens mass is a contributing quantity). In the geometric optics case, the amplification factor becomes invariant over all dimensionless frequency, and the lensed signal thus becomes insensitive to the lens mass. Figure 4b shows the resulting parameter estimates in the case where the data is analysed under the incorrect lens model. Of particular note is that there is little difference between the parameter estimates in this case. This is reflected in the log Bayes Factors: when analysing the data with the true SIS model, the log Bayes Factor comapred with noise is 8495.794 ± 0.241 while for the incorrect NFW case the log Bayes Factor is 8495.065 ± 0.242. This yields a difference in log Bayes Factor of 0.73 -i.e. indicating that, when considering geometric optics alone, the two lens models are insufficiently distinguishable.

Analysis of Real Event GW170809
To further illustrate of the usefulness of Gravelamps and its suitability for analysis of real gravitational-wave event candidates, the real gravitational-wave event GW170809 was next analysed under the hypothesis that this event was lensed by a Point Mass and SIS profile respectively -with prior ranges for each model reflecting the microlensing hypothesis. This event was chosen due to the fact that, of those events analysed from the second observing run -from which the data is available from the Gravitational Wave Open Science Center (GWOSC) (Rich ) -GW170809 was identified as one of the strongest (albeit ultimately rejected) lensing candidates . Consequently, there has also been some additional interest shown in this event as a potential lensed candidate (Broadhurst et al. 2019). Figure 5 shows the result of the Gravelamps runs for GW170809 analysing in the microlensing regime for the point (5a) and SIS (5b) models respectively. In this case, the log Bayes Factors -comparing the hypothesis of modelled signal + noise, versus that of pure noise -for the point and SIS cases are 56.773 ± 0.265 and 57.138 ± 0.258 respectively. We see, therefore, that the signal + noise hypothesis is strongly favoured -consistent, of course, with the fact that GW170809 was identified as a confident gravitational-wave detection on GWTC-1. In this case, however, the lensing parameter estimates are particularly broad as compared with the simulations which may be an indication that the lensing models are ill-fitting to the data, an uninformative posterior in the lensing parameters being what one would expect in the case of an unlensed signal. Thus, while our Gravelamps analysis does marginally favour the SIS lens Figure 2. Subset of results of lens and source model parameter estimation, performed by bilby nested sampling using dynesty, for a GW150914-like GW event microlensed by a 50M mass point lens. Parameters shown are from left-to-right: the chirp mass (M), the mass ratio (q), the luminosity distance (dL), the source frame lens mass, and the source position. As can be seen, Gravelamps is able to successfully recover the model parameters.
model over the point lens model with a log Bayes Factor of 0.365, these results should be taken as merely illustrative of the capabilities of the software rather than any indication of a preference for a particular lens model for GW170809. We emphasise again that Gravelamps does not explicitly take into account the selection effects, multiple image analysis, and more detailed population prior choices that are explored in (The LIGO Scientific Collaboration & The Virgo Collaboration 2021), in order to better assess whether there is evidence supporting the lensing hypothesis for any given gravitational-wave candidate event. Nevertheless, we believe that the example of GW170809, and the previous simulated examples, illustrate the efficacy of Gravelamps for comparing the evidence for different lens models given that a lensed gravitational-wave event has been detected.

CONCLUSION
With the increasing sensitivity of the existing ground-based detectors, and the additional detectors set to join the global network in the future, the detection of a lensed gravitational-wave event within a matter of a few years would appear to distinctly possible. Those searches performed thus far have focused on identifying candidate lensed events;  Figure 3. Subset of results from lens and source model parameter estimation, performed by bilby nested sampling using dynesty, for a GW150914-like GW event macrolensed by a 1000M point lens, analysed as a point and an SIS event respectively. As can be seen, Gravelamps has again successfully recovered the model parameters in the case of the correct model.  . Subset of results from lens and source model parameter estimation, performed by bilby nested sampling using dynesty, for a GW150914-like GW event macrolensed by a 4.4 × 10 7 M SIS lens, analysed as an SIS and an NFW macrolensed event respectively. As can be seen, the performance of the lens parameter estimates has decreased due to geometric optics encoding less information in the signal, however, the source parameter estimates are still constrained however, once such a candidate event has found, the immediately important question becomes what is the astrophysical nature of the lens itself? The package Gravelamps presented here has been designed to help answer that question.
Gravelamps has been designed to be an easy to use and versatile platform, with the flexibility to allow investigation of both microlensing scenarios and the so-called intermediate 'mililensing' case. Gravelamps is particularly adapted to this latter case by having the flexibility to calculate the lensing amplification factors in either geometric or wave optics, or in a hybrid combination of both, over the course of a single analysis run -with full control given to the user as to where each regime may apply. As an initial, illustrative, set of lensing models, Gravelamps presently contains the point mass, Singular Isothermal Sphere, and Navarro, Frenk, White density profile; however it has been designed to easily be extended to include any lens model that the user wishes to investigate.
We have demonstrated the utility of Gravelamps by presenting examples, in both of the wave and geometric optics cases of the results of parameter estimation. Here we demonstrated that one of the major limits of the geometric optics case is that alone it is insufficient to determine the lensing model -at least within the single image analysis performed here. However, scenarios which may partially incorporate geometric optics cases alongside wave optics effects -looking for lensing in the O(1000M ) region for instance -are identifiable.
As noted, Gravelamps is an easily extensible platform and the work presented here is seen as the just the first step in its development. For example, in future we plan to explore its extension to handle multiply lensed signals -for instance the case of a source that simultaneously undergoes both microlensing and macrolensing.