TFGRuizReynésDaniel

Vista previa del material en texto

Facultat de Ciències
Memòria del Treball de Fi de Grau
Estudio numérico de la fusión de un sistema 
binario de agujeros negros
Daniel Ruiz Reynés
Grau de Física
Any acadèmic 2012-13
DNI de l’alumne: 43213705V
Treball tutelat per Sascha Husa
Departament de Física
S'autoritza la Universitat a incloure el meu treball en el Repositori Institucional per a la 
seva consulta en accés obert i difusió en línea, amb finalitats exclusivament acadèmiques 
i d'investigació
 
Paraules clau del treball: 
Agujero negro, Sistema binario, Ondas gravitacionales, Horizonte de sucesos, MareNostrum III.
✓
Índice general
1. Sistema binario de agujeros negros 5
1.1. Agujeros negros . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 5
1.2. Horizonte de sucesos . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 6
1.3. Ondas gravitacionales . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 8
1.4. Sistema binario . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 10
1.5. Sistema binario como fuente de ondas gravitatorias . . . . . . . . . . . . . . . . . . . . . . 10
2. Formalismo 3+1 12
2.1. La separación 3+1 del espacio-tiempo . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 12
2.2. Curvatura extŕınseca . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 14
2.3. Proyección de las ecuaciones de Einstein . . . . . . . . . . . . . . . . . . . . . . . . . . . . 15
2.4. La formulación BSSN . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 16
3. Estudio numérico de un sistema binario de agujeros negros 18
3.1. Ecuación de onda en 1+1 dimensiones . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 18
3.1.1. Datos iniciales . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 19
3.1.2. Método numérico . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 19
3.1.3. Convergencia . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 20
3.2. Sistema binario en 3+1 dimensiones . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 20
3.2.1. Datos iniciales . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 21
3.2.2. Evolución del sistema . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 22
3.2.3. Emisión gravitatoria . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 22
3.2.4. Método numérico . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 23
3.3. Simulaciones con MareNostrum III . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 24
3.3.1. Parámetros y datos iniciales. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 24
3.3.2. Resultados de las simulaciones. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 25
Introducción
La Relatividad General es la teoŕıa formulada por Albert Einstein que describe la gravedad. En esta
teoŕıa la fuerza gravitatoria se entiende como una deformación del espacio-tiempo, esta interpretación
geométrica del espacio-tiempo como una variedad 4-dimensional nos permite entender el universo de for-
ma totalmente diferente a como se hab́ıa entendido anteriormente.
Por muy abstracta que sea la relatividad, ésta ha podido confirmarse en varias ocasiones. La curva-
tura de las trayectorias de luz debido al campo gravitatorio de un cuerpo muy masivo[1], la precesión
del perihelio de Mercurio [2], la pérdida de enerǵıa de sistemas binarios por emisión gravitatoria [3] y el
correcto funcionamiento de los satélites utilizando las leyes de la relatividad general, son algunas de las
evidencias que confirman la teoŕıa.
De la teoŕıa de Einstein surgen asombrosas predicciones como el Big Bang, los agujeros negros, o las
ondas gravitacionales.
Las ondas gravitacionales son deformaciones del espacio-tiempo producidas por cualquier cuerpo masi-
vo acelerado respecto un observador en reposo, o por sistemas no esféricos acelerados. La amplitud de
dichas ondas es muy débil, lo que dificulta su detección. Para poder detectarlas se necesitan fenómenos
de dimensiones astrof́ısicas como la colisión de dos agujeros negros o la explosión de una supernova. Más
adelante nos centraremos en estos aspectos.
La primera evidencia observacional de la existencia de las ondas gravitacionales se debe a Hulse y
Taylor [3]. Descubrieron que el sistema binario PSR B1913+16, formado por un pulsar y una estrella
de neutrones en órbita, sufŕıa una disminución del periodo. Esta pérdida de enerǵıa concuerda con las
pérdidas por emisión gravitatoria predichas por la relatividad general. Desde entonces ha habido algunas
otras evidencias con sistemas mas relativistas [4]. Aún aśı, todav́ıa no se han detectado directamente.
Durante los últimos años se han llevado a cabo grandes proyectos, resultado de la colaboración entre va-
rias universidades, con el objetivo de detectarlas. Estos proyectos han desembocado en grandes detectores
como LIGO[5], VIRGO [6], GEO600 [7], que podŕıan abrir una nueva forma de observar el universo.
Es aqúı donde entra el estudio de los agujeros negros, pues es clave entender como se comporta este
sistema para saber que debemos buscar.
Los primeros en abordar el problema de los agujeros negros de forma numérica fueron Hahn y Lindquist[8]
en 1964, pero no fue hasta mediados de los años setenta cuando Smarr [9] y Eppley [10] realizaron las
primeras simulaciones de la colisión de dos agujeros negros. Hay que tener en cuenta que la capacidad
computacional en estos años era baja y se impońıan simetŕıas para simplificar el problema.
En las dos ultimas décadas, con el aumento de las capacidades computacionales y la construcción del
LIGO, ha habido importantes avances en el estudio de la fusión de agujeros negros. Entre 2005 y 2006
hubo tres avances importantes [11],[12] y [13]. BAM [14], LEAN[15], SpEC [16] son ejemplos de códigos
desarrollados para estudiar sistemas binarios y su producción de ondas gravitacionales.
En este proyecto se pretende estudiar los conceptos de agujero negro y onda gravitacional según la
relatividad general, más concretamente la emisión de ondas gravitacionales producidas por un sistema
binario de agujeros negros, para ello se utilizaran códigos de ordenador con el objetivo de resolver las
ecuaciones de Einstein y analizar los resultados producidos por el sistema binario.
Caṕıtulo 1
Sistema binario de agujeros negros
En este primer caṕıtulo vamos a tratar los aspectos más importantes de un sistema binario de agujeros
negros. El proceso de formación, como describir estos objetos, son algunos de los aspectos que trataremos
de forma descriptiva, como también la emisión gravitacional producida por un sistema binario y las fases
de la fusión.
1.1. Agujeros negros
Se cree que los agujeros negros se forman cuando en una estrella muy masiva se termina el combustible
nuclear y disminuye la presión interna. Esto genera un desequilibrio haciendo colapsar la estrella. Podemos
hacernos un idea de este proceso sin más que estudiar la ecuación Newtoniana de equilibrio hidrostático,
ρ(r)
d2�r
dt2
= −
dP
dr
r̂ −
GM(r)ρ(r)
r2
r̂. (1.1)
Es fácil darse cuenta que si el gradiente de presión no contrarresta el término gravitatorio la aceleración
es no nula y se produce la contracción. Una vez se produce el colapso toda la masa queda concentrada
en una región muy pequeña produciendo un campo gravitatorio muy intenso.
Dependiendo de la masa inicial de la estrella ésta puede evolucionar de diferente forma, una de las
posibilidades es que si la masa es mayor que unas 30 M⊙ se produzca una supernova y acabe formando un
agujero negro. Si la masa esta entre 9 y 30 M⊙ después de la supernova se forma una estrella de neutrones,
y si la masa es menor que unas 9 M⊙ el resultado es una enana blanca. Estas dos ultimas pueden acabar
colapsando a un agujeronegro si su masa excede un cierto ĺımite. Para las estrellas de neutrones este
ĺımite de masa máxima esta entre 1.5-2 M⊙, para enanas blancas es de 1.44 M⊙, éste es el ĺımite de Chan-
drasekhar (ver en [17]). Este ĺımite es consecuencia de que la presión de degeneración de neutrones y de
degeneración electrónica respectivamente no puede soportar la fuerza gravitatoria y se produce el colapso.
Un agujero negro, de acuerdo con la relatividad general, es un región del espacio donde la masa esta
concentrada hasta tal punto que la atracción gravitatoria es tan intensa que nada puede escapar, ni si-
quiera la luz. El ĺımite de esta región es una hipersuperficie llamada horizonte de sucesos y en su interior
el espacio-tiempo tiene una singularidad.
Las leyes que gobiernan la dinámica de un agujero negro son parecidas a las leyes termodinámicas, y si
se tienen en cuenta procesos cuánticos el agujero negro emite radiación de Hawking.
La geometŕıa del espacio para un agujero negro estacionario viene dada por la métrica de Kerr-
Newman en función de la masa M, el momento angular S, y la carga eléctrica Q. Los cuerpos astrof́ısicos
se pueden considerar como cuerpos eléctricamente neutros, por lo tanto la métrica se reduce a la de
Kerr. Esta geometŕıa es resultado de resolver las ecuaciones de Einstein imponiendo simetŕıa axial, que
el sistema se encuentra en el vaćıo, y que se encuentra en un estado estacionario. Esto no es siempre aśı,
6 CAPÍTULO 1. SISTEMA BINARIO DE AGUJEROS NEGROS
podemos tener un agujero negro con un disco de acreción, en este caso el agujero no estaŕıa rodeado de
vaćıo, pero hay que tener en cuenta que éste no estaŕıa en un estado estacionario. Tampoco debe preocu-
parnos el hecho de que un agujero negro no esté en un estado estacionario, donde la métrica de Kerr no
es válida, porque la duración de un proceso no estacionario es tan pequeña que resulta irrelevante, por
estas razones es válido hacer estas asunciones.
La métrica de Kerr en coordenadas de Boyer-Lindquist (t, r, θ,φ) y en un sistema natural de unidades
(donde G=1 y c=1) se escribe aśı:
ds
2 = −
�
1−
2Mr
ρ2
�
dt
2
−
4Mra sin2 θ
ρ2
dφdt+
ρ2
∆
dr
2+ρ2dθ2+
�
r
2 + a2 +
2Mra2 sin2 θ
ρ2
�
sin2 θdφ2, (1.2)
donde ρ2 = r2 + a2 cos2 θ, ∆ = r2 − 2Mr + a2, a = S/M es el parametro de Kerr.
Consideraremos el caso espećıfico en que el agujero no rota (S=0) y de carga neutra (Q=0). Esto nos
lleva a la métrica de Schwarzschild,
ds
2 = −
�
1−
2M
r
�
dt
2 +
�
1−
2M
r
�−1
dr
2 + r2dθ2 + r2 sin2 θdφ2. (1.3)
Podemos ver que para r → ∞ la métrica tiende a la métrica de Minkowski para un espacio-tiempo plano.
También se puede ver que la métrica es singular en r=0 y en r=2M y este último es el radio de Schwarzs-
child o horizonte de sucesos.
Los mejores candidatos para detectar agujeros negros son los sistemas binarios donde al agujero negro
interactúa con el otro astro. Estos sistemas también tienen interés por la emisión de ondas gravitacionales.
El sistema Cygnus X-1[18] es un caso en el que se puede inferir la presencia de una agujero negro. De
hecho, fue la primera vez que se probó su existencia. Este sistema es una binaria de rayos X formado por
un objeto compacto, un agujero negro y una estrella supergigante azul. La materia que absorbe el agujero
negro forma un disco de acreción alrededor de éste de tal forma que la materia, a una temperatura de
millones de kelvins, que está a punto de caer, emite rayos X. Mediante los datos orbitales del sistema se
puede inferir la masa del objeto compacto, es entre 7 y 15M⊙ y como la masa máxima de una estrella de
neutrones es de 2 M⊙ el objeto compacto debe ser un agujero negro.
Otro caso en el que se puede deducir la presencia de una agujero negro es en el centro de las galaxias
[19]. Algunos ejemplos de posibles galaxias con un agujero negro supermasivo en su centro son M31, M32,
M87 NGC 3115, NG C3377, NGC 4258 y NGC 4594 [20].
La mejor evidencia obtenida hasta ahora de la presencia de un agujero negro supermasivo es la obtenida
a partir del estudio del movimiento de estrellas cerca de el centro de la Via Láctea. Las observaciones
indican que estos astros orbitan alrededor de una región oscura. La presencia de un agujero negro super-
masivo es la opción que mejor concuerda con los cálculos, que sugieren la presencia de un un cuerpo de
4.3 millones de masas solares.
1.2. Horizonte de sucesos
En esta sección vamos a desarrollar un poco más detalladamente el concepto de horizonte de sucesos.
Es necesario indagar en este concepto pues el horizonte de sucesos caracteriza el agujero negro, y si
queremos realizar una simulación debemos poder localizarlo y analizarlo.
Como hemos dicho un agujero negro es una región del espacio en la que nada puede escapar, o lo que es
lo mismo que las geodésicas siempre se acaban dirigiendo hacia el interior.
Se define como horizonte de sucesos la superficie imaginaria a partir de la cual ya no hay retorno. Esta
superficie es como una membrana de tal forma que la materia puede entrar pero una vez dentro ya no
hay vuelta atrás. El horizonte de sucesos contiene mucha información geométrica sobre el espacio-tiempo
1.2 Horizonte de sucesos 7
del agujero negro por eso es tan importante.
En un sistema 2+1 dimensional (2 coordenadas espaciales y el tiempo) un sistema binario de agujeros
negros se puede visualizar utilizando el siguiente diagrama.
Figura 1.1: Evolución del horizonte de sucesos con el tiempo en una colisión de agujeros negros (Matzner
et. al [21].
La superficie representa la evolución del horizonte de sucesos con el tiempo.
Si consideramos la intersección H entre el horizonte de sucesos y una hipersuperficie a t constante que
denotaremos como Σ, en términos más sencillos una instantánea, tenemos una superficie cerrada de área
A. El teorema del área formulado por Hawking (ver en [22], [23]) nos impone que esta área A no puede
disminuir con el tiempo,
δA ≥ 0. (1.4)
Para una colisión de un sistema binario el área final debe ser mayor que la suma de las áreas iniciales.
Es interesante analizar el horizonte de sucesos a partir de las geodésicas nulas resultantes de la
geometŕıa del espacio. Para obtener dichas geodésicas podemos aplicar la siguiente condición ds2 = 0.
Nos centraremos en el caso de θ y ϕ constantes de tal forma que podamos estudiar el caso más sencillo.
Si integramos obtenemos la siguiente geodésica nula:
t = r + 2M ln |r − 2M |+ C. (1.5)
En la siguiente imagen podemos ver las diferentes trayectorias que seguiŕıa un rayo de luz visto desde
una distancia infinita donde el espacio es plano. No son más que representaciones de la ecuación de la
geodésica nula para diferentes valores de las condiciones iniciales.
0 1 2 3 4 5 60
2
4
6
8
10
12
14
r�M
t�M
Figura 1.2: Geodésicas nulas de la métrica de Schwarzschild para θ y ϕ constantes.
8 CAPÍTULO 1. SISTEMA BINARIO DE AGUJEROS NEGROS
Se puede ver que este observador percibiŕıa que el cuerpo cae hacia el agujero pero nunca llega a
pasar el horizonte. Esto puede parecer contradictorio, pero debemos tener en cuenta que el sistema de
referencia en cáıda libre es un sistema acelerado y no ocurrirá lo mismo en ambos sistemas, en el sistema
acelerado el cuerpo śı traspasa el horizonte de sucesos. Hay que fijarse en que en el interior, a partir de
r = 2M , las geodésicas se dirigen siempre hacia el interior.
El horizonte de sucesos nos proporciona información sobre las propiedades geométricas del agujero
negro, pero a la hora de realizar una simulación conseguir localizarlo durante todo el proceso es muy
dif́ıcil. Localizarlo después del proceso es interesante para estudiar las propiedades pero no es suficiente
para una simulación numérica.
Ademas la singularidad del espacio-tiempo en el interior del agujero negro debe ser excluida para no
estropear los cálculos. Una manera de abordar este problema es pensar en que la parte interior del
agujero negroestá causalmente desconectada de la parte exterior. Esto quiere decir que como ni la luz
puede escapar del interior de un agujero negro no puede haber eventos en el interior que afecten a la
región exterior. Esto nos permite eliminar la región interior del horizonte de sucesos. Eliminar esta región
nos obliga a saber en todo momento donde está el horizonte de sucesos.
1.3. Ondas gravitacionales
Las ondas gravitacionales son perturbaciones del espacio-tiempo que se propagan. Estas perturbacio-
nes son emitidas por cualquier cuerpo masivo acelerado. La existencia de las ondas gravitacionales se
deduce de la teoŕıa de Einstein. A continuación vamos a obtener las ecuaciones para las ondas gravita-
cionales utilizando la aproximación de campo débil y estudiaremos un poco su comportamiento.
La aproximación de campo débil consiste en suponer que tenemos un espacio plano descrito por la métrica
de Minkowski y una pequeña perturbación añadida. Dicho de otro modo,
gµν = ηµν + hµν , |hµν | � 1. (1.6)
Es de utilidad introducir la traza inversa de la perturbación
hµν ≡ hµν −
1
2
ηabh
α
α
(1.7)
Como tenemos libertad de elección de coordenadas podemos imponer la condición de Lorentz gauge.
∇ah
µν
= 0 (1.8)
Esto se reduce, para el caso de las ecuaciones de Einstein en el vaćıo, a la ecuación de onda
�hµν ≡ ∇α∇αhµν = 0 (1.9)
La solución para esta ecuación son ondas planas
hµν = Aµν exp (ikαx
α), (1.10)
donde el tensor Aµν es el tensor de amplitudes y kα el vector de onda. Se puede obtener, sin más que
sustituir la solución obtenida en la ecuación de onda, la siguiente condición:
kαk
α = 0. (1.11)
Y aplicando la condición de Lorentz gauge tenemos
A
αβ
kβ = 0. (1.12)
Esta última implica ortogonalidad, es decir, que la dirección de propagación es ortogonal a la amplitud
de la perturbación. Se puede obtener que la matriz de amplitudes tiene la forma,
Aµν =


0 0 0 0
0 A+ A× 0
0 A× A+ 0
0 0 0 0

 . (1.13)
1.3 Ondas gravitacionales 9
Se puede ver que tenemos dos estados de polarización, estos dos estados de polarización se pueden en-
tender mejor con la figura que tenemos a continuación.
Figura 1.3: Representación de los dos posibles estados de polarización con propagación en la dirección z
[22, pág. 313].
La polarización + es una oscilación en la dirección x y en la dirección y, en cambi,o la oscilación ×
es una oscilación a lo largo de la recta x = y y x = −y. Con esto ya tenemos como se propagan las ondas
gravitacionales y sus estados de polarización.
El gran reto que propone la existencia de ondas gravitacionales es conseguir realizar una detección de es-
tas ondas. Para ello han sido desarrollados grandes proyectos con esta intención. La idea fundamental que
esta detrás de estos detectores es medir las deformaciones relativas producidas por las ondas gravitaciona-
les para ello se utiliza un interferómetro láser. Este interferómetro consta de dos brazos perpendiculares,
un rayo de luz se separa de tal forma que tengamos dos rayos de luz coherente que viajen por cada brazo
en direcciones perpendiculares, éstos son reflejados en los extremos de cada brazo y finalmente se captan
con un detector. Cualquier deformación en uno de los brazos producirá una diferencia entre los caminos
ópticos de la luz de cada brazo produciendo un desfase en el detector.
En la siguiente figura tenemos un diagrama del interferómetro junto a las deformaciones producidas por
las ondas gravitacionales.
Figura 1.4: Diagrama de las deformaciones producidas por ondas gravitacionales en un interferómetro
láser [24].
Este sistema de detección necesita una sensibilidad muy grande para poder conseguir detectar ondas
gravitacionales. Dicha sensibilidad se consigue con diferentes factores. Uno de los más destacados es la
longitud de los brazos, pero también se necesitan sistemas de amortiguación para evitar el ruido śısmico,
sistemas de refrigeración y ópticas de alta precisión, entre otros sistemas para optimizar el rendimiento.
Estos detectores suponen una capacidad técnica muy elevada en muchos aspectos y un coste económico
imposible de asumir para una única universidad, con lo que son resultado de grandes colaboraciones. A
continuación podemos ver uno de los detectores LIGO en Livingston, Louisiana.
10 CAPÍTULO 1. SISTEMA BINARIO DE AGUJEROS NEGROS
Figura 1.5: LIGO - Laser Interferometer Gravitational Wave Observatory, en Livingston, Louisiana.
1.4. Sistema binario
Un sistema binario de agujeros negros está formado por dos agujeros negros, orbitando uno alrededor
del otro que emiten ondas gravitacionales durante el proceso de fusión y cuyo resultado final es un agujero
negro.
Normalmente se suele dividir el proceso en cuatro fases:
Newtonian Los dos agujeros negros están muy lejos el uno del otro y la emisión de ondas gravitacionales
es demasiado débil para causar la fusión. Hay que tener en cuenta otros procesos independientes del
problema de dos cuerpos como interacción de gas o fricción dinámica, para conseguir que los dos agujeros
negros se acerquen lo suficiente para producir emisión gravitatoria y causar la fusión.
Inspiral En esta parte la emisión de ondas gravitacionales es el proceso dominante. La pérdida de
enerǵıa provoca que los dos agujeros negros se acerquen cada vez más. Esta parte esta controlada por
métodos post-Newtonianos (PN)[25]. Para un sistema binario sin excentricidad y sin rotación inicial las
aproximaciones PN a mayor orden y la aproximación a un cuerpo [26] dan resultados muy parecidos hasta
que nos acercamos a la fusión.
Plunge and merger En esta fase la emisión de ondas gravitacionales se hace lo suficientemente fuerte
como para que la evolución de las órbitas no sea adiabática. Se produce la fusión y el resultado debe ser
un agujero negro. En esta fase es donde es necesario utilizar cálculo numérico para resolver las ecuaciones
de Einstein. Esto se debe a que los métodos PN fallan. Esta fase es muy corta, solo dura alrededor de
dos ciclos de la onda producida.
Ringdown En esta última fase el sistema resultante debe ser un único agujero negro, este se puede
describir como una perturbación del espacio de Kerr. En esta fase la emisión de ondas gravitacionales
puede ser expresada como una superposición de modos cuasi-normales del agujero negro final. Se puede
hacer el cálculo usando métodos perturbativos a partir de las condiciones iniciales.
1.5. Sistema binario como fuente de ondas gravitatorias
La emisión de ondas gravitacionales está estrictamente ligada al proceso de fusión del sistema binario.
Cada fase del proceso tiene una emisión de ondas caracteŕıstica, este aspecto lo trataremos más adelante
con los resultados de las simulaciones, pero es interesante hacerse una idea cualitativa de como evoluciona
la emisión a medida que se desarrolla la fusión. En la siguiente figura podemos ver de forma cualitativa
la emisión de ondas en cada fase del proceso.
1.5 Sistema binario como fuente de ondas gravitatorias 11
Figura 1.6: Representación de la emisión gravitatoria en función de la fase de fusión [27].
Otro punto de interés sobre las ondas gravitacionales es el transporte de enerǵıa. Como sabemos, las
ondas gravitacionales son muy débiles comparadas con las ondas electromagnéticas. Además, a diferen-
cia del electromagnetismo donde no hay momentos monopolares, las ondas gravitacionales no tienen ni
momentos monopolares ni dipolares, lo que implica que deben ser cuadrupolares.
Estas ondas transportan enerǵıa, y un sistema que emite ondas gravitacionales ira perdiendo enerǵıa.
Esta pérdida de enerǵıa puede ser calculada utilizando la aproximación de campo débil, para velocidades
mucho menores que la velocidad de la luz (v � c) y longitudes de onda λ menores que el tamaño del
sistema,
dE
dt
= −
G
5c5
�
i,j
�
d3Qij
dt3
�2
, (1.14)
donde Qij =
�
�(xixj − δijr2/3)d3x es el momento cuadrupolar de masa. Se puede ver como el factorG
5c5 implica una variación de enerǵıa muy pequeña, es decir, ondas débiles.
Caṕıtulo 2
Formalismo 3+1
En este caṕıtulo se pretenden explicar de forma más concreta los aspectos matemáticos necesarios
para abordar el cálculo numérico. La separación 3+1 del espacio-tiempo, la curvatura extŕınseca, las
ecuaciones de Einstein en este formalismo, y la formulación BSSN, son los aspectos que se van a estudiar.
El primer aspecto en el que debemos centrarnos son las ecuaciones de Einstein para la relatividad general.
Estas se expresan de forma covariante con la sutil expresión:
Gµν = 8πTµν . (2.1)
Esta formulación de la relatividad general no hace distinción entre el espacio y el tiempo, el tiempo es
una coordenada más. Esto tiene sus ventajas, pero en nuestro caso no nos proporciona una idea muy
intuitiva de como evoluciona el sistema, por eso vamos a hacer una distinción entre espacio y tiempo, de
tal manera que sea un sistema 3-dimensional que evoluciona con el tiempo.
2.1. La separación 3+1 del espacio-tiempo
Este formalismo está tomado del Alcubierre [28]. Para poder estudiar la evolución de nuestro sistema
binario vamos a separar el espacio del tiempo. Como nuestro espacio es lorentziano podemos tomar
nuestra variedad 4D compuesta por una familia de hipersuperficies a t constante.
Figura 2.1: Familia de hipersuperficies Σ a t constante [22, pág. 30].
En estas hipersuperficies de tres dimensiones analizaremos las cantidades geométricas en Σ inducidas
por la geometŕıa de la variedad M, en concreto la métrica y la curvatura extŕınseca.
Vamos a considerar dos hipersuperficies adyacentes Σt y Σt+dt para hacer las siguientes definiciones.
2.1 La separación 3+1 del espacio-tiempo 13
Figura 2.2: Ilustración de la función lapso y el vector desplazamiento [28, pág. 66].
La métrica 3-dimensional: Ésta es la métrica inducida en cada hipersuperficie, que nos da la noción
de distancia espacial,
dl
2 = γijdx
i
dx
j
. (2.2)
La función lapso: Ésta es la función que nos proporciona el instante de tiempo propio dτ entre dos
superficies, medido por un observador que se mueve en la dirección normal a la hipersuperficie.
dτ = α(t, xi)dt, (2.3)
donde α(t, xi) es nuestra función lapso.
El vector desplazamiento: El vector βi es la velocidad relativa entre un observador Euleriano y las
ĺıneas coordenadas espaciales. También llamado vector desplazamiento,
x
i
t+dt = x
i
t
− β
i(t, xi)dt. (2.4)
La función α y el vector βi también conocidas como funciones de gauge no son únicas, dependen de
las coordenadas elegidas.
Con esto podemos encontrar el elemento de ĺınea en términos de las funciones de gauge, aśı como la
métrica y su inversa,
ds
2 = (−α2 + βiβ
i)dt2 + 2βidtdx
i + γijdx
i
dx
j
, (2.5)
gµν =
�
−α2 + βkβk βi
βj γij
�
, g
µν =
�
−1/α2 βi/α2
βj/α2 γij − βiβj/α2
�
. (2.6)
A partir de aqúı podemos obtener el elemento de volumen 4-dimensional
√
−g = α
√
γ, (2.7)
donde g es el determinante de gµν y γ el de γij .
Es fácil obtener el vector normal a una hipersuperficie. Vendrá dado por:
nµ = (−α, 0), n
µ = (1/α,−βi/α), nµnµ = −1. (2.8)
Este vector normal nos permite obtener la métrica espacial de Σ a partir de la métrica de la variedad M
como
γµν = gµν + nµnν . (2.9)
Esto no es nada más que el operador de proyección sobre la hipersuperficie espacial.
Si consideramos t como el tiempo global asociado al conjunto de hipersuperficies se puede definir la
función lapso como:
α = (−∇t ·∇t)−1/2. (2.10)
14 CAPÍTULO 2. FORMALISMO 3+1
I el vector unitario normal a la hipersuperficie se puede expresar como:
n
µ = −α∇µt. (2.11)
Teniendo en cuenta la ecuación (2.4) podemos obtener el vector desplazamiento
β
i = −α(�n ·∇xi). (2.12)
Definiremos un 4-vector βµ = (0,βi) que es ortogonal a �n. Con esto se puede definir el vector tiempo de
la siguiente forma:
t
µ = αnµ + βµ. (2.13)
Este vector es tangente a las ĺıneas temporales, es decir, es tangente a las ĺıneas a posición constante,
y cumple que la componente en la dirección normal es tµnµ = −α lo que implica,
t
µ
∇µt = 1. (2.14)
Se puede obtener también el vector desplazamiento en términos del vector �t. No es más que la proyección
de �t sobre la hipersuperficie,
βµ = γµνt
ν
. (2.15)
A medida que avanzamos en el tiempo variando el parámetro t nos trasladamos desde la hipersuperficie
Σ0 hasta Σt y la métrica va cambiando a medida que cambiamos de hipersuperficie. La métrica espacial es
función del parámetro t, γµν(t), en definitiva la métrica es dinámica, nos marca la evolución del sistema.
2.2. Curvatura extŕınseca
La curvatura extŕınseca tiene que ver como nuestras hipersuperficies están inmersas en nuestra va-
riedad, es decir, por una parte tenemos la curvatura propia del espacio 3D (curvatura intŕınseca) que
vendrá impuesta por la métrica γij , más concretamente dada por el tensor de Riemann. Por otra parte
la curvatura extŕınseca esta relacionada con como se curva el espacio 3D con el transcurso del tiempo, o
dicho de otra forma, como vaŕıa el campo gravitatorio cada instante.
Antes de todo vamos a definir el operador de proyección sobre las hipersuperficies,
P
α
β
= δα
β
+ nαnβ . (2.16)
Bajando ı́ndices con la métrica gµν es fácil comprobar que Pαβ = γ
α
β
.
Se puede comprobar con facilidad que la proyección de cualquier vector es ortogonal al vector normal.
(Pα
β
v
β)nα = (v
α + nαnβv
β)nα = v
α
nα − nβv
β = 0. (2.17)
Ahora podemos definir la curvatura extŕınseca usando el operador de proyección de la siguiente forma:
Kµν = −P
α
µ
∇αnν = −(∇µnν + nµn
α
∇αnν). (2.18)
El tensor Kµν es de tipo espacial, esto es debido al proyector. Entonces nµKµν = nνKµν = 0. De esta
expresión ya se puede inferir que el tensor es simétrico Kµν = Kνµ.
2.3 Proyección de las ecuaciones de Einstein 15
La siguiente figura nos da una interpretación gráfica de la curvatura extŕınseca, en definitiva es una
medida de como cambia el vector normal bajo transporte paralelo.
Figura 2.3: Interpretación gráfica de la curvatura extrinseca en un hipersuperficie cualquiera [22, pág.
33].
Se puede demostrar que la curvatura extŕınseca se puede expresar en términos de la derivada de Lie
de la siguiente forma:
Kµν = −
1
2
L�nγµν . (2.19)
Usando la expresión (2.13) podemos ver lo siguiente:
Kµν = −
1
2α
Lα�nγµν = −
1
2α
(L�t − L�β)γµν . (2.20)
Teniendo en cuenta que L�t = ∂t podemos reescribir lo anterior,
∂tγij = −2αKij +Diβj +Djβi, (2.21)
donde Di es la derivada covariante espacial associada a γij , o lo que es lo mismo, la proyección espacial
de la derivada covariante total (Dµ = Pαµ ∇α).
Esta ecuación es una ecuación de evolución, nos relaciona la evolución de la métrica espacial con la cur-
vatura extŕınseca, en otras palabras, la curvatura extŕınseca en cierta manera es la curvatura debida a la
evolución temporal. Debemos fijarnos que esta expresión se deduce solo usando conceptos geométricos, no
hace referencia a ningún sistema en concreto, dicho de otra forma es una ecuación puramente cinemática.
La dinámica del sistema vendrá inducida por las ecuaciones de Einstein.
2.3. Proyección de las ecuaciones de Einstein
En este apartado vamos a escribir las ecuaciones de ligadura. Para ello es necesario estar familiarizado
con algunos conceptos como el tensor de Riemann, el tensor y el escalar de Ricci y el tensor de Einstein.
Este último se escribe en función del tensor y el escalar de Ricci como:
Gµν = Rµν −
1
2
gµνR. (2.22)
Para deducir la primera ecuación de ligadura hay que proyectar el tensor de Riemann mediante el operador
de proyección definido anteriormente, además se debe hacer uso de las ecuaciones de Gauss-Codazzi
(desarolladas en [23]). Esto puede encontrarse de forma más detallada en el Alcubierre [28] .
Como hemos dicho, proyectando y aplicando Gauss-Codazzi respectivamente se obtienen las siguientes
expresiones:
P
αµ
P
βν
Rαµβν = 2n
µ
n
ν
Gµν , P
αµ
P
βν
Rαµβν = R+K
2
−KijK
ij
. (2.23)De las que se deduce...
2nµnνGµν = R+K
2
−KijK
ij
. (2.24)
16 CAPÍTULO 2. FORMALISMO 3+1
Utilizando que la densidad de enerǵıa viene dada en términos del tensor de momento-enerǵıa como
ρ = nµnνTµν , podemos reescribir la anterior ecuación, para obtener la primera ligadura, llamada también
ligadura Hamiltoniana:
R+K2 −KijK
ij = 16πρ. (2.25)
Para obtener la segunda ligadura vamos a proyectar el tensor de Einstein,
P
αµ
n
ν
Gµν = P
αµ
n
ν
Rµν . (2.26)
Y mediante la relación Codazzi-Mainardi,
γ
αµ
n
ν
Gµν = D
α
K −DµK
αµ
. (2.27)
Se puede obtener la segunda ecuación de ligadura,
Dj(K
ij
− γ
ij
K) = 8πji, (2.28)
donde jα = −PαµnνTµν es la densidad de momento medida por un observador Euleriano. Esta ligadura
toma el nombre de ligadura de momento.
Las ligaduras (2.25) y (2.28) son ecuaciones que deben cumplirse para todo tiempo y éstas śı nos dan
información sobre el sistema f́ısico. Pueden entenderse análogamente a la divergencia de las ecuaciones
de Maxwell.
2.4. La formulación BSSN
Existen diversas formulaciones de las ecuaciones de evolución. En este trabajo no nos centraremos
en todas ni mucho menos, pero śı comentaremos las más importantes y finalmente nos centraremos en
la formulación de (Baumgarte y Shapiro [29], Shibata y Nakamura [30]) o como nosotros la llamaremos
BSSN.
La primera a la cual debemos hacer referencia es la formulación de Arnowitt, Deser y Misner [31], ADM.
Esta formulación de las ecuaciones de evolución depende de la constricción o ligadura Hamiltoniana ec.
(2.25). York reformuló estas ecuaciones para dar la formulación York-ADM que suprimı́a esta dependen-
cia con la constricción. Estrictamente las dos formulaciones son idénticas y no debeŕıa haber diferencia
entre ambas para la resolución anaĺıtica, pero la cosa cambia cuando hacemos una resolución numérica.
Las pequeñas oscilaciones de la solución pueden provocar que la constricción Hamiltoniana no se cumpla
exactamente de tal forma que la solución puede terminar por ser errónea. Estas dos formulaciones no
han resultado ser muy bien comportadas y la formulación BSSN ha ganado protagonismo debido a su
estabilidad.
El método BSSN introduce nuevas variables, primero introduciremos el factor conforme, el cual nos
permite reescalar la métrica. Se elige de forma que el determinante de la nueva métrica sea como podemos
ver a continuación:
γ̃ij = ψ
−4
γij : γ̃ = 1 ψ
4 = γ1/3. (2.29)
Donde γ y γ̃ son los determinantes respectivos de las dos métricas. En la práctica trabajaremos con el
logaritmo del factor conforme φ = lnψ = 112 lnγ.
Podemos separar la curvatura extŕınseca en su traza y su parte de traza cero como:
Aij = Kij −
1
3
γijK. (2.30)
Usando el factor conforme podemos hacer la siguiente redefinición para la curvatura extŕınseca,
Ãij = ψ
−4
Aij = e
−4φ
Aij . (2.31)
2.4 La formulación BSSN 17
También debemos definir las funciones de conexión conforme como:
Γ̃i = γ̃jkΓi
jk
= −∂j γ̃
ij
. (2.32)
Estas φ, K, γ̃ij , Ãij , Γ̃i, son las variables básicas del formalismo BSSN. A continuación vamos a escribir
las ecuaciones de evolución para estas variables,
d
dt
γ̃
ij = −2αÃij , (2.33)
d
dt
φ = −
1
6
αK, (2.34)
d
dt
Ã
ij = e−4φ
�
−DiD
i
α+ αRij + 4πα[γij(S − ρ)− 2Sij ]
�TF
+ α(KÃij − 2ÃikÃ
k
j
), (2.35)
d
dt
K = −DiDjα+ α(ÃijÃ
ij +K2/3) + 4πα(ρ+ S), (2.36)
d
dt
Γ̃i = γ̃jk∂j∂kβ
i +
1
3
γ
ij
∂j∂kβ
k
− 2Ãij∂jα+ 2α[Γ̃
i
jk
Ã
jk + 6Ãij∂jφ−
2
3
γ̃
ij
∂jK − 8πj̃
i], (2.37)
donde d
dt
= ∂t − L�β ,j̃
i = ψ4ji Di es la derivada covariante respecto la métrica espacial, Rij es el tensor
de Ricci asociado a γij y el super indicie TF indica que es la parte de traza cero del parentesis.
Con la formulación BSSN evitamos que aparezca la inestabilidad, pero además el sistema se comporta
mucho mejor que el método York-ADM y ADM. Esto fue comprobado por comparación directa de ambas
simulaciones por Baumgarte y Shapiro en la Ref. (32 del Art. Alcubierre).
Caṕıtulo 3
Estudio numérico de un sistema
binario de agujeros negros
En este caṕıtulo vamos a centrarnos en la resolución numérica de las ecuaciones para un sistema bina-
rio. Para obtener una mejor comprensión del funcionamiento vamos a empezar de manera introductoria
con la resolución de la ecuación de onda en un caso de 1+1 dimensiones, esto nos ayudara a entender
diferentes aspectos importantes en la simulación.
3.1. Ecuación de onda en 1+1 dimensiones
Vamos a empezar con la resolución de la ecuación de ondas para un sistema de una dimensión espacial
y una dimensión temporal para un espacio plano donde el elemento de ĺınea es el siguiente:
ds
2 = −dt̃2 + dx̃2. (3.1)
La ecuación de onda en términos del operador de Alembert en un sistema de unidades donde la velocidad
de la luz es c=1, se escribe
�φ = −∂
2φ
∂ t̃2
+
∂2φ
∂x̃2
= 0. (3.2)
Introduciendo el desplazamiento no nulo, de la misma forma que en (2.4), que no es más que un cambio
de coordenadas, tenemos:
dt = dt̃, dx = dx̃− βdt̃. (3.3)
Reescribiendo el elemento de ĺınea en las nuevas coordenadas se puede encontrar la métrica,
gµν =
�
−1 + β2 β
β 1
�
, g
µν =
�
−1 β
β 1− β2
�
. (3.4)
Para la resolución de la ecuación de ondas es útil reformular la ecuación en términos de derivadas de
primer orden, para ello utilizaremos la ecuación generalizada,
�φ = 1√
−g
∂µ[
√
−gg
µν
∂νφ] = 0, (3.5)
donde g es el determinante de la métrica. Desarrollando el sumatorio tenemos,
�φ =∂t[gtt∂tφ+ gtx∂xφ] + ∂x[gxt∂tφ+ gxx∂xφ] =
∂t[−∂tφ+ β∂xφ] + ∂x[β∂tφ+ (1− β
2)∂xφ] = 0.
(3.6)
Se puede comprobar que si se introducen las variables ψ = ∂xφ y π = ∂tφ − β∂xφ se puede obtener el
siguiente sistema de ecuaciones en términos de derivadas de primer orden.
∂tπ = ∂x(ψ + βπ), ∂tψ = ∂x(π + βψ) y ∂t = π + βψ. (3.7)
3.1 Ecuación de onda en 1+1 dimensiones 19
Estas ecuaciones van a permitirnos integrar en el tiempo y hallar la solución, para ello lo primero que
se necesita es el estado inicial del sistema, dicho de otro modo, las condiciones iniciales. Podemos darnos
cuenta de que hemos seguido un proceso similar al formalismo 3+1 pues hemos obtenido una ecuación
de evolución en función de la geometŕıa del espacio.
3.1.1. Datos iniciales
Las condiciones iniciales para el caso que queremos resolver pretenden ser sencillas. Para ello hemos
elegido una gaussiana, esta función nos permitirá ver de forma clara la propagación en el tiempo. A
continuación tenemos el estado inicial para todo x a tiempo cero.
π(0, x) = 0
φ(0, x) = A exp(−(x− x0)
2
/σ
2))
ψ(0, x) = −2
x− x0
σ2
φ(0, x)
(3.8)
Debemos tener en cuenta que las condiciones iniciales deben cumplir ψ = ∂xφ.
3.1.2. Método numérico
Una vez elegidas las condiciones iniciales ya podemos integrar en el tiempo. Para ello utilizaremos
el método Runge-Kutta 4 (RK4). Este algoritmo utiliza el valor de la derivada para proporcionarnos el
siguiente punto, lo que implica que debemos conocer el valor de la derivada en todo momento. Como
obtenerla lo veremos en seguida, pero antes vamos a centrarnos en el funcionamiento del algoritmo.
Definimos S(f) como la derivada de la función, tal como se ve a continuación:
∂tf = S(f). (3.9)
A partir de aqúı se calcula el punto n a partir del n− 1 de la siguiente forma:
k1 = S(t
n−1
, f
n−1),
k2 = S(t
n−1 +
∆t
2
, f
n−1 +
∆t
2
k1),
k3 = S(t
n−1 +
∆t
2
, f
n−1 +
∆t
2
k2),
k4 = S(t
n−1 +∆t, fn−1 +∆tk3),
f
n = fn−1 +
∆t
6
(k1 + k2 + k3 + k4) +O(∆t
5).
(3.10)
Hay que notar que el algoritmo debe integrar las funciones π, ψ, y φ, por lo que deberá ir alternando la
función a integrar.
Para calcular S(f) necesitamos calcular la derivada espacial de una función que nos vendrá dada (ver
ecuaciones 3.7). Para hacer esta derivada numéricamente utilizaremos diferencias finitas, en concreto en
nuestro caso de cuarto orden
h
�
i
=
hi−2 − 8hi−1 + 8hi+1 − hi+2
12∆x
, (3.11)
donde h es la función a derivar, que puede ser ψ+ βπ o π+ βψ. Obtener la derivada en todos los puntos
dela malla como hemos indicado presenta el problema que los puntos de los extremos (i=1,2,n-1,n) no
tienen vecinos y no se puede calcular la derivada. Esto se soluciona fácilmente imponiendo periodicidad,
es decir, f0 = fn, f−1 = fn−1, fn+1 = f1 y fn+2 = f2, donde n es el número total de puntos.
Una vez hecha la simulación podemos representar la solución en función del tiempo y comprobar si
los resultados son coherentes. A continuación tenemos un gráfico donde se puede ver la que la gaussiana
20
CAPÍTULO 3. ESTUDIO NUMÉRICO DE UN SISTEMA BINARIO DE AGUJEROS
NEGROS
que hab́ıamos impuesto como condición inicial se propaga en el espacio en ambas direcciones a medida
que pasa el tiempo.
Figura 3.1: Representación de φ en función del tiempo y el espacio.
3.1.3. Convergencia
Para tener la certeza de que nuestros cálculos no son incorrectos podemos hacer un análisis de con-
vergencia. Como nuestro algoritmo es de cuarto orden tenemos que la solución puede expresarse como
una suma del valor exacto más un cierto error de cuarto orden.
f = f0 + e∆t
4 +O(∆t5), (3.12)
donde f es el valor obtenido, f0 el valor exacto, e una constante de proporcionalidad y ∆t corresponde a
la discretización. Si se compara las soluciones obtenidas utilizando dos discretizaciones de ∆t/4, ∆t/2 y
∆t como se ve a continuación:
gcon =
f3 − f2
f2 − f1
=
e∆t4
�
1− 116
�
+O(∆t5)
e∆t4
�
1
16 −
1
162
�
+O(∆t5)
= 16 +O(∆t5), (3.13)
obtenemos que el cociente debe ser siempre 16, esto nos proporciona una herramienta para detectar po-
sibles fallos en el código, y hacer una estimación del error. Realizando un test de convergencia para un
tiempo intermedio obtenemos el siguiente gráfico donde podemos observar como mayoritariamente para
todos los puntos vale alrededor del valor especificado. Se puede observar como las zonas de los bordes y
donde se encuentran las gaussianas son las más problemáticas.
0 100 200 300 400 500x
0.5
1.0
1.5
2.0
Φ�x�
0 100 200 300 400 500x
5
10
15
20
gcon�x�
Figura 3.2: Representación de φ y gcon respectivamente a t=75.
3.2. Sistema binario en 3+1 dimensiones
La resolución de las ecuaciones de Einstein para la fusión de un sistema binario de agujeros negros,
debe resolverse numéricamente, pues no puede abordarse de otra forma debido a su complejidad. BAM[14]
3.2 Sistema binario en 3+1 dimensiones 21
es uno de los códigos desarrollados para resolver este problema y el que utilizamos aqúı para presentar
los resultados a los que nos dedicaremos más adelante. El código BAM utiliza el método de punción en
movimiento traducido del inglés ”moving puncture method”[14]. Este método aprovecha los grados de
libertad de las coordenadas para evitar la singularidad. Los datos iniciales poseen topoloǵıa de agujero
de gusano de Brill-Lindquist, ésta es asintóticamente plana en los extremos y compactada e identificada
por ri. Nos referiremos a las coordenadas de la singularidad ri como punción (puncture). Las punciones
nos permitirán tener localizados los agujeros negros en todo momento.
3.2.1. Datos iniciales
Los datos iniciales vienen dados por la métrica tridimensional de Riemann γij , y la curvatura Kij
inducida por la geometŕıa del espacio-tiempo en la hipersuperficie Σ, relacionados por la expresión
Kµν = −
1
2α
Lα�nγµν , (3.14)
donde �n es un vector tipo tiempo normal a la hipersuperficie y α la función lapso presentada en el caṕıtulo
2.
En el instante inicial se eligen los valores de mi, �xi, �pi y �Si, éstos determinarán la geometŕıa del espacio-
tiempo. Hay que decir que no cualquier momento inicial es válido, deben obtenerse con métodos post-
newtonianos. Una vez se tienen los datos iniciales se elige para el instante inicial un espacio plano rees-
calado mediante el factor conforme ψ,
γij = ψ
4
δij . (3.15)
Para determinar el factor conforme correcto debe obtenerse a partir de las ecuaciones de ligadura, esto nos
garantizará que cumple las ecuaciones de Einstein. Se puede obtener la parte derecha de las ecuaciones
de ligadura (2.25) y (2.28) en función del factor conforme y el lado izquierdo con los datos iniciales mi,
�xi, �pi y �Si, con esto podemos determinar ψ. Esta solución toma el nombre de solución de Bowen-York.
La solución toma la forma,
ψ = ψBL + u y ψBL = 1 +
N�
i=1
mi
2ri
, (3.16)
donde u viene determinado por una ecuación eĺıptica en R3 y es C∞ por todo excepto en la punción que
es C2. Las variables mi y ri corresponden al agujero negro i-éssimo y mi parametriza la masa de cada
agujero negro. En general la masa ADM viene dada por,
Mi = mi

