-->

Etiquetas

El diagrama de Hubble y las supernovas

El fenómeno de la expansión del universo, que separa unos de otros todos aquellos objetos cósmicos que no están unidos por la fuerza de la gravitación, es un hecho conocido a partir de las observaciones de Hubble en 1929. El ritmo de esta expansión viene regido por la llamada ley de Hubble, que relaciona el valor de corrimiento al rojo observado en el espectro de las galaxias lejanas, con la distancia que nos separa de ellas.

La ley de Hubble suele mostrarse mediante un diagrama que representa este corrimiento al rojo, en forma de magnitud con unidades de velocidad (\(c \, z\)) en \(Km \, s^{-1}\) con respecto a las distancias medidas en Megaparsec. Claramente la mayor dificultad para elaborar un diagrama de este tipo es conocer las distancias que nos separan de objetos tan lejanos. En esta determinación, los objetos conocidos como supernovas juegan, tal como veremos, un papel fundamental. Una vez construido este diagrama se aprecia una relación aproximadamente lineal entre ambas magnitudes, siendo la constante de proporcionalidad \(H_0\) la llamada constante de Hubble. El sufijo $ _0$ indica que nos referimos aquí a su valor actual, ya que este puede variar con el tiempo en función del modelo de evolución del Universo que se adopte. El mejor valor actual de dicha constante es \(H_0 = 67.8 \; Km \; s^{-1} \; Mpc^{-1}\).

Mi propósito en este artículo es construir un diagrama de Hubble y obtener un valor aproximado de la constante de Hubble.

Importaciones y referencias

Los objetos con los que vamos a construir el diagrama de Hubble serán las supernovas de tipo Ia obtenidas en la encuesta de supernovas del Sloan Digital Sky Survey (SDSS). Puede consultarse la página principal de esta encuesta. Un resumen con enlaces a la documentación más relevante para este ejercicio puede consultarse aquí

En particular, esta es la tabla de supernovas encontradas por el SDSS en esta encuesta.

Un artículo en wikipedia sobre las supernovas y la utilidad de las supernovas de tipo Ia en el cálculo de distancias cósmicas puede consultarse aquí.

Este es un poster (en inglés) a un nivel muy accesible conn un propósito similar al de este post, realizar un diagrama de Hubble con datos de la encuesta de supernovas del SDSS.

In [1]:
%matplotlib inline

from __future__ import division

import re    # Expresiones regulares
import numpy as np
import matplotlib.pyplot as plt
import pandas as pd
import quantities as pq
import scipy.stats as s

# Generar un cuadro con versiones de las librerías utilizadas en este notebook
#https://github.com/jrjohansson/version_information
%load_ext version_information
%version_information re, numpy, matplotlib, pandas, quantities, scipy
Out[1]:
SoftwareVersion
Python2.7.9 64bit [GCC 4.4.7 20120313 (Red Hat 4.4.7-1)]
IPython2.3.1
OSLinux 3.13.0 44 generic x86_64 with debian jessie sid
re2.2.1
numpy1.9.1
matplotlib1.4.2
pandas0.15.2
quantities0.10.1
scipy0.15.1
Sat Jan 24 17:14:54 2015 CET
In [2]:
# Definición de nuevas constantes:
earth_radius = 6378160.0 * pq.m
moon_radius = 1740000.0 * pq.m
sun_radius = 695000000.0 *pq.m
earth_mass =  5.97219e24 * pq.kg
sun_mass = 332946 * earth_mass
Hubble_constant = 67.80 * (pq.km/pq.s)/(1e6*pq.pc)

# Definición de nuevas unidades
Mly = pq.UnitQuantity('megalight_year', 1e6*pq.ly, symbol = 'Mly')
Mpc = pq.UnitQuantity('megaparsec', 1e6*pq.pc, symbol = 'Mpc')

Un poco de teoria

La expansión homogenea e isótropa del universo

Una característica de la expansión del universo es que es homogenea (igual en todos los puntos) e isótropa (igual en todas las direcciones). Una expansión de estas características conserva la forma de los objetos. Por ejemplo, un triángulo se convertiría, pasado un cierto timpo, en un triángulo mayor pero semejante al primero. La razón de semejanza, llamada factor de escala dependerá del instante \(t\) considerado, y por ello se la designa por \(a(t)\).

Consideremos pues dos galaxias, y sea \(\vec r(t_0)\) el vector que los une en el instante \(t_0\), que podemos suponer por ejemplo que es el momento actual. En otro instante \(t\) el vector que une estos mismos dos objetos será:

