In [1]:
import numpy as np
import matplotlib.pyplot as plt

import cosmology as cosm

# Some plot settings
rcParams['font.size'] = 12
#rcParams['text.usetex'] = True

In [2]:
## UNITS

C      = 3.e10 # cm/s
GPC    = 1.e9 * 3.e18 # cm
ARCSEC = 5.e-6 # radians
ARCMIN = 3.e-4 # radians
YEAR   = 3.e7 # seconds

In [3]:
## Values to use an example

ZS     = 2.0
COSMO  = cosm.Cosmology()

In [4]:
print "Omega_Lambda:", COSMO.l
print "Omega_Matter:", COSMO.m

source_distance = cosm.DChi(ZS, 0.0) # Gpc
print "Distance to a z=1 source [Gpc]:", cosm.DChi(1.0, 0.0)
print "Distance to a z=" + np.str(ZS)[0:2] + " source [Gpc]:", source_distance
print "Distance to a z=4 source [Gpc]:", cosm.DChi(4.0, 0.0)


Omega_Lambda: 0.7
Omega_Matter: 0.3
Distance to a z=1 source [Gpc]: 3.08571227298
Distance to a z=2. source [Gpc]: 4.83792413272
Distance to a z=4 source [Gpc]: 6.69718960559

The basic formula and some calculations

Corrales & Paerels (2012): $$ \Delta \chi_{\rm max} \approx 4.5\ {\rm kpc}\ \chi_{Gpc} E^{-2}_{keV} a_{0.1}^{-2} $$

And in general (Corrales & Paerels 2012, Eq 1 of Miralde-Escude 1999): $$ \Delta \chi = \chi_{OS} \alpha^2 (1-x) $$