1 + ui +
�
i �=j
mj
2dij

 , (3.17)
donde ui es u para cada punción, y dij la distancia entre la punción i y j. Esta expresión no tiene por
que coincidir con la masa del horizonte aparente MAH,i, pero coincide para los casos sin esṕın. Para los
casos con esṕın la masa viene dada por Christodoulou [32].
M
2
i
= M2
AH,i
+
S2
i
4M2
AH,i
, MAH,i =
�
A
16π
. (3.18)
Para completar los datos iniciales debemos añadir el valor del vector desplazamiento y la función lapso
en el instante inicial,
β
i = 0, (3.19)
α = 1 o ψ−2(t = 0), (3.20)
22
CAPÍTULO 3. ESTUDIO NUMÉRICO DE UN SISTEMA BINARIO DE AGUJEROS
NEGROS
3.2.2. Evolución del sistema
Las ecuaciones de evolución para nuestro sistema son como vimos, las obtenidas en la formulación
BSSN (2.33-2.37). En el instante inicial las variables de la formulación BSSN φ, γ̃ij , Ãij , K y Γ̃i toman
el valor:
φ = lnψ, (3.21)
Ãij = ψ
−6
Āij , (3.22)
Γ̃i = −∂j γ̃
ij
. (3.23)
Las variables γ̃ij y K no han sido cambiadas.
Las ecuaciones de evolución con las que se integran las anteriores variables son las escritas a continuación,
que se diferencian de (2.33-2.37) en los términos de densidad de enerǵıa y momento.
d
dt
γ̃
ij = −2αÃij , (3.24)
d
dt
φ = −
1
6
αK, (3.25)
d
dt
Ã
ij = e−4φ
�
−DiD
i
α+ αRij + 4πα[γij(S − ρ)− 2Sij ]
�TF
+ α(KÃij − 2ÃikÃ
k
j
), (3.26)
d
dt
K = −DiDjα+ α(ÃijÃ
ij +K2/3), (3.27)
d
dt
Γ̃i = γ̃jk∂j∂kβ
i +
1
3
γ
ij
∂j∂kβ
k
− 2Ãij∂jα+ 2α[Γ̃
i
jk
Ã
jk + 6Ãij∂jφ−
2
3
γ̃
ij
∂jK], (3.28)
donde d
dt
= ∂t − L�β ,j̃
i = ψ4ji Di es la derivada covariante respecto la métrica espacial, Rij es el tensor
de Ricci asociado a γij y el super indicie TF indica que es la parte de traza cero del paréntesis.
Se deben añadir dos condiciones para obtener una evolución estable:
1. La imposición de la ligadura det(γ) = 1 y Tr(Ãij) = 0. Éstas deben imponerse cuando se calculan
las partes de la derecha de las ecuaciones, y al final de cada paso de la evolución.
2. Cuando Γ̃i aparece sin diferenciar se usa Γ̃i = −∂j γ̃ij expĺıcitamente, si no es aśı se utiliza Γ̃i.
3.2.3. Emisión gravitatoria
Uno de los objetivos de la simulación es obtener información sobre la emisión de ondas gravitacionales.
En este apartado vamos a seguir la aproximación de Newman-Penrose, donde el escalar de Weyl Ψ4 (o
escalar de Newman-Penrose)[14] mide la radiación gravitacional transversal hacia fuera en un espacio
asintóticamente plano. No vamos a detallar aqúı su definición porque requiere de las variables del forma-
lismo York-ADM que no hemos tratado en este trabajo.
Nos será de interés a continuación tener la proyección del escalar de Newman-Penrose sobre los armónicos
esféricos. Proyectando Ψ4 sobre Y
−2
lm
tenemos la contribución de Y −2
lm
que denotamos como Alm.
Alm = �Y
−2
lm
,Ψ4� =
� 2π
0
�
π
0
Ψ4Ȳ
−2
lm
sin θdθdφ. (3.29)
Se pueden definir los armónicos esféricos en términos de las funciones de Wigner como:
Y
s
lm
(θ,φ) = (−1)s
�
2l + 1
4π
d
l
m(−s)(θ)e
imφ
, (3.30)
donde las funciones de Wigner toman la forma
d
l
ms
θ =
C2�
t=C1
(−1)t[(l +m)!(l −m)!(l + s)!(l − s)!]1/2
(l +m− t)!(l − s− t)!t!(t+ s−m)!
�
cos
θ
2
�2l+m−s−2t �
sin
θ
2
�2t+s−m
. (3.31)
Los ĺımites de suma C1 y C2 se definen como C1 = max(0,m− s) y C2 = min(l +m, l − s). Para l = 2
y s = −2 tenemos
3.2 Sistema binario en 3+1 dimensiones 23Y
−2
2−2 =
�
5
64π
(1− cos θ)2e−2iφ, Y −22−1 = −
�
5
16π
sin θ(1− cos θ)e−iφ, (3.32)
Y
−2
22 =
�
5
64π
(1 + cos θ)2e2iφ, Y −221 = −
�
5
16π
sin θ(1 + cos θ)eiφ, (3.33)
Y
−2
20 =
�
15
32π
sin2 θ. (3.34)
Hecho esto podemos centrarnos en la enerǵıa que pierde el sistema. Para ello se debe integrar el
escalar de Newman-Penrose en todas las direcciones t → ∞ desde hasta un tiempo t para obtener toda
la radiación emitida. Además hay que hacer el ĺımite r → ∞ para obtener la emisión a distancias largas.
dE
dt
= ĺım
r→∞
�
r2
16π
�
Ω
����
�
t
−∞
Ψ4dt̃
����
2
dΩ
�
(3.35)
El interés de introducir los armónicos esféricos es que podemos desarrollar la expresión anterior como
una suma de las contribuciones de cada armónico esférico. Esta es la expansión de Newman-Penrose.
dE
dt
= ĺım
r→∞