\[\vec r(t) = a(t) \, \vec r(t_0)\]

La velocidad de alejamiento de las dos galaxias será pues (representando con un punto la derivada con respecto al tiempo):

\[\vec v_r(t) = \dot a(t) \, \vec r(t_0) = \dot a(t) \frac{\vec r(t)}{a(t)}\]

Es decir:

\[ \vec v_r = \frac{\dot a}{a} \vec r \]

Y llamando \(H(t) = \dot a(t) / a(t)\):

\[\vec v_r = H \vec r\]

O bien, eliminando la notación vectorial: \(v = H \, r\)

Es decir, la velocidad de alejamiento mutuo de las galaxias es proporcional a la distancia que las separa. El factor de proporcionalidad, \(H\) puede variar con el tiempo, pero en cada instante es la mismo en todos los puntos del universo. \(H = H(t)\) es conocido como parámetro de Hubble y como hemos mencionado no es constante puesto que puede cambiar con el tiempo. Pero lo más importante es tener en cuenta que esta relación entre velocidades y distancias es una propiedad fundamental, exacta, consecuencia matemática de una expansión homogenea e isótropa del universo, y no debe confundirse con la ley de Hubble que veremos a continuación y que tiene un carácter observacional.

Corrimiento cósmico hacia el rojo

Cuando observamos la luz emitida en cualquier longitud de onda por un objeto alejado, se constata una diferencia entre la longitud de onda medida por el observador y la de la luz emitida por el objeto. Esta discrepancia se conoce como corrimiento al rojo y se caracteriza por la cantidad sin unidades \(z\) definida como:

\[z = \frac{\lambda_{\text{obs}}-\lambda_{\text{emit}}}{\lambda_{\text{emit}}} = \frac{\Delta \lambda}{\lambda_{emit}}\]

Que también se expresa como:

\[ z + 1 = \frac{\lambda_{obs}}{\lambda_{emit}} \]

Ahora solo voy a considerar el efecto de corrimiento al rojo que es consecuencia de la expansión del universo. Si imaginamos este como una cuadrícula en expansión, estaríamos en el caso en que tanto la fuente emisora como el observador ocupan nodos de esta cuadrícula. Es decir, suponemos que no existen movimientos peculiares propios de cada objeto (en realidad siempre existirán, pero dadas las grandes distancias que se consideran, su efecto será despreciable).

Pues bien, en este caso puede demostrarse que la longitud de onda instantanea \(\lambda(t)\) es proporcional al factor de escala, y se puede escribir:

\[\lambda = \text{cte} \cdot a \Rightarrow \lambda \propto a\]

Es decir, al expandirse el espacio, la longitud de onda se expande en la misma proporción, y midiendo la longitud de onda recibida y comparandola con la emitida, es posible saber cuanto se ha expandido el universo mientras que la luz ha estado viajando hasta llegar a nosotros. Es como si la expansión del universo "estirara" las ondas de luz que viajan a su través.

Por ejemplo, si en el instante \(t_{emit}\) se emite un pulso luminoso y nosotros lo recibimos en el instante \(t_{ahora}\), podremos escribir:

\[1 + z = \frac{\lambda_{ahora}}{\lambda_{emit}} = \frac{a(t_{ahora})}{a(t_{emit})}\]

A este corrimiento al rojo, que no es una manifestación del efecto Doppler, sino que es consecuencia exclusivamente de la expansión del universo, y está, como se ha visto, íntimamente relacionado con el factor de escala de la expansión, se le llama corrimiento al rojo cosmológico, y es al que nos referiremos en lo sucesivo en este post.

Ley de Hubble

Como ya hemos comentado, en 1929 Hubble formuló su famosa ley que afirma la proporcionalidad entre la cantidad de corrimiento al rojo medida en el espectro de la luz que observamos proveniente de las galaxias, y su distancia obtenida por diversos métodos, aunque principalmente mediante la utilización del módulo de distancia resultante de la comparación entre la magnitud absoluta y la magnitud aparente del objeto:

\[ z = \frac{H_0}{c}d \]

En su momento, Hubble interpretó el corrimiento al rojo de las galaxias como una manifestación del efecto Doppler como consecuencia de su alejamiento de nosotros. Como por otro lado, para valores de z pequeños el efecto Doppler admite la aproximación \(v \approx cz\), siendo \(v\) la velocidad de alejamiento de la fuente, la ley de Hubble se escribe a menudo como \(v = H_0 d\), si bien, para dejar claro el método de obtención de esta velocidad, es preferible escribir la ley de Hubble como:

