lensingGW - a Python package for lensing of gravitational waves

Copyright © 2020 Giulia Pagano

lensingGW is a software package to simulate lensed gravitational waves in ground-based interferometers from arbitrary compact binaries and lens models.

Its algorithm allows to resolve strongly lensed images and microimages simultaneously, such as the images resulting from hundreds of microlenses embedded in galaxies and galaxy clusters.

Please submit bugs/suggestions etc to giulia.pagano@outlook.com.


Installation & requirements

  • lensingGW requires a forked version of Lenstronomy. The following commands download it in the current folder and install it

    $ git clone https://github.com/gipagano/lenstronomy.git
    $ cd lenstronomy
    $ python setup.py install
    
  • lensingGW can then be installed as follows:

    $ cd ..
    $ git clone https://gitlab.com/gpagano/lensinggw.git
    $ cd lensinggw
    $ python setup.py install
    
  • It requires python>=3.7

Quickstart

This section illustrates lensingGW’s core functionalities through worked-out examples. They demonstrate the solver usage and the simulation of the lensed and unlensed gravitational signals.

The lens model is specified through a list of lens profile names and a list of related parameters. It is handled through a forked version of Lenstronomy (original version here): any lens profile implemented therein can be used.

Ready-to-use python scripts implementing the following examples are provided along with each subsection and at the end of the section.

Solve the lens model

Radians

We first define the lens configuration: we use radians, which is lensingGW’s default. As an illustration, we consider two point lenses with Einstein radii \(\theta_{E_1}\) and \(\theta_{E_2}\), whose parameters we have defined in kwargs_lens_list

# lens model
lens_model_list  = ['POINT_MASS', 'POINT_MASS']
kwargs_lens_list = [{'center_x': 6.9e-11,'center_y': 0.0, 'theta_E': thetaE1} , {'center_x': -6.9e-11,'center_y': 0.0, 'theta_E': thetaE2}]

We then specify the solver settings

# solver setup
solver_settings = {'SearchWindowMacro': 4*thetaE1,
                   'SearchWindow'     : 4*thetaE2}
The solver implements the two-step-procedure presented in Pagano et al. 2020. The minimal setup to specify consists in the extension of the search region for the macromodel, SearchWindowMacro and for the complete model, SearchWindow.
Finally, we import the solver routine and solve for the image positions
from lensinggw.solver.images import microimages

# solve for the image positions with the two-step procedure
Img_ra, Img_dec, MacroImg_ra, MacroImg_dec, pixel_width  = microimages(source_pos_x    = 1.39e-11,
                                                                       source_pos_y    = 1.20e-10,
                                                                       lens_model_list = lens_model_list,
                                                                       kwargs_lens     = kwargs_lens_list,
                                                                       **solver_settings)
We obtain the right ascensions and the declinations of the images of the complete model, of the macromodel and the pixel width of the last iteration, in the same units as the input.
The description of the each iteration step can be enabled through the Verbose option in the solver settings.
solver_settings.update({'Verbose' : True})

The verbosity is a useful diagnostic in non-converging cases, as it allows to identify the problematic part of the iteration.

Image properties

Magnifications, time delays and Morse indices can be accessed through the related modules

from lensinggw.utils.utils import TimeDelay, magnifications, getMinMaxSaddle

# time delays, magnifications, Morse indices
tds = TimeDelay(Img_ra, Img_dec,
                source_pos_x = 1.39e-11, source_pos_y = 1.20e-10,
                zL = 0.5 , zS = 2.0,
                lens_model_list, kwargs_lens_list)

mus = magnifications(Img_ra, Img_dec, lens_model_list, kwargs_lens_list)
ns  = getMinMaxSaddle(Img_ra, Img_dec, lens_model_list, kwargs_lens_list)

By default, the solver considers the first profile in lens_model_list as macromodel and applies the two-step procedure to all the macroimages. It sets the same precision requirement on both types of images, as well as the same number of pixels and the same grid sizes for the iterations. These items can be customized through the solver settings. A detailed description of the available options can be found in the solver technical notes.

Arbitrary units

lensingGW supports arbitrary coordinate units, defined as \(x_{\mathrm{a.u.}} = x_{\mathrm{radians}}/ \alpha\), where \(\alpha\) is a constant. To use arbitrary units, all input must be given in the scaled coordinate system and the scaled flag must be enabled in the routines involving coordinates or conversions to physical units.
The previous example can be modified to use \(\alpha = \theta_{E}\), where \(\theta_{E}\) is the Einstein radius of the total mass: replace each quantity in radians with its scaled equivalent and modify the solver settings and the TimeDelay input as follows
# solver setup
solver_settings = {'Scaled'           : True,             # indicate that the input is in scaled units
                   'ScaleFactor'      : thetaE,           # and the scale factor
                   'SearchWindowMacro': 4*thetaE1/thetaE, # scaled input
                   'SearchWindow'     : 4*thetaE2/thetaE} # scaled input
# time delays, magnifications, Morse indices
tds = TimeDelay(Img_ra, Img_dec,                                                # already in scaled units
                source_pos_x = 1.39e-11/thetaE, source_pos_y = 1.20e-10/thetaE, # scaled input
                zL = 0.5 , zS = 2.0,
                lens_model_list, kwargs_lens_list,
                scaled       = solver_settings['Scaled'],      # indicate that the input is in scaled units
                scale_factor = solver_settings['ScaleFactor']) # and the scale factor

The magnification and Morse indices computation doesn’t need modifications, provided that all radians input has been scaled. This includes the lens parameters and the source position.

Multi-component macromodel

By default, lensingGW considers the first profile in lens_model_list as macromodel. However, users can indicate the profiles that form the macromodel by specifying the matching indices in the lens list. This allows to indicate a different profile or multi-component macromodels.
For example, we can define a 2-component macromodel constituted by both the point lenses of the previous examples by adding the following specification in the solver settings
solver_settings.update({'MacroIndex' : [0,1]})

If the number of specified indices corresponds to the number of lens profiles, all the deflectors are considered as part of the macromodel and only the macroimages are investigated.

Only macromodel

If the lens model is composed of a macromodel and a micromodel (such as multiple microlenses), users can require the study of the macromodel only by enabling an appropriate flag in the solver settings

solver_settings.update({'OnlyMacro': True})
This functionality, combined with the multi-component choice, allows to define the lens model once and for all and then investigate the impact of each component on the resulting images by acting on the solver settings.
For example, one can
  • enable OnlyMacro and change the lens profile/profile group indicated as macromodel - assesses the impact of the macromodel components on the images of the macromodel

  • disable OnlyMacro and change the lens profile/profile group indicated as macromodel - assesses the impact of the macromodel components on the images of the complete model

Simulate gravitational waves

Configuration files

If gravitational-wave (GW) signals are requested, a configuration file specifying the binary parameters and the waveform settings has to be provided to lensingGW. A template configuration file with explanations of each entry and a ready-to-use configuration file are provided in examples/ini_files on GitLab.

Note on angles

GWs from compact binaries depend on a set of angles. In addition to the binary right ascension and declination, the waveform depends on the inclination of the system, the polarization angle and on the phase of coalescence.
The standard choice in gravitational-wave astronomy is to indicate such angles in radians. To avoid confusion, the same choice has been adopted within lensingGW. This only affects the angles specified in the waveform configuration file: although any routine of lensingGW supports the definition of arbitrary units for the source coordinates, the source right ascension and declination in the configuration file must by expressed in radians.

Interferometers’ noise curves

Power spectral densities (PSDs) have to be provided for each interferometer in which the waveform is requested. Advanced LIGO and Advanced Virgo design PSDs released by the LIGO/Virgo Collaborations are provided in examples/psds on GitLab.

Get the unlensed gravitational waves

To generate unlensed waveforms, lensingGW’s signal model has to be initialized through the configuration file

from lensinggw.waveform.waveform import gw_signal

# read the waveform parameters
config_file = 'ini_files/waveform_config.ini'

# instantiate the waveform model
waveform_model = gw_signal(config_file)

Polarizations, strains and their frequencies can be accessed through lensingGW’s unlensed routines

# compute the unlensed waveform polarizations, strains in the requested detectors and their frequencies
freqs, hp_tilde, hc_tilde, strain_dict = waveform_model.unlensed_gw()

# and their signal-to-noise-ratios
SNR_dict = waveform_model.unlensed_snr()

The strain and SNR dictionaries allow to access information of each detector individually

# access an unlensed strain
sH1 = strain_dict['H1']

Get the lensed gravitational waves