r2
16π
������
�
t
−∞
�
l,m
Almdt̃
������
2

 (3.36)
Esto permite calcular cada modo por separado. El modo l = 2, m = ±2 es el modo el cual contribuye en
mayor medida a la emisión de enerǵıa por radiación gravitatoria.
Se puede escribir dicho modo utilizando propiedades de simetŕıa A22 como:
A22 = �Y
−2
22 ,Ψ4� =
� 2π
0
�
π/2
0
(Ψ4Ȳ
−2
22 + Ψ̄4Ȳ
−2
2−2) sin θdθdφ (3.37)
A2−2 se puede obtener fácilmente cambiando m = 2 por m = −2 y viceversa, las propiedades de simetŕıa
le afectan de forma similar.
3.2.4. Método numérico
En general el funcionamiento básico del código BAM se basa en diferencias finitas para la diferencia-
ción espacial y Runge-Kutta 4 para la integración temporal, de la misma forma que en el caso de 1+1
dimensiones que vimos en la sección anterior.
Para obtener mayor eficiencia en el proceso se utiliza una versión del refinamiento de malla adaptativo
tipo Berger-Oliger [33]. Este método de refinamiento se basa en una malla cartesiana la cual tiene dife-
rentes niveles de refinamiento, como se explica a continuación.
Tenemos L niveles de refinamiento etiquetados por l = 0, ..., L− 1. Cada nivel de refinamiento consta de
una o más mallas con espaciado entre puntos hl para el nivel l. Niveles sucesivos están relacionados por
un factor de 1/2, de forma que hl = h0/2l, y hmax = h0 y hmin = hL−1.
Estas mallas con diferentes niveles de refinamiento están anidadas como se muestra en la figura.
Figura 3.3: Ejemplo de una mallas con 3 niveles de refinamiento[22, pág. 213].
24
CAPÍTULO 3. ESTUDIO NUMÉRICO DE UN SISTEMA BINARIO DE AGUJEROS
NEGROS
Además cuando dos mallas quedan superpuestas se crea una malla rectangular lo más pequeña posible
que las contenga a las dos. Debemos recordar que son mallas cúbicas con N3
l
puntos donde Nl es el número
de puntos en una dirección, solo esto ya requiere un coste computacional elevado.
Este sistema de refinamiento nos permite elegir una región de interés en la cual necesitamos más precisión.
Las punciones de los agujeros negros son regiones en las cuales es necesario hacer un refinamiento de la
malla para obtener mayor precisión. Es muy importante localizar la posición de la punción (xi
punc
) para
saber dónde hay que refinar la malla. Para ello se hace uso del método de Crank-Nicholson (ICN) utilizado
por Campanelli en [12].
∂tx
i
punc
= −βi(xj
punc
) (3.38)
Esta ecuación de advección con desplazamiento se obtiene de buscar dónde el factor conforme ψ diverge.
La amplitud de la emisión de ondas gravitatorias producidas por el sistema binario cae con la distancia
como 1/r lo que implica que se necesita un mayor refinamiento de la malla para obtener una medida pre-
cisa, pero hay que tener en cuenta que aumentar el nivel de refinamiento implica un coste computacional
mayor, y las punciones ya consumen una gran cantidad de recursos. Para obtener con precisión la emisión
gravitatoria se debe refinar la malla alĺı donde se hace la extracción de la emisión, sin embargo esto se
ve limitado por el coste computacional.
3.3. Simulaciones con MareNostrum III
En este apartado vamos a tratar los temas necesarios para la realización de las simulaciones, y sobre
todo en presentar los resultados obtenidos.
Como bien hemos comentado anteriormente la simulación de un sistema binario de agujeros negros es
un tarea computacionalmente muy costosa. Para salvar este obstáculo es necesario utilizar potentes or-
denadores e invertir mucho tiempo de espera. Las simulaciones aqúı presentadas han sido realizadas
con MareNostrum III, el superordenador más potente de España y uno de los mas potentes de Europa.
Está situado en el Barcelona Supercomputing Center (BSC).
3.3.1. Parámetros y datos iniciales.
Para empezar el proceso de simulación se deben especificar los parámetros de entrada con archivos de
parámetros, los cuales serán el input del código BAM[14]. Éste producirá los datos que posteriormente
serán analizados con Mathematica. Aqúı presentaremos 4 simulaciones que denotaremos como S0.85,
S0.5, S0.85 80 y S0.5 80. A continuación tenemos dos tablas, la primera con las condiciones iniciales de
nuestro sistema para cada simulación y la segunda con algunos de los parámetros de cada simulación.
Cuadro 3.1
Condición inicial m �rinicial (x, y, z) �pinicial (px, py, pz) �Sinicial (Sx, Sy, Sz)
S0.5 y S0.5 80 0.5 (0,5.5,0) (0.087,0,0) (0,0,0.125)
S0.85 y S0.85 80 0.5 (0,5,0) (0.091,0,0) (0,0,0.213)
Cuadro 3.2
Parámetro S0.5 S0.85 S0.5 80 S0.85 80
∆t/∆xi 0.5 0.5 0.4 0.4
∆xi
max
48 48 38.4 44.2
∆xi
min
11.7× 10−3 11.7× 10−3 18.8× 10−3 10.8× 10−3
3.3 Simulaciones con MareNostrum III 25
Una vez especificados los parámetros de entrada vamos a centrarnos en el análisis de los resultados.
Antes de todo hay que comentar que las simulaciones S0.85 y S0.5 finalizaron antes de alcanzar la fusión.
Este tipo de errores es poco frecuente, pero sucede en algunas ocasiones. Estas dos últimas simulaciones
han sido añadidas porque han sido de utilidad para calibrar las demás. Otro error se produjo durante la
simulación S0.85 80. Algunos datos del final de la simulación no fueron registrados.
3.3.2. Resultados de las simulaciones.
Vamos a centrarnos ahora en los resultados obtenidos representados en diversos gráficos, todos ellos
en el sistema natural de unidades.
A continuación tenemos la trayectoria de ambos agujeros negros para las dos simulaciones realizadas
donde se puede apreciar como las órbitas de ambos agujeros negros decaen por la emisión gravitatoria de
manera simétrica debido a que las masas son iguales. Se puede apreciar también observando la concen-
tración de ĺıneas que las emisiones de enerǵıa no son iguales.
�4 �2 0 2 4
�4
�2
0
2
4
x
y
Simulación S0.85_80.
BH2
BH1
�4 �2 0 2 4
�4
�2
0
2
4
x
y
Simulación S0.5_80.
BH2
BH1
Figura 3.4: Trayectorias de los agujeros negros.
En los siguientes gráficos podemos ver la velocidad angular del sistema para las diferentes simulaciones
en diferentes intervalos de tiempo. El primero presenta todo el proceso y el segundo se centra en la fase
inicial. Se puede apreciar como ésta va aumentando a medida que transcurre el tiempo y los dos cuerpos
están cada vez mas cerca.
0 500 1000 1500 2000
0.00
0.02
0.04
0.06
0.08
0.10
t�M�
V
el
oc
id
ad
an
gu
la
r S0.85_80
S0.5_80
S0.85
S0.5
0 10 20 30 40 50 60 70
0.014
0.016
0.018
0.020
0.022
0.024
0.026
t�M�
V
el
oc
id
ad
an
gu
la
r
S0.85_80
S0.5_80
S0.85
S0.5
Figura 3.5: Representación de la velocidad angular del sistema para las cuatro simulaciones en diferentes
intervalos de tiempo.
26
CAPÍTULO 3. ESTUDIO NUMÉRICO DE UN SISTEMA BINARIO DE AGUJEROS
NEGROS
En estos dos gráficos podemos ver la velocidad de la punción de cada agujero negro, la cual está clara-
mente relacionada con los gráficos anteriores, también se va incrementando con el transcurso del tiempo.
Es interesante fijarse en que hay pequeñas oscilaciones durante todo el proceso. Esto se debe a la excen-
tricidad y a oscilaciones en las coordenadas libres de la métrica.
0 500 1000 1500 20000.00
0.05
0.10
0.15
0.20
0.25
t �M�
Vel
oc
id
ad
pu
nc
ió
n
BH1
r: S0.85
r: S0.5
r: S0.85_80
r: S0.5_80
0 500 1000 1500 20000.00
0.05
0.10
0.15
0.20
0.25
t �M�
V
el
oc
id
ad
pu
nc
ió
n
BH2
r: S0.85
r: S0.5
r: S0.85_80
r: S0.5_80
Figura 3.6: Representación de la velocidad de la punción en función del tiempo.
A continuación tenemos representado el radio medio y el radio máximo, en la primera y segunda fila
respectivamente, para cada agujero negro. Podemos ver como al iniciarse la simulación el tamaño del
horizonte crece repentinamente, a continuación oscila para luego estabilizarse y decrecer paulatinamente
hasta llegar a la fusión. Este crecimiento espontáneo se debe a que inicialmente imponemos un vector
desplazamiento nulo y al inicio de la simulación las coordenadas tienen que adaptarse.
0 500 1000 1500 20000.0
0.2
0.4
0.6
0.8
t �M�
av
g�r AH
�
BH1
r: S0.85
r: S0.5
r: S0.85_80
r: S0.5_80
0 500 1000 1500 20000.0
0.2
0.4
0.6
0.8
t �M�
av
g�r AH
�
BH2
r: S0.85
r: S0.5
r: S0.85_80
r: S0.5_80
0 500 1000 1500 2000
0.2
0.4
0.6
0.8
t �M�
M
ax
�r AH�
BH1
max: S0.85
max: S0.5
max: S0.85_80
max: S0.5_80
r: S0.85
r: S0.5
r: S0.85_80
r: S0.5_80
0 500 1000 1500 2000
0.2
0.4
0.6
0.8
t �M�
M
ax
�r AH�
BH2
max: S0.85
max: S0.5
max: S0.85_80
max: S0.5_80
r: S0.85
r: S0.5
r: S0.85_80
r: S0.5_80
Figura 3.7: Representación del radio medio y máximo del horizonte aparente en función del tiempo.
3.3 Simulaciones con MareNostrum III 27
Los siguientes gráficos, expresados en porcentaje, nos dan una idea mas intuitiva de la deformación
de cada agujero negro y su evolución temporal. .
�
��
���
�
�����
���
��
�
�
���
�
����
�
�
����
����
�����
���
��
�
���
������
����
����
��
�
�
�
�
�
�
�
�
�
�
��
��
��
��
���
�
�
�����
�
��
�������������������������������������
��������
��
���
�
�
�
�����
����
�
�
�
��������
0 500 1000 1500 20000
5
10
15
t �M�
M
ax
�r AH�r
A
H
��1�
BH1
� r: S0.85
� r: S0.5
� r: S0.85_80
� r: S0.5_80
�
��
���
�
�����
�����
�
�
���
�
�
����
���
�
�
�����
���
�
���
���
��
��
���
�
��
�����
��
��
�
�
�
�
�
�
�
�
�
�
��
��
��
��
���
�
�
�����
�
�����������������������������������������
������
����
��
�
�
�����
�����
�
�
��������
0 500 1000 1500 20000
5
10
15
t �M�
M
ax
�r AH�r
A
H
��1�
BH2
� r: S0.85
� r: S0.5
� r: S0.85_80
� r: S0.5_80
Figura 3.8: Representación de la deformación del horizonte aparente en tanto por ciento.
En estos gráficos podemos ver diferentes magnitudes de los agujeros negros, una de ellas es la masa, que
permanece prácticamente constante durante todo el proceso. Esta última se puede comparar con el radio
medio del horizonte y con su área, también representadas.
0 10 20 30 40 50 60 70
0.0
0.5
1.0
1.5
2.0
t�M�
Simulación S0.5
Area en coord. BH2
Radio BH2
Masa BH2
Area en coord. BH1
Radio BH1
Masa BH1
0 10 20 30 40 50 60 70
0.0
0.2
0.4
0.6
0.8
t�M�
Simulación S0.85
Area en coord. BH2
Radio BH2
Masa BH2
Area en coord. BH1
Radio BH1
Masa BH1
0 500 1000 1500 2000
0
1
2
3
4
t�M�
Simulación S0.5_80
Area en coord. BH2
Radio BH2
Masa BH2
Area en coord. BH1
Radio BH1
Masa BH1
0 200 400 600 800 1000 12000.0
0.2
0.4
0.6
0.8
1.0
1.2
t�M�
Simulación S0.85_80
Area en coord. BH2
Radio BH2
Masa BH2
Area en coord. BH1
Radio BH1
Masa BH1
Figura 3.9: Representación de la masa, radio medio, y área del horizonte aparente para cada agujero
negro.
28
CAPÍTULO 3. ESTUDIO NUMÉRICO DE UN SISTEMA BINARIO DE AGUJEROS
NEGROS
A continuación tenemos la representación del momento angular de rotación para toda la evolución.
Se puede apreciar que permanece prácticamente constante durante todo el proceso y finalmente crece
debido a que el momento lineal se transforma en momento de rotación después de la fusión. Solo hay
representada una de las simulaciones, dado que la segunda no está completa.
En los otros dos gráficos podemos ver la misma representación para ambas simulaciones en intervalos de
tiempo que permiten apreciar mejor las pequeñas oscilaciones durante el proceso. Estas oscilaciones son
debidas como ya hemos mencionado a la difusividad de la simulación, un aspecto que no se ha tratado aqúı.
0 500 1000 1500 2000
0.5
0.6
0.7
0.8
0.9
t�M�M
om
en
to
an
gu
la
rd
e
ro
ta
ci
ón
Simulación S0.5_80
BH2
BH1
0 500 1000 1500 2000
0.498
0.500
0.502
0.504
t�M�M
om
en
to
an
gu
la
rd
e
ro
ta
ci
ón
Simulación S0.5_80
BH2
BH1
0 200 400 600 800 1000 1200
0.8505
0.8510
0.8515
0.8520
0.8525
0.8530
t�M�M
om
en
to
an
gu
la
rd
e
ro
ta
ci
ón
Simulación S0.85_80
BH2
BH1
Figura 3.10: Representación del momento angular de rotación en función del tiempo. La primera gráfica
muestra toda la evolución mientras que las otras dos muestran el intervalo antes de la fusión.
Este último gráfico es la representación de la emisión gravitatoria y es uno de los más importantes.
Es clave saber como será dicha emisión para que los grandes detectores puedan detectar estas ondas.
Si nos fijamos en el gráfico se detecta inmediatamente el instante de la fusión donde la mayor parte de
la radiación gravitatoria es emitida. Se puede ver como el sistema binario emite ondas gravitacionales,
las órbitas van disminuyendo y la radiación empieza a aumentar progresivamente en intensidad hasta
alcanzar la fusión, posteriormente la emisión decae hasta cero.
�1500 �1000 �500 0
�0.10
�0.05
0.00
0.05
0.10
t�M�
r�
4
Emisión gravitatoria. Simulación S0.5_80.
Figura 3.11: Representación de la emisión gravitatoria en función del tiempo.
3.3 Simulaciones con MareNostrum III 29
Se ha realizado un desplazamiento temporal de la simulación porque tener la fusión en t = 0 facilita
la comparación de las emisiones gravitatorias.
Com bien hemos comentado, la simulación S0.85 80 no esta finalizada, por lo que no se puede presentar
aqúı la emisión gravitatoria. Los siguientes gráficos corresponden a otras simulaciones con los mismos
parámetros f́ısicos que S0.5 80 y S0.85 80 cuyos resultados no fueron satisfactorios debido a un fallo en
la simulación, pero de los cuales śı se ha podido obtener la emisión gravitatoria. Estos resultados son
presentados a pesar de diferir ligeramente de los anteriores para poder hacer una comparación entre las
dos emisiones gravitatorias para diferentes datos iniciales.
�1500 �1000 �500 0
�0.10
�0.05
0.00
0.05
0.10
t�M�
r�
4
Emisión gravitatoria. Simulación S0.5_80�.
�1500 �1000 �500 0
�0.10
�0.05
0.00
0.05
0.10
t�M�
r�
4
Emisión gravitatoria. Simulación S0.85_80�.
Figura 3.12: Representación de la emisión gravitatoria en función del tiempo.
A continuación tenemos una tabla con los datos finales de la simulación S0.5 80 y otros datos de
interés.
Simulación xi,inicial tfinal Mfinal Mperdida vmáx BH1 vmáx BH2 Sfinal
S0.5 80 5.5 2230.8 0.932 0.069 0.128 0.128 0.84
Aunque estos datos son útiles para manejar las variables más cómodamente, no nos dan una idea f́ısica
de la magnitud de la fusión. Para ello vamos a obtener un par de datos en unidades más intuitivas. Estos
datos corresponden a un sistema de masa M , por lo que debemos especificar una masa total. La masa
del agujero negro para el sistema Cygnus X-1[18] era entre 7 y 15 M⊙ por lo que tomaremos un sistema
de masa total 10 M⊙. Teniendo en cuenta esto, un sistema de 10 M⊙ pierde por emisión gravitatoria
aproximadamente 0.7 M⊙. Lo más sorprendente es que el proceso completo tiene una duración de 0.2 s,
lo cual denota la violencia del proceso.
Conclusiones
En este trabajo se pretend́ıa entender dos de las predicciones más importantes de la relatividad ge-
neral: los agujeros negros y las ondas gravitacionales. Estudiando un sistema que las engloba a las dos
hemos podido acercarnos a conocer en mayor medida ambos fenómenos.
Aunque el objetivo de este trabajo es de carácter académico, hemos podido estudiar uno de los grandes
problemas de la f́ısica actual y aprender diferentes aspectos necesarios para el estudio de un problemade tales caracteŕısticas. Bien hay que decir que este problema puede ser desarrollado con mucha más
profundidad y propone posibles investigaciones en un futuro utilizando casos más complejos.
En los caṕıtulos anteriores hemos tratado conceptos f́ısicos como el concepto de onda gravitacional, el
de agujero negro, su formación y las fases que constituyen la fusión. Pero también hemos tratado nuestro
sistema desde un punto de vista más matemático para presentar la formulación BSSN de las ecuaciones
de evolución. Y no debemos olvidar la componente numérica que ha sido necesaria para la realización
de las simulaciones, tanto la utilización de potentes ordenadores, entre ellos MareNostrum III, como en
el análisis de datos realizado con Mathematica. Estos aspectos hacen que este trabajo acabe siendo un
trabajo completo que permite aprender conceptos de diferentes campos.
Respecto a los resultados de las simulaciones, ha habido dificultades en el proceso, entre ellas simu-
laciones inacabadas, pero aún aśı los resultados presentados nos han permitido entender mejor cómo se
comporta este sistema durante la fusión, cómo la trayectoria se ve afectada por la emisión, cómo vaŕıa
en consecuencia la velocidad, cómo se deforma el agujero negro, la masa perdida por el sistema, y lo más
importante cómo es la emisión gravitatoria, qué tipo de señal debe esperarse en los detectores de ondas
gravitacionales.
Un sistema binario de agujeros negros en proceso de fusión es sin duda un cataclismo. Una cantidad
enorme de enerǵıa se emite en forma de ondas gravitatorias y cruza el universo a la espera de ser detectada
a años luz de distancia, en los grandes detectores, aqúı, en la Tierra.
Bibliograf́ıa
[1] Clifford M. Will. The confrontation between general relativity and experiment. Living Reviews in
Relativity, 9(3), 2006.
[2] G. M. Clemence. The relativity effect in planetary motions. Rev. Mod. Phys., 19:361–364, Oct 1947.
[3] R.A. Hulse and J.H. Taylor. Discovery of a pulsar in a binary system. Astrophysical Journal,
195:L51–L53, Jan 1975.
[4] M. Burgay, N. D’Amico, A. Possenti, R. N. Manchester, A. G. Lyne, B. C. Joshi, M. A. McLaughlin,
M. Kramer, J. M. Sarkissian, F. Camilo, V. Kalogera, C. Kim, and D. R. Lorimer. An increased
estimate of the merger rate of double neutron stars from observations of a highly relativistic system.
Nature, 426:531–533, December 2003.
[5] Alex Abramovici, William E. Althouse, Ronald W. P. Drever, Yekta Gürsel, Seiji Kawamura, Frede-
rick J. Raab, David Shoemaker, Lisa Sievers, Robert E. Spero, Kip S. Thorne, Rochus E. Vogt, Rainer
Weiss, Stanley E. Whitcomb, and Michael E. Zucker. Ligo: The laser interferometer gravitational-
wave observatory. Science, 256(5055):325–333, 1992.
[6] B. Caron and et al. The virgo interferometer. Classical and Quantum Gravity, 14(6):1461, 1997.
[7] H. Lück and et al. The geo 600 gravitational wave detector. Classical and Quantum Gravity,
19(7):1377, 2002.
[8] Susan G. Hahn and Richard W. Lindquist. The two-body problem in geometrodynamics. Annals of
Physics, 29(2):304 – 331, 1964.
[9] L. L. Smarr. The structure of general relativity with a numerical illustration: The collision of two
black holes. PhD thesis, Texas Univ., Austin., 1975.
[10] Kenneth Eppley. Evolution of time-symmetric gravitational waves: Initial data and apparent hori-
zons. Phys. Rev. D, 16:1609–1614, Sep 1977.
[11] Frans Pretorius. Evolution of binary black-hole spacetimes. Phys. Rev. Lett., 95:121101, Sep 2005.
[12] M. Campanelli, C. O. Lousto, P. Marronetti, and Y. Zlochower. Accurate evolutions of orbiting
black-hole binaries without excision. Phys. Rev. Lett., 96:111101, Mar 2006.
[13] John G. Baker, Joan Centrella, Dae-Il Choi, Michael Koppitz, and James van Meter. Gravitational-
wave extraction from an inspiraling configuration of merging black holes. Phys. Rev. Lett., 96:111102,
Mar 2006.
[14] Bernd Brügmann, José A. González, Mark Hannam, Sascha Husa, Ulrich Sperhake, and Wolfgang
Tichy. Calibration of moving puncture simulations. Phys. Rev. D, 77:024027, Jan 2008.
[15] Ulrich Sperhake. Binary black-hole evolutions of excision and puncture data. Phys. Rev. D,
76:104015, Nov 2007.
32 BIBLIOGRAFÍA
[16] B. Szilágyi, L. Lindblom, and M. A. Scheel. Simulations of binary black hole mergers using spectral
methods. Phys. Rev. D, 80:124010, Dec 2009.
[17] B. Kiziltan. Reassessing the Fundamentals: On the Evolution, Ages and Masses of Neutron Stars.
Universal Publishers, 2011.
[18] Jeffrey.E. McClintock and Ronald.A. Remillard. Black hole binaries. 2003.
[19] Andrew King. Black holes, galaxy formation, and the mbh-σ relation. The Astrophysical Journal
Letters, 596(1):L27, 2003.
[20] J. Kormendy and D. Richstone. Inward Bound-The Search For Supermassive Black Holes In Galactic
Nuclei.
[21] Stuart L. Shapiro L. Smarr W.-M. Suen Saul A. Teukolsky Richard A. Matzner, H. E. Seidel and
J. Winicour. Geometry of a black hole collision. Science, 270:941–947, 1995.
[22] T. Baumgarte and S. Shapiro. Numerical Relativity. Solving Einstein’s Equations on the Computer.
Cambridge University Press, 2010, Chap. 7.
[23] R. M. Wald. General relativity. The University of Chicago Press, Chicago, 1984.
[24] David Reitze. Chasing gravitational waves. Nature Photonics, 2(10):582–585, 2008.
[25] L. Blanchet. Gravitational Radiation from Post-Newtonian Sources and Inspiralling Compact Bina-
ries. Living Reviews in Relativity, 5:3, April 2002.
[26] A. Buonanno and T. Damour. Effective one-body approach to general relativistic two-body dynamics.
Phys. Rev. D, 59:084006, Mar 1999.
[27] Thomas W. Baumgarte and Stuart L. Shapiro. Binary black hole mergers. Physics Today, 64(10):32–
37, 2011.
[28] Miguel Alcubierre. Introduction to 3+1 Numerical Relativity. Oxford Universty Press, 2008, Chap.
2.
[29] Thomas W. Baumgarte and Stuart L. Shapiro. Numerical integration of einstein’s field equations.
Phys. Rev. D, 59:024007, Dec 1998.
[30] Masaru Shibata and Takashi Nakamura. Evolution of three-dimensional gravitational waves: Har-
monic slicing case. Phys. Rev. D, 52:5428–5444, Nov 1995.
[31] R. Arnowitt, S. Deser, and C. W. Misner. Republication of: The dynamics of general relativity.
General Relativity and Gravitation, 40:1997–2027, September 2008.
[32] Demetrios Christodoulou. Reversible and irreversible transformations in black-hole physics. Phys.
Rev. Lett., 25:1596–1597, Nov 1970.
[33] Marsha J Berger and Joseph Oliger. Adaptive mesh refinement for hyperbolic partial differential
equations. Journal of Computational Physics, 53(3):484 – 512, 1984.
	Sistema binario de agujeros negros
	Agujeros negros
	Horizonte de sucesos
	Ondas gravitacionales
	Sistema binario
	Sistema binario como fuente de ondas gravitatorias
	Formalismo 3+1
	La separación 3+1 del espacio-tiempo
	Curvatura extrínseca
	Proyección de las ecuaciones de Einstein
	La formulación BSSN
	Estudio numérico de un sistema binario de agujeros negros
	Ecuación de onda en 1+1 dimensiones
	Datos iniciales
	Método numérico
	Convergencia
	Sistema binario en 3+1 dimensiones
	Datos iniciales
	Evolución del sistema
	Emisión gravitatoria
	Método numérico
	Simulaciones con MareNostrum III
	Parámetros y datos iniciales.
	Resultados de las simulaciones.