\[ c \, z = H_0 \, d \]

Sobre esto cabe hacer una serie de precisiones:

  • A diferencia de la ley teórica \(v = H(t) \, r\), la ley de Hubble \(c \, z = H_0 \, d\) es puramente experimental, siendo el resultado de múltiples observaciones. De hecho, esta relación lineal deja de cumplirse para altos valores de z

  • La distancia \(r\) que aparece en la relación teórica entre velocidades de recesión y distancias, es la distancia propia, exacta, que nos separa del objeto, medida en un sistema de referencia inercial del observador. Es decir, sería la distancia que mediriamos poniendo una a continuación de otra reglas que no se expanden con el universo.

  • Por el contrario, \(d\) (observese el cambio de notación) es la distancia estimada por diferentes métodos de observación. Lo más frecuente es que d sea la "distancia de luminosidad" \(d_L\) calculada mediante el módulo de distancia que por ejemplo se vio en este post. La distancia de luminosidad proporciona una buena aproximación a la distancia propia \(r\) para objetos "próximos", pero se desvia de ella en el caso de los objetos más lejanos.

  • La ley \(v = H(t) \, r\) relaciona velocidades y distancias en el momento "actual", es decir, son velocidades y distancias ahora, en el momento de hacer la observación, medidas en el tiempo propio del observador. En cambio, \(z\) y \(d\) están referidos al momento en que la luz fue emitida por el objeto observado.

  • Debido a la coincidencia de la velocidad \(v\) con la magnitud \(c \, z\) para valores pequeños de \(z\), podemos escribir: \(H(t_{actual}) = H_0\). La constante de Hubble es el valor que toma, en el instante actual, el parámetro de Hubble.

Una excelente explicación de todas estas sutilezas en relación con la ley de Hubble se puede encontrar en las notas del curso de astronomía extragaláctica de Mark Whittle.

Tras esta introducción teórica paso a obtener el diagrama de Hubble correspondiente a las supernovas de tipo Ia de la encuesta de supernovas del SDSS.

Cálculo de la distancia a una supernova

Tal como se ha comentado, centraré mi atención en las supernovas de tipo Ia. Sus principales características con vistas a su utilización en la elaboración de un diagrama de Hubble son las siguientes:

  • Su magnitud absoluta máxima en los filtros azul (B) y visual (V) en el sistema fotomérico UBV de Johnson es:

    \[M_U \approx M_V = -19.3 \pm 0.3\]

  • La magnitud aparente \(m\) se determina examinando su curva de luz que es una curva que representa la evolución de su magnitud \(m\) en función del tiempo.

  • Al ser (aproximadamente) fija su magnitud absoluta \(M\), se dice que las supernovas de tipo Ia son "candelas estándar".

  • Una vez disponemos del dato de magnitud relativa \(m\), podemos calcular la distancia utilizando el módulo de distancia, que expresando la distancia \(d_{Mpc}\) en Megaparsec se formula como:

    \[\mu = m - M = 5 \, \log_{10}d_{Mpc} + 25\]

In [3]:
# Ejemplo: supongamos que en la curva de luz de la supernova Ia
# medimos una magnitud aparente con el filtro V de 19.5
m = 19.5
M = -19.3
k = (m-M)/5 - 5
d = 10**k
print "La distancia de luminosidad es de %d Mpc" %d
La distancia de luminosidad es de 575 Mpc

Para calcular la magnitud máxima pueden aplicarse correcciones adicionales, siendo probablemente la más importante la Phillips relationship, que no vamos a utilizar en este post.

Lectura de las curvas de luz

Para progresar en nuestro estudio necesitaremos las curvas de luz de las supernovas de tipo Ia identificadas por la encuesta del SDSS. Estas han sido publicadas en el siguiente artículo: "The Sloan Digital Sky Survey-II: Photometry and Supernova Ia Light Curves from the 2005 data" de Holtzman et al 2008, AJ, 136, 2306.

Realmente lo que nos van a interesar más son las propias curvas de luz en formato electrónico, las cuales pueden descargarse de la propia página del SDSS Supernova Survey

Así pues, lo primero que debemos hacer es descargar todos estos ficheros en un subdirectorio de nuestro directorio de trabajo. En mi caso este subdirectorio va a llamarse "SDSS_light_curves". Lo que sigue a continuación asume que ya se han descargado las curvas de luz en dicho directorio.