Keep in mind that the minimal x value is: $$ x_{min} \approx \frac{\alpha}{2 \sigma} $$ where $$ \sigma = \frac{10.4'}{a(0.1\ \mu{\rm m})\ E({\rm keV})} $$


In [5]:
def charsig(a, E):
    """ returns the characteristic scattering angle, in arcseconds
     a: [um] Grain size
     E: [keV] Photon energy"""
    charsig0 = 1.04 * 60 # arcsec per 1 um grain size per keV
    return charsig0 / a / E

def delta_time(alpha, x, zs=ZS):
    # alpha: [arcsec] Observation angle
    time_os = cosm.DChi(zs, 0.0) * (GPC/C) / YEAR # years
    return time_os * (alpha*ARCSEC)**2 * (1-x)

Set x=0 to get the coefficient (T):

$$ \Delta t = T (1-x) $$

In [6]:
# Time delay coefficient for 30 arcseconds
print delta_time(30, 0.0) # years


362.844309954

This agrees with my pen and paper calculation, good


In [7]:
alpha = np.logspace(-1,3,300)
plt.plot(alpha, delta_time(alpha, 0.0), 'k-')

plt.loglog()
plt.xlabel("Observation angle [alpha]")
plt.ylabel("Time delay coefficient [years]")

# Overplot the characteristic halo size
plt.axvline(charsig(0.1,1.0), ls='--', color='r')


Out[7]:
<matplotlib.lines.Line2D at 0x106b75490>

Fade time for a total halo ($\alpha = \sigma$)

The furthest you can look back along a sight line is $x_{min} \approx \alpha/2\sigma$

The largest observable halo angle is $\sim \sigma$, so the total time it takes for a halo to fade is: $$ t_{fade} \approx \frac{\chi_{OS}}{c} \frac{\sigma^2}{2} $$


In [8]:
def tfade(a=0.1, E=1.0, zs=ZS):
    time_os = cosm.DChi(zs, 0.0) * (GPC/C) / YEAR # years
    sig_rad = charsig(a, E) * ARCSEC # radiansg
    return time_os * sig_rad**2 / 2.0

In [9]:
print "Tfade for 1 keV halo from 0.1 um grains [years]:", tfade(a=0.1, E=1.0)
print "Tfade for 1 keV halo from 1.0 um grains [years]:", tfade(a=1.0, E=1.0)


Tfade for 1 keV halo from 0.1 um grains [years]: 78490.4811292
Tfade for 1 keV halo from 1.0 um grains [years]: 784.904811292

In [10]:
ztest  = np.linspace(0.5,4.0,20)
tfade_max = [tfade(a=0.1, zs=z) for z in ztest]
tfade_min = [tfade(a=1.0, zs=z) for z in ztest]

plt.plot(ztest, tfade_max, 'k-', lw=2, label='0.1 um grains')
plt.plot(ztest, tfade_min, 'k-', lw=2, alpha=0.5, label='1.0 um grains')

plt.semilogy()

plt.legend(loc='upper left', frameon=False)
plt.xlabel("Redshift of source ($z$)")
plt.ylabel("Halo lifetime [years]")


Out[10]:
<matplotlib.text.Text at 0x106b96a10>

Fade time as a function of angle

Using the fact that $x_{min} \approx \alpha/2\sigma$, I can calculate the time delay for the longest ($\chi_{max}$) value at a particular sight line ($\alpha$).

$$ \Delta t_{max} = \frac{\chi_{OS}}{c} \alpha^2 (1 - \frac{\alpha}{2 \sigma}) $$

In [11]:
def tmax(alpha, a=0.1, E=1.0, zs=ZS):
    """ Calculate the maximum time delay associated with 
        observation angle and typical grain population
        alpha: [arcmin] Observation angle
        a: [micron] Grain size
        E: [keV] Photon energy
        zs: Redshift of source"""
    time_os   = cosm.DChi(zs, 0.0) * (GPC/C) / YEAR
    alpha_rad = alpha * ARCSEC
    xterm     = 1 - alpha/(2*charsig(a,E))
    return time_os * alpha_rad**2 * xterm

In [12]:
plt.plot(alpha, tmax(alpha, a=0.1), 'k-', lw=2, label='0.1 um grains')
plt.plot(alpha, tmax(alpha, a=1.0), 'k-', lw=2, alpha=0.5, label='1.0 um grains')

plt.loglog()
plt.ylim(1.0,1.e5)
plt.xlim(1.0,1.e3)

plt.legend(loc='upper left', frameon=False)
plt.xlabel("Observation angle [arcsec]")
plt.ylabel("Lookback time [years]")

# Overplot location of characteristic scattering angle
plt.axvline(charsig(0.1,1.0), ls='--', color='k')
plt.axvline(charsig(1.0,1.0), ls='--', color='0.5')
plt.axhline(1.e3, ls=':', color='r')


Out[12]:
<matplotlib.lines.Line2D at 0x106d874d0>

Pretty plot for HEAD 2014 poster


In [13]:
## Pretty plot settings
#from astroML.plotting import setup_text_plots
#setup_text_plots(fontsize=8, usetex=True)

#import matplotlib as mpl
#mpl.rc('axes', labelsize='large')
#mpl.rc('xtick', labelsize='large')
#mpl.rc('ytick', labelsize='large')

In [14]:
fig = plt.figure(figsize=(6,2))

plt.plot(alpha, tmax(alpha, a=0.1), 'k-', lw=2, label='0.1 um grains')
plt.plot(alpha, tmax(alpha, a=1.0), 'k-', lw=2, alpha=0.5, label='1.0 um grains')

plt.loglog()
plt.ylim(1.0,1.e5)
plt.xlim(0.1,1.e3)

plt.legend(loc='upper left', frameon=False)#, fontsize=12)
plt.xlabel("Observation angle [arcsec]", size=15)
plt.ylabel("$\delta t$ [years]", size=15)

#plt.savefig('/Users/lia/Dropbox/HEAD2014/time_delay.pdf', format='pdf')


Out[14]:
<matplotlib.text.Text at 0x106e44890>

In [15]:
print "Time delay at 10'':", tmax(30.0, a=0.1, zs=2)


Time delay at 10'': 354.122090965

In [ ]: