A first user case for GRaLE

First set up all relevant imports:

from grale import GWEvent, LensingCalculator
from astropy.cosmology import FlatLambdaCDM
import astropy.units as u
import matplotlib.pyplot as plt
import numpy as np

GW Event properties

You first create the Gravitational Wave Event, e.g.: a binary Neutron star merger, with masses: \(m_1 = 1.5 M_\odot\) and \(m_2 = 1.2 M_\odot\), using GWEvent().

event = GWEvent(m1=1.3, m2=1.5)

Now, you can acces all the functions of GWEvent() : chirp_mass, true_redshift and redshift_range.

Lets first calculate the Chirp mass:

chirp_mass = event.chirp_mass()
print(f"Chirp Mass: {chirp_mass:.3f} M_sun")

Chirp Mass: 1.215 M_sun

You can also acces the detector frame chirp mass just by calling it (small extra, not necessary)

event.M0

1.2150360414642816

Now you can calculate the redshift of the source, if you now the source frame chirp mass

Note: If you don’t know it, you first have to calculate it using the Luminosity distance. However this package is meant to test

source frame Chirp mass ranges, if, for example the second component of the merger is unknown, as in the case of GW170817.

z_lensed = event.true_redshift(M=1.2)
print(f"Redshift of the lens: {z_lensed:.3f}")

Redshift of the lens: 0.013

If, as above said, the second component of the binary merger is unknown, we can let the masses that caused the merger vary, e.g \(m_1\) between the neutron star mass and \(m_2\) be a different star. This is a rather weird approach to this problem but we decided to go this way.

If you let the stars vary, you can use the redshift_rangefunction. You can define a mass difference delta and you need to specify the redshift of the lens.

This function prints you immediately an overview of the calculated parameters.

This function calculates the Chirp mass, the redshift range and the redshift range under a lensing scenario ( so the source has to be at a larger redshift than the lens.)

# Get redshift results over those mass ranges
z_lens = 0.0098

redshift_results = event.redshift_range(delta=0.5, step=0.01, z_lens=z_lens)

m1 ∈ [0.80, 1.80]
m2 ∈ [1.00, 2.00]
Chirp mass range: [0.78, 1.65]
Plausible redshift range: [0.0000, 1.1234]
Lensed redshift range (z ≥ 0.0098): [0.0101, 1.1234]

Lensing Calculator

Once the lensing event properties are calculates, you can continue with the properties of the Lensing Galaxy.

Important: This is done under the assumption of an Isothermal Sphere!

Now you have to set up the lens properties using the LensingCalculator, for this you need a cosmological model cosmo(e.g Flat Lambda CDM), the luminosity distance of the galaxy D_mu1, the velocity dispersion sigma, and the angular offset of the source to the lens theta_offset

# SETUP
cosmo = FlatLambdaCDM(H0=67.9, Om0=0.3065)
distance_mu1 = 40 * u.Mpc
sigma = 160
theta_offset = 10.07

lens = LensingCalculator(
    cosmo=cosmo,
    D_mu1=distance_mu1,
    sigma=sigma,
    theta_offset=theta_offset
)

In this, all needed Luminosity distances can be calculated.

You can calculate the possible magnification of the source by the lens:

d_lum = lens.luminosity_distance(z_lens)
print(f"Luminosity Distance: {d_lum:.3f} Mpc")

Luminosity Distance: 43.594 Mpc Mpc

magnification = lens.magnification(d_lum)
print(f"Magnification: {magnification:.3f}")

Magnification: 1.188

Then you can alculat the Einstein Radius of the lens:

# Calculate the angular diameter distances
# First: angular diameter distance from the lens to the source
DLS = lens.angular_diameter_distance_z1z2(z_lens, z_lensed)

# Second: angular diameter distance from the observer to the source
DS = lens.angular_diameter_distance(z_lensed)

# Calculate einstein radius
Re = lens.einstein_radius(DLS, DS)
print(f"Einstein Radius: {Re:.3f} arcsec")

Einstein Radius: 0.160 arcsec arcsec

Now you can also calculate the magnifying power that the lens has on the source:

mu_geo = lens.magnifying_power(Re)
print(f"Geometric Magnification: {mu_geo:.3f}")

Geometric Magnification: 1.016

If you want to work with the ranges as before, you can just use the compute_over_redshift_range function and everything above will be done for you:

You need to define the redshift range or extract it from the prior calculated properties:

# Extracted from above:
z_lensed = redshift_results["plausible_redshifts_lensed"] #get an array of redshifts

# Use the lensing calculator to compute the results over a range of redshifts
# This will return a dictionary with the results
lens_results = lens.compute_over_redshift_range(z_lensed, z_lens=z_lens)

z ∈ [0.0101, 1.1234]
DS range: [43.85, 1737.50] Mpc
DLS range: [1.11, 1717.17] Mpc
Magnification range of source: [1.25, 38355.84]
Einstein radius range: [0.019, 0.729] arcsec
Magnification range of lens: [1.00, 1.08]