In [4]:
# Lista con nombres de todos los ficheros en el subdirectorio
lc_files = !ls ../SDSS_light_curves
lc_files = list(lc_files.n.split('\n'))
lc_files = lc_files[1:]    # Eliminamos de la lista el fichero ReadMe
In [5]:
#Una muestra de los nombres
print lc_files[0:3]
['SN10028.snphot-5.04.sum', 'SN10096.snphot-5.04.sum', 'SN10106.snphot-5.04.sum']

Abriendo uno cualquiera de estos ficheros con un editor de texto se observa que las 18 primeras filas contienen una descripción de los diferentes campos de los datos. La línea 19 es una línea de cabecera con los nombres de las columnas de datos, y luego hay una serie de líneas con datos separados por espacios en blanco. Se va a leer a continuación el primero de los ficheros creando un dataframe de Pandas con los campos que nos interesan. Una vez creado el dataframe con los datos del primer fichero, recorreré todos los ficheros en el subdirectorio agregando a nuestro dataframe los datos de la "light curve" de cada supernova.

In [6]:
filename = lc_files[0]
df = pd.read_table('../SDSS_light_curves/'+filename, 
     skiprows=18, sep=r'\s+') #Separación por espacios
df.columns
Out[6]:
Index([u'#FLAG', u'MJD', u'FILT', u'MAG', u'MERR', u'MSKYERR', u'MGALERR', u'FLUX', u'FLUXERR', u'SKYERR', u'GALERR', u'NPRE', u'TELE', u'RUN', u'STRIP'], dtype='object')
In [7]:
# Utilizaré solamente las cuatro primeras columnas
df = df.ix[:,0:4]

# Cambiando el nombre de la primera columna
df = df.rename(columns={'#FLAG':'FLAG'})
df.head()
Out[7]:
FLAG MJD FILT MAG
0 0 53616.363706 1 25.659
1 32 53616.360389 2 22.739
2 0 53616.361218 3 27.020
3 32 53616.362877 4 21.498
4 4 53616.362048 0 25.049

A continuación se extraerá el valor del redshift "Z" que se encuentra en la segunda línea del fichero:

In [8]:
with open('../SDSS_light_curves/'+filename) as f:
    f.readline()
    scndline = f.readline()
resp = re.search(r'[\d.]+', scndline)  # Usamos una expresión regular
Z = resp.group()
print scndline
print type(Z), Z
# Z:  0.06533 ZERR:  0.00016 SDSS TYPE:  120

<type 'str'> 0.06533

Y añadimos a nuestro dataframe una columna con el valor de Z, convirtiendolo a número flotante

In [9]:
df['Z'] = float(Z)

También tendremos que incluir en el dataframe el id de la supernova, el cual puede extraerse del nombre del fichero:

In [10]:
print filename
SN10028.snphot-5.04.sum

In [11]:
# El id de la supernova lo vamos a obtener del nombre del fichero
# mediante una expresión regular
resp = re.search(r'\w+', filename)
SN = resp.group()
print SN
SN10028

A continuación, se va a añadir una columna al dataframe con el id de la supernova

In [12]:
df['SN'] = SN
In [13]:
df.head() # Una muestra de los datos
Out[13]:
FLAG MJD FILT MAG Z SN
0 0 53616.363706 1 25.659 0.06533 SN10028
1 32 53616.360389 2 22.739 0.06533 SN10028
2 0 53616.361218 3 27.020 0.06533 SN10028
3 32 53616.362877 4 21.498 0.06533 SN10028
4 4 53616.362048 0 25.049 0.06533 SN10028

Este es el modelo de dataframe que queremos. Ahora vamos a recorrer los nombres de ficheros restantes agregando filas.

In [14]:
for filename in lc_files[1:]:
    dftemp = pd.read_table('../SDSS_light_curves/'+filename, 
                           skiprows=18, sep=r'\s+')
    dftemp = dftemp.ix[:,0:4]
    dftemp = dftemp.rename(columns={'#FLAG':'FLAG'})
    with open('../SDSS_light_curves/'+filename) as f:
        f.readline()
        scndline = f.readline()
    resp = re.search(r'[\d.]+', scndline)
    Z = resp.group()
    dftemp['Z'] = float(Z)
    resp = re.search(r'\w+', filename)
    SN = resp.group()
    dftemp['SN'] = SN
    df = pd.concat([df, dftemp], ignore_index=True)
