Kyle Miller
# import stuff
import matplotlib.pyplot as plt
import numpy as np
import pandas as pd
import seaborn as sns
import veusz
from IPython.display import HTML, IFrame, Video
import astropy.coordinates as coord
import astropy.units as u
from astropy import constants
from astropy.coordinates import SkyCoord
import ugradio
import fruitbat
import astroquery
from astroquery.gaia import Gaia
from astropy.table import QTable
# output control
import warnings
def function_that_warns():
"deprecated", DeprecationWarning)
warnings.warn(
with warnings.catch_warnings():
"ignore")
warnings.simplefilter(# this will not show a warning function_that_warns()
The first fast radio burst (FRB) was found in archival data from the Parkes observatory. The FRB occurred in 2001 but was noticed in 2007. It took a few more years before this mysterious new transient was recognized as a unique type of astrophysical event when more were reported.
This radio telescope was the main receiver for the Apollo 11 landing TV transmission, incidentally. Another interesting tidbit of Parkes history was the enigmatic phenomenon of “perytons”:
Today the resolution of this radio telescope has been eclipsed by larger single dishes and interferometers. The increase in resolution is not always helpful with FRBs, however. It has been established that FRBs are generally extragalactic in origin because there is no clustering on the galactic plane. Moreover, they appear to have no identifiable pattern and so are assumed to be isotropic much like cosmology writ large. While there has only been a moderate number of FRB detections it has been estimated that about 1000 of these events happen per day across the entire sky. This would be at all values of for the redshift, which is analogous to distance modulo the metric of spacetime.
This image was presented by Casey Law who has implemented a transient detection pipeline at the VLA. As the VLA is used by many scientists the computer back-end of the telescope can be simultaneously combed for transients. The program is called “realfast” and is an example of the emerging field of time domain astronomy.
Its a high energy burst in the radio band, typically found from 300 MHz to 8 GHz [Zhang2020].
# read in frb data and show columns
= pd.read_csv('frbcat_csv_noindex')
frb_catalog = ['backend', 'bandwidth [MHz]', 'beam [arcmin]', 'beam_rotation_angle [deg]',
units'beam_semi_major_axis [arcmin]', 'beam_semi_minor_axis [arcmin]', 'bits_per_sample [bits/sample]',
'centre_frequency [MHz]', 'channel_bandwidth [MHz]', 'circular_poln_frac',
'circular_poln_frac_error', 'circular_poln_frac_frmt',
'circular_poln_frac_frmt_err', 'count', 'dec', 'decj',
'dispersion_smearing [ms]', 'dm [pc/cm^3]', 'dm_error', 'dm_frmt', 'dm_frmt_err',
'dm_index', 'dm_index_error', 'dm_index_frmt', 'dm_index_frmt_err',
'fluence', 'fluence_error_lower', 'fluence_error_upper', 'fluence_frmt',
'fluence_frmt_err_down', 'fluence_frmt_err_up', 'flux [Jy]',
'flux_error_lower', 'flux_error_upper', 'flux_frmt',
'flux_frmt_err_down', 'flux_frmt_err_up', 'frb_name', 'gain [K/Jy]',
'galactic_electron_model', 'gb [deg]', 'gl [deg]', 'linear_poln_frac',
'linear_poln_frac_error', 'linear_poln_frac_frmt',
'linear_poln_frac_frmt_err', 'mw_dm_limit', 'notes_author_x',
'notes_author_y', 'notes_last_modified_x', 'notes_last_modified_y',
'notes_note_x', 'notes_note_y', 'npol', 'o_data_link',
'pub_description', 'pub_link', 'pub_reference', 'pub_type', 'ra', 'raj',
'receiver', 'redshift_host [z]', 'redshift_inferred', 'rm', 'rm_error',
'sampling_time [ms]', 'scattering [ms]', 'scattering_error', 'scattering_frmt',
'scattering_frmt_err', 'scattering_index', 'scattering_index_error',
'scattering_index_frmt', 'scattering_index_frmt_err',
'scattering_model', 'scattering_timescale [ms]', 'snr [S/N]', 'snr_frmt',
'spectral_index', 'spectral_index_error', 'spectral_index_frmt',
'spectral_index_frmt_err', 'telescope', 'tsys [K]', 'type', 'unnamed: 0',
'utc', 'width [ms]', 'width_error_lower', 'width_error_upper', 'width_frmt',
'width_frmt_err_down', 'width_frmt_err_up']
= units frb_catalog.columns
# list the telescopes that have reported FRBs
= pd.unique(frb_catalog['telescope'].values).shape[0]
num_telescopes 'telescope'].values) pd.unique(frb_catalog[
array(['askap', 'utmost', 'srt', 'gbt', 'apertif', 'chime', 'effelsberg',
'vla', 'dsa-10', 'fast', 'parkes', 'pushchino', 'arecibo'],
dtype=object)
# show average stats for each telescope for FRBs
'telescope').mean() frb_catalog.groupby(
13 rows x 79 columns
# plot flux vs. SNR
set(context="poster", style="whitegrid",
sns.={'grid.alpha':.5, 'grid.linestyle':'dotted', 'grid.color':'k',
rc"lines.linewidth":2.5, "axes.edgecolor":'k',
'legend.facecolor':'w', 'text.usetex':False, 'mathtext.fontset':'cm',
'legend.fontsize':'x-small', 'legend.framealpha':0.88,
'axes.formatter.use_mathtext':True, 'figure.dpi':180, 'font.family':'monospace',
"figure.facecolor":'none', 'axes.facecolor':'none'})
with warnings.catch_warnings():
"ignore")
warnings.simplefilter(=frb_catalog, x='flux [Jy]', xlim=(-10,180), palette=sns.color_palette("hls", 13),
sns.jointplot(data='snr [S/N]', ylim=(-40,420), hue='telescope', space=0, height=11, ratio=7, s=155); y
# plot sampling time vs. SNR
with warnings.catch_warnings():
"ignore")
warnings.simplefilter(=frb_catalog, x='sampling_time [ms]', xlim=(-.5,5.5),
sns.jointplot(data=sns.color_palette("hls", 13), y='snr [S/N]', ylim=(-40,420), hue='telescope',
palette=0, height=11, ratio=7, s=155); space
# plot bandwidth vs. SNR
with warnings.catch_warnings():
"ignore")
warnings.simplefilter(=frb_catalog, x='bandwidth [MHz]', xlim=(-10,420),
sns.jointplot(data=sns.color_palette("hls", 13),
palette='snr [S/N]', ylim=(-40,420), hue='telescope',
y=0, height=11, ratio=7, s=155); space
# plot histogram
# convert the utc times to pandas_datetimes
'utc'] = pd.to_datetime(frb_catalog['utc'])
frb_catalog['utc'].sort_values()[0]
frb_catalog[=(11,7), dpi=180)
plt.figure(figsize=pd.Timestamp('2001-07-24'), label="First FRB Event", alpha=.5, color='g');
plt.axvline(x=pd.Timestamp('2007-11-01'), label="First FRB Report", alpha=.5, color='b');
plt.axvline(x=pd.Timestamp('2013-07-05'), label="Second FRB Report", alpha=.5, color='r');
plt.axvline(x= plt.hist(frb_catalog['utc'].sort_values(), cumulative=True, bins=20, histtype='stepfilled', alpha=1,
bins ='gray', zorder=2.5)
color'# of FRBs reported'); plt.grid(True, linestyle='dotted', alpha=.5);
plt.ylabel(='upper left', fontsize='x-small'); plt.legend(loc
# get coordinates for all FRB observations
= []
observations for i,obs_name in enumerate(frb_catalog['frb_name']):
try:
observations.append(SkyCoord.from_name(obs_name))except:
# print(f"Unable to find {obs_name} in database. Entering manually.")
try:
'raj'][i],frb_catalog['decj'][i],unit=(u.hourangle,u.deg)))
observations.append(SkyCoord(frb_catalog[except:
print(f"Unable to enter {obs_name} manually.")
pass
pass
# get local stars from Gaia database
= '''SELECT TOP 4096 ra, dec, l, b
query_text FROM gaiadr2.gaia_source
ORDER BY random_index
'''
# retrieve Gaia data with astroquery and save to fits file
# job = Gaia.launch_job(query_text)
# gaia_data = job.get_results()
# gaia_data.write('gaia_data.fits')
= QTable.read('gaia_data.fits')
gaia_data
# make skycoord objects for stars
= []
stars for i,star in enumerate(gaia_data):
0], gaia_data[i][1])) stars.append(SkyCoord(gaia_data[i][
# get star coordinates ready for plotting
= []
stars_ra_rad = []
stars_dec_rad = []
stars_l_rad = []
stars_b_rad for star in stars:
180*u.deg).radian)
stars_ra_rad.append(star.ra.wrap_at(
stars_dec_rad.append(star.dec.radian)180*u.deg).radian)
stars_l_rad.append(star.galactic.l.wrap_at(180*u.deg).radian) stars_b_rad.append(star.galactic.b.wrap_at(
# get FRB coordinates ready for plotting
= []
ra_rad = []
dec_rad = []
galactic_lon = []
galactic_lat for obs in observations:
180 * u.deg).radian)
ra_rad.append(obs.ra.wrap_at(
dec_rad.append(obs.dec.radian)180*u.deg).radian)
galactic_lon.append(obs.galactic.l.wrap_at(180*u.deg).radian) galactic_lat.append(obs.galactic.b.wrap_at(
# make radio source skycoord objects
= ['cas A', 'crab nebula', 'M42', 'cygnus A', 'Sagittarius A*', 'Sagittarius B',
radio_sources 'M31', 'virgo A', 'SGR 1935+2154', 'FRB 121102']
= []
sources for source in radio_sources:
try:
sources.append(SkyCoord.from_name(source))except:
print(f"Unable to find {source} in database.")
pass
# get radio source coordinates ready for plotting
= []
src_ra_rad = []
src_dec_rad = []
src_galactic_lon = []
src_galactic_lat for src in sources:
180 * u.deg).radian)
src_ra_rad.append(src.ra.wrap_at(
src_dec_rad.append(src.dec.radian)180*u.deg).radian)
src_galactic_lon.append(src.galactic.l.wrap_at(180*u.deg).radian) src_galactic_lat.append(src.galactic.b.wrap_at(
# plot equatorial
# set up formatting for both plots
'font', size=10) #controls default text size
plt.rc('axes', titlesize=10) #fontsize of the title
plt.rc('axes', labelsize=10) #fontsize of the x and y labels
plt.rc('xtick', labelsize=10) #fontsize of the x tick labels
plt.rc('ytick', labelsize=10) #fontsize of the y tick labels
plt.rc('legend', fontsize=10) #fontsize of the legend
plt.rc(
=(11,7), dpi=180)
plt.figure(figsize111, projection="hammer") # options: "mollweide", "lambert", "aitoff", "hammer"
plt.subplot("FRBs (ra,dec)")
plt.title(True, linestyle='dotted')
plt.grid('o', markersize=2, alpha=.1, color='gray',label="Gaia star")
plt.plot(stars_ra_rad, stars_dec_rad,for i,s in enumerate(sources):
'o', markersize=11, alpha=.75, label=radio_sources[i])
plt.plot(src_ra_rad[i], src_dec_rad[i], '+', markersize=7, alpha=0.55, color='k', label="FRB")
plt.plot(ra_rad, dec_rad, 'o', markersize=7, alpha=0.25, color='k')
plt.plot(ra_rad, dec_rad, ='small',loc='upper right'); plt.legend(fontsize
# plot galactic
=(11,7), dpi=180)
plt.figure(figsize111, projection="hammer") # options: "mollweide", "lambert", "aitoff", "hammer"
plt.subplot("FRBs (l,b)")
plt.title(True, linestyle='dotted')
plt.grid('o', markersize=2, alpha=.1, color='gray', label="Gaia star")
plt.plot(stars_l_rad, stars_b_rad,for i,s in enumerate(sources):
'o', markersize=11, alpha=.75, label=radio_sources[i])
plt.plot(src_galactic_lon[i], src_galactic_lat[i], '+', markersize=7, alpha=0.55, color='k', label="FRB")
plt.plot(galactic_lon, galactic_lat, 'o', markersize=7, alpha=0.25, color='k')
plt.plot(galactic_lon, galactic_lat, ='small',loc='upper left'); plt.legend(fontsize
Alfven waves, which are transverse waves on magnetic field lines, have zero dispersion. In general, dispersion relations relate the wavelength/wavenumber () to its frequency . The dispersion relation form
shows how the velocity of the wave is a function of wavelength/frequency. Typically, the well-known relation is for electromagnetic waves. This can result in an emission with a bandwidth getting spread out or dispersed and requires a medium for electromagnetic waves, as the example of a prism with an index of refraction yielding a dispersion relation. The other way to represent a dispersion relation more geometrically is
where is the wavenumber and is the angular frequency . This allows the group and phase velocities to be more readily accessible:
The following image shows a single FRB and is a good visualization of dispersion, i.e., phase velocity being dependent upon frequency. The image is from https://github.com/caseyjlaw/FRB121102/blob/master/demo_FRB121102.ipynb and is showing data that was used for the first localization of a FRB with the VLA [Chatterjee2017].
# FRB chirp -- terminology
'<video width="887" height="499" controls> \
HTML( <source src="videos/FRB_chirp.mp4" type="video/mp4">\
</video>')
# Whistlers -- local electromagnetic observations
'<video width="800" height="600" controls>\
HTML( <source src="videos/whistlers.mp4" type="video/mp4">\
</video>')
Lightning like examples have been shown before in the solar system and farther away [Aleksi__2014]. This is an example from Jupiter [Imai2019].
Dispersion measure is one of the values reported for the FRBs. For pulsars this quantity is defined as the integral of free electron density over the line of sight path to the pulsar
where is the free electron density and would be an infinitesimal line element for the path , which should be defined by a metric for extragalactic observations. For FRBs the integral typically is written as [Zhang2020]
where the variability of is parameterized and cosmic expansion factored in with the redshift factor .
# show dispersion measure distribution
=(11,7), dpi=180)
plt.figure(figsize= plt.hist(frb_catalog['dm [pc/cm^3]'].astype('float'), cumulative=False, bins=25,
bins ='stepfilled', alpha=1, color='gray')
histtype'Dispersion Measure [pc/cm^3]')
plt.xlabel('# of FRBs'); plt.ylabel(
# calculation
= (constants.c*.001*u.s).to(u.km) l_frb
A FRB was localized to a known magnetar on April 28th, 2020. Repeating FRBs are much easier to localize as they can be observed with higher resolution telescopes that are too precise for initial detection. This might bias the evidence towards repeating FRBs, which would be biased towards periodic sources like pulsars. The “width” (length in time) of any given FRB is on the order of milliseconds as the units above might indicate. This fact has been used to extrapolate a length scale by multiplying this time by the speed of light to get the result [Zhang2020]:
where the length scale ends up being roughly 300 km. This is about a micro-au (!), so if this analysis is to be believed then it certainly points towards compact objects like white dwarfs, black holes, and neutron stars.
The width of FRBs appears to follow the power law behavior of scattering in cold plasma [Zhang2020]
Follow up observations at other wavelengths have imaged the magnetar in question [Bailes_2021].
This was associated with FRB200428. This is one of ~20 known magnetars within the Milky Way.
# make fruitbat objects to calculate redshift
= []
fb_objs for dm in frb_catalog['dm [pc/cm^3]']:
float(dm))) fb_objs.append(fruitbat.Frb(
# generate cosmology data
= ['WMAP5', 'WMAP7', 'WMAP9', 'Planck13', 'Planck15', 'Planck18', 'EAGLE']
cosmologies = []
redshifts for cosmo in cosmologies:
float(fb_objs[i].calc_redshift(cosmology=cosmo)) \
redshifts.append( np.array( [ for i in range(len(fb_objs)) ]) )
= pd.DataFrame(data=redshifts, index=cosmologies)
rsdf = rsdf.stack()
rsdfs = 'redshifts'
rsdfs.name =rsdfs.reset_index()
df_tidy# df_tidy['redshifts'] = np.log(df_tidy['redshifts'])
# plot redshift for each FRB for each cosmology
set(context="poster", style="whitegrid",
sns.={'grid.alpha':.5, 'grid.linestyle':'dotted', 'grid.color':'k',
rc"lines.linewidth":2.5, "axes.edgecolor":'k',
'legend.facecolor':'w', 'text.usetex':False, 'mathtext.fontset':'cm',
'legend.fontsize':'x-small', 'legend.framealpha':0.88,
'axes.formatter.use_mathtext':True, 'figure.dpi':180, 'font.family':'monospace',
"figure.facecolor":'none', 'axes.facecolor':'none'})
with warnings.catch_warnings():
"ignore")
warnings.simplefilter(= sns.jointplot(data=df_tidy, x='level_1', y='redshifts', hue='level_0', kind='kde',
splt =sns.color_palette('cubehelix', 7),
palette=0, height=11, ratio=7, s=155)
space'FRB', 'z') splt.set_axis_labels(
[Zhang2020] Zhang Bing, ``The physical mechanisms of fast radio bursts’’, Nature, vol. 587, number 7832, pp. 45-53, Nov 2020. online
[Chatterjee2017] Chatterjee C. J. and et al., ``A direct localization of a fast radio burst and its host’’, Nature, vol. 541, number 7635, pp. 58-61, Jan 2017. online
[Aleksi__2014] Aleksić J., Ansoldi S., Antonelli L. A. et al., ``Black hole lightning due to particle acceleration at subhorizon scales’’, Science, vol. 346, number 6213, pp. 1080–1084, Nov 2014. online
[Imai2019] Imai Masafumi and et al., ``Evidence for low density holes in Jupiter’s ionosphere’’, Nature Communications, vol. 10, number 1, pp. 2751, Jun 2019. online
[Bailes_2021] Bailes M, Bassa C G, Bernardi G et al., ``Multi-frequency observations of SGR J1935+2154’’, Monthly Notices of the Royal Astronomical Society, vol. , number , pp. , Mar 2021. online