In addition to the waveform configuration file and to the noise curves, the simulation of lensed signals requires to specify a lens model and a set of images, from which the amplification factor is computed.
Images can be found as described in the solver section, where the lens model initialization is also shown.
Assuming that a lens_model_list and a kwargs_lens_list have been defined and that the Img_ra and Img_dec arrays store the images’ right ascensions and declinations, one can modify the previous example as follows to get the lensed quantities
# compute the lensed waveform polarizations, strains in the requested detectors and their frequencies
waveform_model.lensed_gw(Img_ra, Img_dec,
                         source_pos_x = 1.39e-11, source_pos_y = 1.20e-10,
                         zL = 0.5 , zS = 2.0,
                         lens_model_list,
                         kwargs_lens_list)
# and their signal-to-noise-ratios
waveform_model.lensed_snr(Img_ra, Img_dec,
                          source_pos_x = 1.39e-11, source_pos_y = 1.20e-10,
                          zL = 0.5 , zS = 2.0,
                          lens_model_list,
                          kwargs_lens_list)

Arbitrary units

If the input is in arbitrary units, the scaling has to be notified to the lensed waveform routines by enabling the scaled flag and indicating the scale factor. The previous lines can be modified to use \(\alpha = \theta_{E_1}\) as follows

# compute the lensed waveform polarizations, strains in the requested detectors and their frequencies
waveform_model.lensed_gw(Img_ra, Img_dec,
                          .
                          .
                         scaled = True,
                         scale_factor = thetaE1)
# and their signal-to-noise-ratios
waveform_model.lensed_snr(Img_ra, Img_dec,
                           .
                           .
                          scaled = True,
                          scale_factor = thetaE1)

where omitted lines are unchanged.

User-defined cosmology

By default, a flat cosmology with \(H_0=69.7, \Omega_0=0.306, T_{cmb0}=2.725\) is considered to convert the lens and source redshifts \(z_L\) and \(z_S\) to angular diameter distances. However, users can indicate a different one by specifying an instance of the astropy cosmology class to TimeDelay, lensed_gw and to lensed_snr through the cosmo keyword argument
# time delays, magnifications, Morse indices
TimeDelay(Img_ra, Img_dec,
             .
             .
            cosmo = asropy_cosmo_instance)
# compute the lensed waveform polarizations, strains in the requested detectors and their frequencies
waveform_model.lensed_gw(Img_ra, Img_dec,
                          .
                          .
                         cosmo = asropy_cosmo_instance)
# and their signal-to-noise-ratios
lensed_SNR_dict = waveform_model.lensed_snr(Img_ra, Img_dec,
                                             .
                                             .
                                            cosmo = asropy_cosmo_instance)

where omitted lines are unchanged.

Technical notes

lensingGW has a modular structure which allows to import each routine individually. Detailed package description is available at the lensingGW package section of the documentation. Here, we provide a few technical notes on the packages that require special attention.

The amplification factor package

This package computes the geometrical optics amplification from a set of images and lensing potentials.

Unlike in the light lensing, gravitational waves from compact binaries are transient signals. Thus, only the images that superimpose into the detectors should be fed to the amplification routine. This amounts to considering groups of images whose pairwise time delays are shorter than the unlensed signal duration. The time delay of each image can be checked through the TimeDelay module in the utils package.

The solver package

The solver package finds the image positions of lensed compact sources by specifying a number of settings. It implements the two-step procedure presented in Pagano et al. 2020: the minimal solver setup requires to specify the extension of the search windows for the macromodel and for the complete model. However, each part of the iterative procedure is customizable by enabling the \(\texttt{optimization}\) mode.
Tuning of the solver settings is also recommended when no images are found. As the number of pixels and grid size of each iteration determine how well the lens potential is resolved, users may have to adapt them to the specific model under consideration.

This section contains the complete list of solver settings that can be tuned by the users along with their descriptions and default values. For the sake of clarity, here they are grouped in:

  • standard settings – specifications that are always considered for the two-step procedure. If users don’t specify input for these settings, default values are assumed

  • optimization-related settings – settings that are considered only when the \(\texttt{optimization}\) mode is active. Users who wish to use these settings should enable the \(\texttt{optimization}\) mode.

  • further specifications – additional features.

The settings are specified through a dictionary of keywords arguments, as shown in the solver usage examples.

Standard settings

For most cases, the \(\texttt{optimization}\) mode is not required. The following settings apply the solving algorithm in its standard form.

Option

Type

Default

Description

Scaled

bool

False

Specifies if the input is given in arbitrary units. If
\(False\), radians are assumed

ScaleFactor

float

1

Constant factor \(\alpha\) to convert from radians to
arbitrary units, as per \(x_{a.u.} = x_{radians}/\alpha\)

MacroIndex

list

[0]

Indices matching the profiles in
lens_model_list that form the macromodel

SearchWindowMacro

float

needs user input

Size of the first macromodel grid

PixelsMacro

int

\(10^3\)

Number of pixels of the first macromodel
grid

OverlapDistMacro

float

\(10^{-15}\rm{rad}\)

Distance below which macroimages are
considered overlaps

ImgIndex

list
or
None

None

Indices of the macroimages to be considered for
the two-step procedure. If None,
all macroimages are considered

SearchWindow

float

needs user input

Same as SearchWindowMacro, but for
the complete model

Pixels

int

\(10^3\)

Same as PixelsMacro, but for
the complete model

OverlapDist

float

\(10^{-15}\rm{rad}\)

Same as OverlapDistMacro, but for
the complete model

PrecisionLimit

float

\(10^{-20}\rm{rad}\)

Precision of the solutions in the
source plane

Optimization

bool

False

Enables the \(\texttt{optimization}\) mode

Verbose

bool

False

Prints detailed diagnostic

Entries for SearchWindow and SearchWindowMacro must necessarily be specified by the user.
Additional tuning is possible through the \(\texttt{optimization}\) mode, which is particularly useful in non-converging cases.

Optimization-related settings

When the \(\texttt{optimization}\) mode is active, users can tune the following additional parameters

Option

Type

Default

Description

OptimizationWindowMacro

float

2

Multiplying factor \(\gamma\) for the size
of the macromodel iteration grids.
The next window size \(w_s\) is
\(w_s = \gamma p_s\), where \(p_s\) is the pixel size
of the current iteration

OptimizationPixelsMacro

int

\(30\)

Number of pixels of the macromodel
iteration grids

MinDistMacro

float
or
None

None

Minimum ray-shooted distance
from the source \(d^s_0\) that pixels of the
first macromodel grid must satisfy
to be iterated over.
If None, all candidate pixels
are iterated over

ImprovementMacro

float
or
None

None

Improvement \(\delta\) on the ray-shooted
distance \(d^s_0\) that pixels of the
macromodel grids must satisfy to be
iterated over.
Pixels are considered for the \(n^{th}\)
iteration if their ray-shooted
distance at the \({n-1}^{th}\) iteration
is \(d^s_{n-1} < d^s_0 \cdot \delta^{n-1}\),
where \(0 < \delta \leq1\).
If None, \(d^s_{n-1}\equiv d^s_0\)
for all the iterations

OptimizationPrecisionLimitMacro

float

\(10^{-20}\rm{rad}\)

Precision of the macromodel
solutions in the source plane

OptimizationWindow

float

2

Same as OptimizationWindowMacro,
but for the complete model

OptimizationPixels

int

\(30\)

Same as OptimizationPixelsMacro,
but for the complete model

MinDist

float
or
None

None

Same as MinDistMacro, but for the
complete model

Improvement

float
or
None

None

Same as ImprovementMacro,
but for the complete model

OptimizationPrecisionLimit

float

\(10^{-20}\rm{rad}\)

Same as
OptimizationPrecisionLimitMacro,
but for the complete model

further specifications

In addition to the above-mentioned options, users can specify

Option

Type

Default

Description

OnlyMacro

bool

False

Solves for the macromodel only

NearSource

bool

False

Enables a further screening for
macroimages close to the unlensed
image position, with a customized
window size and number of pixels

SearchWindowNearSource

float
or
None

needs user input

Same as SearchWindowMacro, but for
the inspection near the unlensed
image position.
It requires \(\texttt{NearSource}= True\)

PixelsNearSource

int

\(10^3\)

Same as PixelsMacro, but for
the inspection near the unlensed
image position.
It requires \(\texttt{NearSource}= True\)

The \(\texttt{NearSource}\) option is compatible with the \(\texttt{optimization}\) mode. When both are active, the \(\texttt{optimization}\) settings applied to the screening near the unlensed image position are the same as the ones specified for the macromodel. The SearchWindowNearSource must necessarily be specified when this option is active.

Attribution and source code

Indices and tables