In [15]:
df.info()
<class 'pandas.core.frame.DataFrame'>
Int64Index: 16708 entries, 0 to 16707
Data columns (total 6 columns):
FLAG    16708 non-null int64
MJD     16708 non-null float64
FILT    16708 non-null int64
MAG     16708 non-null float64
Z       16708 non-null float64
SN      16708 non-null object
dtypes: float64(3), int64(2), object(1)
memory usage: 913.7+ KB

In [16]:
df[100:110] # Muestra de los datos
Out[16]:
FLAG MJD FILT MAG Z SN
100 32 53704.265479 1 18.955 0.06533 SN10028
101 1184 53704.262162 2 18.825 0.06533 SN10028
102 32 53704.262991 3 19.032 0.06533 SN10028
103 160 53704.264650 4 19.148 0.06533 SN10028
104 32 53704.263820 0 20.735 0.06533 SN10028
105 8 53616.395290 1 25.260 0.07775 SN10096
106 0 53616.391972 2 26.469 0.07775 SN10096
107 32 53616.392802 3 23.493 0.07775 SN10096
108 32 53616.394460 4 21.202 0.07775 SN10096
109 128 53616.393631 0 26.155 0.07775 SN10096

Tal como se indica en la cabecera de los ficheros de datos de curvas de luz, las medidas con un flag > 1024 no son fiables. Vamos a limpiar nuestro dataframe quitando los valores con un flag \(\geq\) 1024

In [17]:
b = df['FLAG'] < 1024
df = df[b]
df.info()
<class 'pandas.core.frame.DataFrame'>
Int64Index: 15838 entries, 0 to 16707
Data columns (total 6 columns):
FLAG    15838 non-null int64
MJD     15838 non-null float64
FILT    15838 non-null int64
MAG     15838 non-null float64
Z       15838 non-null float64
SN      15838 non-null object
dtypes: float64(3), int64(2), object(1)
memory usage: 866.1+ KB

Representación de la curva de luz

Aunque no es esencial para el propósito de representar un diagrama de Hubble, voy generar una representación gráfica de la curva de luz de una de las supernovas cuyos datos se han leido, lo que ayudará a comprender mejor todo el proceso.

In [18]:
# Creación de una lista con todos los identificativos de Supernovas
list_SN = list(df['SN'].unique())

# Se elige uno de ellos
SNid = list_SN[30]

# Seleccionar solo valores del filtro 1 = 'g'
b = (df['SN'] == SNid) & (df['FILT']==1)
df_g = df[b].sort(columns = 'MJD')
# Y a continuación solo valores del filtro 2 = 'r'
b = (df['SN'] == SNid) & (df['FILT']==2)
df_r = df[b].sort(columns = 'MJD')

Ahora se va a representar gráficamente la curva de luz correspondiente a los filtros "g" y "r"

In [19]:
f, (ax1, ax2) = plt.subplots(1,2, sharex=True)
f.set_size_inches(12, 5)
f.suptitle(SNid, fontsize=20)
f.subplots_adjust(hspace=0)
ax1.plot(df_g['MJD'], df_g['MAG'],'o-g')
ax1.invert_yaxis()
ax1.set_title('Filtro g', fontsize=15)
ax1.set_xlabel('MJD')
ax1.set_ylabel('MAG')
ax2.plot(df_r['MJD'], df_r['MAG'],'o-r')
ax2.invert_yaxis()
ax2.set_title('Filtro r', fontsize=15)
ax2.set_xlabel('MJD')
ax2.set_ylabel('MAG');

En esta gráfica se aprecia muy bien el aumento de brillo aparente de la supernova hasta alcanzar un máximo a partir del cual decrece más lentamente. Es en este máximo cuando su magnitud absoluta alcanza un valor conocido que es el que nos permitirá estimar la distancia a que se encuentra de nosotros.

Las fechas están en forma MJD (Modified Julian Date). Las unidades son días, por lo que podemos también hacernos idea de la duración del estallido de una supernova. El aumento brusco del brillo aparente seguido de su caida menos abrupta es claramente apreciable en un periodo de solo un mes. El periodo de tiempo total registrado es de 90 días.

Obtención de las magnitudes aparentes mínimas

Una vez en posesión del dataframe con los datos de todas las curvas de luz, vamos a construir el diagrama de Hubble para estas supernovas. Una determinación rigurosa de las distancias requeriría comparar las curvas de luz disponibles con un conjunto de curvas "patrón" y aplicar correcciones para determinar la magnitud absoluta máxima. Me limitaré a crear un dataframe seleccionando la observación del mayor brillo aparente (es decir, valor mínimo de la magnitud aparente correspondiente) en los filtros "g" y "r", el valor de Z y el identificativo de la supernova.

In [20]:
b = df['FILT'] == 1    # Seleccionar el filtro g
dftemp1 = df.ix[b,['MAG', 'Z', 'SN']]
dftemp1.columns = ['g', 'Z', 'SN']
grouped = dftemp1.groupby('SN') # Agrupar por id de supernova
df_g = grouped.min() #Seleccionar valores mínimos

b = df['FILT'] == 2    # Seleccionar filtro r
dftemp2 = df.ix[b,['MAG', 'SN']]
dftemp2.columns = ['r', 'SN']
grouped = dftemp2.groupby('SN') # Agrupar por id de supernova
df_r = grouped.min()  #Seleccionar valores mínimos de la magnitud

# Finalmente se concatenan las columnas
# de los dos dataframes anteriores
# Pandas se encargará de alinear los valores
# correspondientes a una misma supernova
df_gr = pd.concat([df_g, df_r], axis=1)
df_gr['g'] = df_gr['g'].astype(np.float)
df_gr['r'] = df_gr['r'].astype(np.float)
df_gr = df_gr[['Z', 'g', 'r']] #Reordenar las columnas
df_gr.head()    # ver una muestra
Out[20]:
Z g r
SN
SN10028 0.06533 18.292 18.354
SN10096 0.07775 20.323 19.904
SN10106 0.14720 20.844 20.665
SN1032 0.12975 20.410 20.205
SN10434 0.10430 19.432 19.369

Corrección de la extinción

El fenómeno de la extinción interestelar es el obscurecimiento de la luz que nos llega proveniente de los objetos celestes, debido a su absorción y dispersión más o menos intensa por las nubes de polvo cósmico existentes en el medio interestelar. En consecuencia los objetos celestes nos parecen menos brillantes, y por lo tanto más lejanos, de lo que correspondería a la distancia en que están situados realmente.

Por lo tanto, debemos aplicar una corrección a los valores de las magnitudes aparentes en cada uno de los filtros, restando en cada caso una constante \(A_{\lambda} > 0\) cuyo efecto será obtener una magnitud menor, es decir, correspondiente a un objeto más brillante, que sería la magnitud aparente que observaríamos si no existiera la extinción. A esta nueva magnitud le podremos aplicar la fórmula del módulo de distancia para obtener una mejor aproximación de la distancia que lo separa de nosotros.

Veamos un ejemplo: la supernova SN722 presenta su menor magnitud aparente en la banda del filtro \(g\) con un valor de g = 19.394. La corrección de extinción para su galaxia anfitriona en las longitudes de onda de este filtro es \(A_g\) = 0.10. Estimar su distancia, sin y con la aplicación de la corrección.

In [21]:
df_gr.ix['SN722',:]
Out[21]:
Z     0.08658
g    19.39400
r    19.19400
Name: SN722, dtype: float64
In [22]:
M = -19.3
m = 19.394
k = (m-M)/5 - 5
d = 10**k    # esta será la distancia en Mpc
print "distancia sin corrección: %.2f Mpc" %d
A = 0.10
m -= A
k = (m-M)/5 - 5
d = 10**k
print "distancia con corrección: %.2f Mpc" %d
distancia sin corrección: 548.02 Mpc
distancia con corrección: 523.36 Mpc

Podemos estar seguros de que estas correcciones no han sido ya aplicadas, porque en la cabecera de cada fichero de curvas de luz figura el comentario "NO extinction correction applied to mags or fluxes"

El problema que tenemos que resolver ahora es obtener las correcciones de extinción para cada supernova y banda. Estos datos figuran en la cabecera de cada fichero, en la línea 11 (contando desde 1) hay una línea del tipo:

# SFD extinctions [ugriz] = [ 0.13 0.10 0.07 0.05 0.04 ]

De modo que vamos a leer cada fichero y obtendremos los valores de la corrección para cada supernova

In [23]:
lc_files = !ls ../SDSS_light_curves
lc_files = list(lc_files.n.split('\n')) # Crear lista con nombres de ficheros
lc_files = lc_files[1:]    # Eliminar de la lista el fichero ReadMe
In [24]:
#Crear lista con los nombres de las supernovas
# Obtenidos de la lista de nombres de los ficheros
SN = lambda filename: re.search(r'\w+', filename).group()

#Crear sendas listas con los valores del index y las columnas
# del dataframe a crear con las extinciones 
# para cada SN y cada filtro
index = [SN(filename) for filename in lc_files]
columns = ['Au', 'Ag', 'Ar', 'Ai', 'Az']

#Crear un dataframe vacio pero con la estructura que necesitamos
df_ext = pd.DataFrame(index=index, columns=columns)
df_ext.head()
Out[24]:
Au Ag Ar Ai Az
SN10028 NaN NaN NaN NaN NaN
SN10096 NaN NaN NaN NaN NaN
SN10106 NaN NaN NaN NaN NaN
SN1032 NaN NaN NaN NaN NaN
SN10434 NaN NaN NaN NaN NaN
In [25]:
# Finalmente se recorren todos los ficheros
# Extrayendo los valores de las extinciones
# y poblando el dataframe df_ext

for filename in lc_files:
    id = SN(filename)  # Extraer id de la supernova
    with open('../SDSS_light_curves/'+filename) as f:
        for linea in f:
            if 'extinctions' in linea: # Esta es la linea con las extinciones
                l_campos = linea.split()
                
                # Escribir extinciones en el dataframe
                df_ext.ix[id] = [float(a) for a in l_campos[6:11]]
                break
In [26]:
# Comprobamos que todo parece correcto
df_ext.head()
Out[26]:
Au Ag Ar Ai Az
SN10028 0.13 0.1 0.07 0.05 0.04
SN10096 0.14 0.1 0.07 0.06 0.04
SN10106 0.38 0.28 0.2 0.15 0.11
SN1032 0.38 0.28 0.2 0.15 0.11
SN10434 0.48 0.36 0.26 0.2 0.14

Una vez se dispone de todas las correcciones de extinción, se van a corregir las magnitudes aparentes de los filtros \(g\) y \(r\) en el dataframe "df_gr" que era el que contenía los valores mínimos de la magnitud aparente en cada uno de estos dos filtros para cada supernova. La corrección consistirá simplemente en restar de cada magnitud aparente el valor de la corrección para ese filtro y supernova.

In [27]:
# Hacer una copia "profunda" del dataframe original
# Es decir, con copia de los datos
df_gr_ext = df_gr.copy()

# Aplicar las correcciones de extinción
df_gr_ext['g'] = df_gr['g'] - df_ext['Ag']
df_gr_ext['r'] = df_gr['r'] - df_ext['Ar']
In [28]:
# Se puede comprobar en una muestra
# que los valores obtenidos son los correctos
df_gr_ext.head()
Out[28]:
Z g r
SN
SN10028 0.06533 18.192 18.284
SN10096 0.07775 20.223 19.834
SN10106 0.14720 20.564 20.465
SN1032 0.12975 20.13 20.005
SN10434 0.10430 19.072 19.109

Corrección debida al sistema fotométrico utilizado

Hay aún otra corrección que voy a llevar a cabo. Es debida a que se van a utilizar los valores de magnitud absoluta máxima M = -19.3 que es un valor válido para los filtros B o V del sistema fotométrico UBV, pero aquí estamos utilizando el sistema de filtros del SDSS que es diferente.

Una manera de tener en cuenta este hecho es hacer una conversión entre las magnitudes del sistema "ugriz" del SDSS y las magnitudes en el sistema UBV. Existen diversos algoritmos, que pueden consultarse en esta página. Voy a utilizar por ejemplo la transformación Lupton (2005), aunque no tengo claro si esta es la transformación más adecuada en este caso. Por otro lado, llevar a cabo esta conversión presenta la ventaja adicional de que en el cálculo de por ejemplo la magnitud aparente en el filtro "V" intervienen las magnitudes aparentes de los filtros "g" y "r" del SDSS, con lo cual estaremos tomando una especie de "promedio" entre ambos filtros.

He aquí la fórmula que se va a utilizar:

\[V = g - 0.5784 (g - r) - 0.0038\]

Ahora ampliaré el dataframe df_gr_ext anterior, que ya llevaba las correcciones de extinción incorporadas, creando una nueva columna calculada con la fórmula anterior:

In [29]:
df_gr_ext['V'] = df_gr_ext['g'] - 0.5784 * (df_gr_ext['g']-df_gr_ext['r']) - 0.0038
df_gr_ext.head()
Out[29]:
Z g r V
SN
SN10028 0.06533 18.192 18.284 18.24141
SN10096 0.07775 20.223 19.834 19.9942
SN10106 0.14720 20.564 20.465 20.50294
SN1032 0.12975 20.13 20.005 20.0539
SN10434 0.10430 19.072 19.109 19.0896

Obtención del diagrama de Hubble

Determinación de las velocidades de recesión

Examinemos en primer lugar el rango de valores de \(z\) de las supernovas en la encuesta del SDSS

In [30]:
df_gr_ext['Z'].min(), df_gr_ext['Z'].max()
Out[30]:
(0.01306, 0.42199999999999999)

Esto nos dice que existe al menos una observación con un valor de z = 0.42, es decir, si utilizamos en este caso la aproximación \(v = z \, c\) tendríamos un objeto que se desplazaría a un 42% de la velocidad de la luz, lo cual constituye ya una velocidad relativista. Esto significa que el uso de esta aproximación no es válido en casos como este si el resultado lo interpretamos como la velocidad real de recesión de la supernova con respecto de nosotros.

No obstante, tal como hemos visto, la constante de Hubble, para valores elevados del parámetro z se sigue definiendo como:

\[H_0 = \frac{c \, z}{d} \]

Utilizaré por tanto, como es habitual, \(c \, z\) en la realización del diagrama de Hubble, si bien esta magnitud no debe interpretarse en general como la velocidad real de recesión de la supernova, sino tan solo como una aproximación a esta velocidad para \(z\) pequeño. La relación exacta entre la velocidad de recesión debida a la expansión del universo, y la cantidad \(z\) de corrimiento al rojo, depende del modelo de universo que se considere, y requiere el uso de la Relatividad General. Al respecto puede consultarse el articulo Expanding Confusion: common misconceptions of cosmological horizons and the superluminal expansion of the Universe.

In [31]:
z = df_gr_ext['Z']

#Series de Pandas sin unidades, con velocidades en Km/s
v = z * pq.c.rescale(pq.km/pq.s).magnitude

Cálculo de las distancias

Para las distancias vamos a utilizar la expresión del módulo de distancia cuando las distancias \(d\) se miden en Mpc:

\[\mu = m - M = 5 \, \log_{10}d + 25\]

Tomando como magnitud absoluta de cualquier supernova de tipo Ia: \(M = -19.3\) y siendo \(m\) el valor de la magnitud aparente máxima (el menor valor) tomado del dataframe:

In [32]:
M = -19.3
#m = df_gr_ext['V']
m = df_gr_ext['V'].astype(np.float)
k = (m-M)/5 - 5

# d es una Series de Pandas con la distancia en Mpc
d = 10**k

Obtención de \(H_0\) y el diagrama de Hubble

Ahora ajustaré por mínimos cuadrados una recta a los valores de distancias y velocidades

In [33]:
slope, intercept, r_value, p_value, std_err = s.linregress(d,v)
print 'El coeficiente de determinación es r^2 = %.2f' %r_value**2
print 'H0 = %.2f' %slope
El coeficiente de determinación es r^2 = 0.86
H0 = 60.71

In [34]:
f, ax = plt.subplots(figsize=(12,8))

ax.scatter(d, v, edgecolors='none')
ax.set_xlim(0, 2500)
ax.set_ylim(0,140000)


ax.set_title('Diagrama de Hubble encuesta supernovas tipo Ia SDSS-II',
             fontsize=20)
ax.set_xlabel('Distancias en Mpc')
ax.set_ylabel('$c \, z \; Km \, s^{-1}$')
ax.xaxis.label.set_fontsize(15)
ax.yaxis.label.set_fontsize(20)

ax.plot([0, 3000], [intercept, slope * 3000 + intercept], '-k', lw=2)
ax.text(1000, 80000,'Valor obtenido $H0 = %.2f \; Km \; s^{-1} \; Mpc^{-1}$'
        %slope, fontsize=15, rotation = 39)
ax.plot([0, 3000], [intercept, 67.80 * 3000 + intercept], ':r', lw=2)
ax.text(500, 130000,'Mejor valor actual $H0 = 67.80 \; Km \; s^{-1} \; Mpc^{-1}$',
        fontsize=15, rotation = 41, color='r')
ax.grid();

Como puede verse, se ha obtenido un valor menor que el mejor valor actual en aproximadamente un 10%, lo cual considero que puede ser una buena aproximación dado el método tan básico que se ha empleado en la obtención de la magnitud aparente correspondiente al momento de máxima radiación de las supernovas estudiadas. A este respecto puede ser interesante leer la azarosa historia de la determinación de \(H_0\) en The Ups and Downs of the Hubble Constant. ¡Ha costado unos 75 años obtener un valor de \(H_0\) con una incertidumbre menor del 10%!

Bien, pues aquí concluye este artículo. ¡Hasta la próxima!