NIVERSITA DEGLI ` STUDI DI PARMA - infn.it delle componenti del tensore di curvatura di Weyl...

UNIVERSITA DEGLI STUDI DI PARMA

FACOLTA DI SCIENZE MATEMATICHE FISICHE NATURALI

CORSO DI LAUREA SPECIALISTICA IN FISICA TEORICA

ESTRAZIONE DEL SEGNALE DI ONDE GRAVITAZIONALI DASIMULAZIONI NUMERICHE DELLE EQUAZIONI DI EINSTEIN

Relatore: Chiarissimo Prof. Roberto DE PIETRI

Candidato: Giovanni CONTI

ANNO ACCADEMICO 2010/2011

Nature uses only the longest threads to weave her patterns, so that each smallpiece of her fabric reveals the organization of the entire tapestry.

Richard P. Feynman

Il presente documento e stato redatto con LATEX, AMS- EulerPer l’analisi e l’elaborazione dati e stato utilizzato MATLAB R2008b

SOMMARIO

La seguente tesi verte sull’estrazione del segnale di onde gravitazionali dasimulazioni numeriche delle equazioni di Einstein. Il progetto e stato diviso indue parti: la prima ha previsto lo studio e l’analisi di metodi attualmenteproposti in letteratura per l’estrazione del segnale di onde gravitazionaliemesso da sorgenti autogravitanti in relativita generale, mentre la secondaha riguardato lo studio dell’implementazione su calcolatore dei metodianalizzati e la comprensione di eventi fisici e non, dovuti alla discretizzazio-ne dello spazio-tempo. Il lavoro di questa parte e consistito nella rianalisicritica dei risultati discussi in [1].

Lo studio analitico dei metodi attualmente disponibili si e basato sulformalismo di Newman-Penrose o l’identificazione delle perturbazionidello spazio-tempo di Schwarzschild con la decomposizione in multipoli.Il secondo e basato sul metodo tradizionale, sviluppato originariamenteda Regge e Weeler [4], Zerilli [5, 6], Moncrief [9] e altri per lo studioperturbativo delle oscillazioni di buchi neri non rotanti. Negli anni recenti,tuttavia, ha preso il sopravvento, nella comunita di relativita numerica,l’estrazione dei segnali dell’informazione di onde gravitazionali in terminidelle componenti del tensore di curvatura di Weyl rispetto ad una tetradenulla, usando quello che viene conosciuto come formalismo di Newman-Penrose.

Lo studio e proseguito analizzando il formalismo 3+1, mezzo che permet-te una rappresentazione piu intuitiva delle equazioni di Einstein e che aiutaa pensare ad un evoluzione dinamica del campo gravitazionale nel tempo,il che rappresenta un valido aiuto nella scrittura di codice computazionalee una piu profonda comprensione delle equazioni stesse.

Per le simulazioni svolte e stato utilizzato il sistema di calcolo Cactus/-CATIE/Whisky per la computazione in 3+1 dimensioni dell’equazioni diEinstein accoppiate ad un fluido perfetto. I test svolti al calcolatore hannopermesso un confronto diretto dei segnali delle onde permettendo unacomparazione tra i due metodi di estrazione e una valutazione sui pregi edifetti di entrambe.

iii

INDICE

1 Estrazione di onde gravitazionali 11.1 Introduzione . . . . . . . . . . . . . . . . . . . . . . . . . . . 11.2 Perturbazione della Schwarzschild . . . . . . . . . . . . . . 2

1.2.1 Espansione in multipoli . . . . . . . . . . . . . . . . . 21.2.2 Perturbazioni Pari . . . . . . . . . . . . . . . . . . . . 41.2.3 Perturbazioni Dispari . . . . . . . . . . . . . . . . . . 61.2.4 TT Gauge . . . . . . . . . . . . . . . . . . . . . . . . . 7

1.3 Gli scalari di Weyl . . . . . . . . . . . . . . . . . . . . . . . . . 81.3.1 La tetrade nulla di Newman-Penrose . . . . . . . . . 9

1.4 Energia e momento delle onde gravitazionali . . . . . . . . . 101.4.1 La media di Isaacson . . . . . . . . . . . . . . . . . . . 111.4.2 Energia e momento irradiati . . . . . . . . . . . . . . . 13

2 La sorgente delle onde 152.1 Stella di neutroni . . . . . . . . . . . . . . . . . . . . . . . . . 152.2 Soluzione di TOV . . . . . . . . . . . . . . . . . . . . . . . . . 162.3 EOS . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 172.4 Stabilita . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 18

3 Formalismo 3+1 e formulazione BSSNOK 213.1 Il formalismo 3+1 . . . . . . . . . . . . . . . . . . . . . . . . . 213.2 Equazioni ADM . . . . . . . . . . . . . . . . . . . . . . . . . . 243.3 La formulazione BSSNOK . . . . . . . . . . . . . . . . . . . . 253.4 Condizioni di Gauge . . . . . . . . . . . . . . . . . . . . . . . 27

3.4.1 Maximal slicing . . . . . . . . . . . . . . . . . . . . . . 283.4.2 Condizione di slicing iperbolico . . . . . . . . . . . . 283.4.3 Elliptic shift condition . . . . . . . . . . . . . . . . . . 303.4.4 Driver Conditions . . . . . . . . . . . . . . . . . . . . 32

4 Il Software 354.1 CACTUS e Einstein Toolkit . . . . . . . . . . . . . . . . . . . 354.2 Runge-Kutta e MOL . . . . . . . . . . . . . . . . . . . . . . . 51

v

Indice

5 Simulazioni 535.1 Introduzione . . . . . . . . . . . . . . . . . . . . . . . . . . . . 535.2 Il sistema fisico e la sua evoluzione . . . . . . . . . . . . . . . 545.3 CACTUS-CARPET-WHISKY:codice non-lineare 3D per la

relativita generale . . . . . . . . . . . . . . . . . . . . . . . . . 545.4 Gauges . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 565.5 Evoluzione della materia . . . . . . . . . . . . . . . . . . . . . 565.6 Estrazione delle onde con lo scalare di Newmann-Penrose Ψ4 575.7 Estrazione delle onde tramite multipoli . . . . . . . . . . . . 585.8 I dati iniziali . . . . . . . . . . . . . . . . . . . . . . . . . . . . 595.9 Simulazioni . . . . . . . . . . . . . . . . . . . . . . . . . . . . 60

6 Conclusioni 77

A Spin-weighted spherical harmonics 79

B Codice 81B.1 Configurazione . . . . . . . . . . . . . . . . . . . . . . . . . . 81B.2 Thornlist . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 82B.3 Parametri iniziali . . . . . . . . . . . . . . . . . . . . . . . . . 85

Bibliografia 93

vi

ELENCO DELLE FIGURE

2.1 Rappresentazione schematica del quadrato delle frequenzein funzione della densita centrale di una stella, per i modi 0,1, 2. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 19

3.1 Due adiacenti superfici spacelike. . . . . . . . . . . . . . . . . 22

4.1 Schema di funzionamento dell’Einstein Toolkit . . . . . . . . 364.2 Come appare una griglia con diversi livelli di raffinamento

nella parte centrale. . . . . . . . . . . . . . . . . . . . . . . . . 43

5.1 Alcuni test iniziali riferiti alla tabella 5.2 . . . . . . . . . . . . 605.2 Variazione di ρmax nel tempo per la tabella (5.3). . . . . . . . . . 645.3 Variazione di ρmax nel tempo per la tabella (5.4). . . . . . . . . . 645.4 Profilo densita iniziale e finale per il Test 3. In verde la configura-

zione iniziale. . . . . . . . . . . . . . . . . . . . . . . . . . . . . 655.5 Confronto densita iniziali per i test della tabella (5.3). . . . . . . . 655.6 Controllo dei constraints del momento e hamiltoniani per le

tabelle (5.3) e (5.4) . . . . . . . . . . . . . . . . . . . . . . . . . 665.7 Controllo dell’ EOS per l’ultimo step di dell’evoluzione. . . . 675.8 Confronto Ψe ottenute dalla tabella (5.3) con il metodo di Abra-

hams e Price. . . . . . . . . . . . . . . . . . . . . . . . . . . . . 675.9 Confronto Ψe ottenute dalla tabella (5.4) con il metodo di Abra-

hams e Price. . . . . . . . . . . . . . . . . . . . . . . . . . . . . 685.10 Confronti Ψe ottenute per la tabella (5.3) al variare del raggio di

estrazione. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 685.11 Confronti Ψe ottenute per la tabella (5.4) al variare del raggio di

estrazione. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 695.12 Andamento lineare dell’ampiezza del junk al variare del

raggio di estrazione. . . . . . . . . . . . . . . . . . . . . . . . 705.13 Controllo delle H := h+ − ih× ottenute con il metodo degli

scalari di Weyl. . . . . . . . . . . . . . . . . . . . . . . . . . . 705.14 Ampiezza perturbazione λ = 0.1. . . . . . . . . . . . . . . . . . 715.15 Ampiezza perturbazione λ = 0.05. . . . . . . . . . . . . . . . . . 71

vii

Elenco delle figure

5.16 Ampiezza perturbazione λ = 0.01. . . . . . . . . . . . . . . . . . 725.17 Ampiezza perturbazione λ = 0.001. . . . . . . . . . . . . . . . . 725.18 Ampiezza perturbazione λ = 0.1. . . . . . . . . . . . . . . . . . 735.19 Ampiezza perturbazione λ = 0.05. . . . . . . . . . . . . . . . . . 745.20 Ampiezza perturbazione λ = 0.01. . . . . . . . . . . . . . . . . . 745.21 Ampiezza perturbazione λ = 0.001. . . . . . . . . . . . . . . . . 75

viii

ELENCO DELLE TABELLE

5.1 Proprieta di equilibrio del modello A0. Da sinistra a destrale colonne riportano: densita centrale a riposo, densita totaledi energia centrale, massa gravitazionale, raggio, compattezza. 59

5.2 Test preliminari. . . . . . . . . . . . . . . . . . . . . . . . . . . 605.3 Test eseguiti con la massima risoluzione, l’intervallo sulla

griglia piu piccolo e 0.250 . . . . . . . . . . . . . . . . . . . . 615.4 Test eseguiti con risoluzione inferiore, l’intervallo sulla gri-

glia piu piccolo e 0.375 . . . . . . . . . . . . . . . . . . . . . . 615.5 Frequenze dominanti delle H per i vari modelli perturbati. . 73

ix

NOTAZIONI E CONVENZIONI

La segnatura dello spazio tempo utilizzata e (−1,+1,+1,+1).Indici greci (α, β, γ, δ) indicano lo spazio-tempo quadridimensiona-

le e vanno da 0 a 3, mentre indici latini scorrono tra 1 e 3. Si sfrutta laconvenzione standard di somma sugli indici ripetuti.

Viene utilizzato il sistema standard di unita c = G = M = 1.

xi

CAPITOLO 1

ESTRAZIONE DI ONDE GRAVITAZIONALI

Indice1.1 Introduzione . . . . . . . . . . . . . . . . . . . . . . . . . 11.2 Perturbazione della Schwarzschild . . . . . . . . . . . 2

1.2.1 Espansione in multipoli . . . . . . . . . . . . . . 21.2.2 Perturbazioni Pari . . . . . . . . . . . . . . . . . 41.2.3 Perturbazioni Dispari . . . . . . . . . . . . . . . 61.2.4 TT Gauge . . . . . . . . . . . . . . . . . . . . . . 7

1.3 Gli scalari di Weyl . . . . . . . . . . . . . . . . . . . . . 81.3.1 La tetrade nulla di Newman-Penrose . . . . . . 9

1.4 Energia e momento delle onde gravitazionali . . . . . 101.4.1 La media di Isaacson . . . . . . . . . . . . . . . . . 111.4.2 Energia e momento irradiati . . . . . . . . . . . . 13

La storia dell’evoluzione insegna che l’universo non hamai smesso di essere creativo o inventivo.

KARL POPPER

1.1 Introduzione

Le onde gravitazionali sono perturbazioni nella geometria dello spaziotempo che si propagano alla velocita della luce. Sono una delle piu impor-tanti predizioni della relativita generale. Sebbene non siano ancora staterilevate direttamente, ci sono forti evidenze indirette sulla loro esistenzaderivanti dallo studio di sistemi binari di pulsar come PSRB 1913-16 o ilpiu recente PSRJ 037-3039, il cui periodo orbitale cambia in accordo conquanto previsto dalla relativita generale come conseguenza di emissionidi onde gravitazionali. Questo tipo di radiazione puo portare via da unsistema isolato energia ed impulso, oltre ad informazioni sul sistema stesso.

1

1.2. Perturbazione della Schwarzschild

Ci sono due importanti approcci per l’estrazione di onde gravitazionalida simulazioni numeriche, il primo basato sulla teoria di perturbazione diuna metrica di Schwarzschild sviluppabile in multipoli pensata originaria-mente da Regge e Wheeler [4], Zerilli [5, 6] e il secondo basato sull’estra-zione in termini di componenti del tensore di curvatura di Weyl rispetto aduna tetrade nulla, conosciuto come formalismo di Newman-Penrose.

1.2 Perturbazione della Schwarzschild

Consideriamo una metrica della forma

gµν = g(0)µν + hµν (1.1)

dove |hµν | 1 e una piccola perturbazione e g(0)µν e la metrica di Sch-

warzschild nelle coordinate standard

g(0)µν dx

µdxν = −(

1− 2M

r

)dt2 +

(1− 2M

r

)−1

dr2 + r2dΩ2. (1.2)

Stiamo considerando quindi una metrica di background stazionaria (unvettore di Killing timelike) e statica (vettore di Killing timelike ortogonalealle ipersuperfici t = cost. ), ovvero che descrive un oggetto massivo sferi-camente simmetrico non rotante, ad esempio una pulsar oscillante e nonrotante. Il teorema di Birkhoff ci assicura che questa e l’unica soluzione delleequazioni di Einstein per un vuoto sfericamente simmetrico. Riusciamo ascrivere la metrica in questo modo perche pensiamo alla varieta considera-ta come il prodotto di due sottovarietaM = M2 × S2, con M2 spazio diMinkowski di dimensione 2 e S2 che indica una due sfera. Fogliando lospazio con una sottovarieta massimamente simmetrica posso scrivere lametrica come

ds2 = gABdxAdxB + r2Ωabdx

adxb (1.3)

Indichero da ora in poi le coordinate su M2 con indici latini maiuscoli,mentre con indici latini minuscoli indichero le coordinate di S2. Ωab ela metrica della due-sfera unitaria. Con DA e Da denotero le derivatecovarianti in M2 e S2.

1.2.1 Espansione in multipoli

E’ conveniente espandere la parte di metrica perturbata hµν in multipoliusando le armoniche sferiche Y (l,m)(θ, ϕ) e separare questi termini in basealle loro proprieta rispetto a trasformazioni di parita e trasformazioni dicoordinate sulla sfera unitaria. Rispetto a trasformazioni di parita (θ, ϕ)→(π − θ, π + ϕ), ovvero riflessioni rispetto all’origine, diremo che i multipolisono pari o polari se trasformano con (−1)l oppure le chiameremo dispario assiali se trasformano con (−1)l+1. Se moltiplichiamo la perturbazione

2


per la corrispondente armonica sferica, le perturbazioni pari rimangonoinvariate mentre quelle dispari cambiano segno.

Sotto trasformazioni di coordinate sulla sfera unitaria, le componenti dihµν cambiano come scalari, vettori e tensori. Dobbiamo quindi introdurrenuovi enti matematici a partire dalle armoniche sferiche. Le componentihAB trasformano come tre scalari e quindi le possiamo sviluppare in terminidelle usuali armoniche sferiche. Le componenti hAb trasformano comedue vettori e hab come un tensore. I vettori armonici sono definiti come ilgradiente delle armoniche sferiche

Y l,ma := DaY

l,m, (1.4)

mentre i vettori armonici dispari sono definiti come

X l,ma := −εbaDb (1.5)

dove ε e il tensore di Levi-Civita su una due sfera.Tensori armonici pari possono essere costruiti in due differenti modi, sia

moltiplicando la metrica della due sfera per le armoniche sferiche scalari,sia prendendo la derivata covariante seconda delle armoniche sferichescalari:

ΩabYl,m, DaDbY

l,m. (1.6)

Sorge un problema pero, queste ultime funzioni non formano un set linear-mente indipendente. Invece di sfruttare∇a∇bY l,m possiamo usufruire deltensore armonico di Zerilli-Mathews cosi definito

Z l,mab := DaDbYl,m +

1

2l(l + 1)ΩabY

l,m (1.7)

Quest’ultimo e un tensore a traccia nulla. Anche i tensori armonici disparisono a traccia nulla, ma fortunatamente possono essere scritti in un unicomodo

X l,mab =

1

2

(DaX

l,mb +DbX

l,ma

). (1.8)

Da notare che essendo Y 0,0 una costante i vettori armonici saranno nonnulli solo per l ≥ 1 mentre i tensori armonici per l ≥ 2.

Riassumendo possiamo scrivere la perturbazione della metrica, consi-derando il settore pari nel modo seguente

(hl,mµν )e =

(H l,mABY

l,m H l,mA Y l,m

b

HAYl,mb r2(K l,mΩabY

l,m +Gl,mZ l,mab )

)(1.9)

e il settore dispari

(hl,mµν )o =

(0 hl,mA X l,m

b

hl,mA X l,mb hl,mX l,m

ab )

)(1.10)

dove in generale i coefficienti (H l,mAB , H

l,mA ,K l,m, Gl,m, hl,mA , hl,m) sono fun-

zioni di r e t.

3


1.2.2 Perturbazioni Pari

I coefficenti dello sviluppo sono chiaramente dipendenti dalle coordi-nate ed in particolare cambiano per una trasformazione di gauge dellaforma

xµ → xµ + ξµ. (1.11)

con |ξ|<< 1. Un cambiamento di questo tipo conduce a

hµν → hµν −∇µξν −∇νξµ. (1.12)

Per visualizzare il perche di questa relazione possiamo pensare a due va-rieta che chiamero “spazio-tempo di background”Mb e “spazio-tempofisico”Mp e un diffeomorfismo ϕ : Mb −→ Mp. Considerando un dif-feomorfismo, dico che queste due varieta sono la “stessa cosa”. Possiamodefinire la perturbazione come la differenza tra il pull-back della metricafisica e la metrica di background:

hµν = (ϕ∗g)µν − g(0)µν . (1.13)

Da questa definizione non si evince nessuna ragione per cui le componentihµν debbano essere piccole: tuttavia se i campi gravitazionali suMp sonodeboli, allora per qualche diffeomorfismo avremo |hµν | 1. Delimitiamol’attenzione a quei diffeomorfismi per cui questo e vero. hµν ubbidira alleequazioni di Einstein linearizzate.

Consideriamo ora un campo vettoriale ξµ su Mb, questo campo cipermettera di descrivere il grande numero di diffeomorfismi, che ripettanola condizione sopra detta, tra le due varieta. Questo campo vettorialegenerera una famiglia ad un parametro di diffeomorfismi: ψε :Mb −→Mb.

Possiamo definire una famiglia ad un parametro di perturbazioni para-metrizzate da ε:

hεµν = [(ϕ ψε)∗g]µν − g(0)µν (1.14)

= [ψε∗(ϕ∗g)]µν − g(0)µν . (1.15)

4


La seconda eguaglianza deriva dal fatto che il pullback di una composizionee dato dalla composizione dei pullbacks in ordine inverso. Inserendo questarelazione in (1.13), troviamo:

hεµν = ψε∗(h+ g(0))µν − g(0)µν (1.16)

= ψε∗(hµν) + ψε∗(g(0)µν )− g(0)

µν . (1.17)

Usiamo adesso il fatto che ε e piccolo: in questo caso ψε∗(hµν) sara ugualea hµν al piu basso ordine mentre gli altri due termini danno la derivata diLie:

hεµν = ψε∗(g(0)µν ) + ε

(ψε∗(g

(0)µν )− g(0)

µν

ε

)(1.18)

= hµν + ε£ξg(0)µν (1.19)

= hµν + 2ε∂(µξν). (1.20)

Se prendiamo |ξµ| 1 possiamo scrivere

hµν = hµν −∇µξν −∇νξµ. (1.21)

Possiamo costruire combinazioni gauge-invarianti dei coefficienti del-lo sviluppo della perturbazione (H l,m

AB , Hl,mA ,K l,m, Gl,m, hl,mA , hl,m). Il pro-

cedimento e semplice, consideriamo una trasformazione di gauge pariξµ = (ξA, ξa), che puo essere espansa in termini di scalari e vettori armonici

ξA =∑l,m

El,mA Y l,m, ξa =∑l,m

El,mY l,ma . (1.22)

Una volta calcolati i simboli di Christoffel per la nostra metrica e le derivatecovarianti dei vettori sopra (1.22) possiamo calcolare l’espressione (1.21)da cui si ottengono le variazioni dei coefficienti dello sviluppo. Per i det-tagli consultare [2]. Questi risultati possono essere utilizzati per costruirequantita gauge-invarianti, due di queste combinazioni invarianti sono [13].

K l,m := K l,m +1

2l(l + 1)Gl,m − 2

rrAεl,mA , (1.23)

H l,mAB := H l,m

AB −DAεl,mB −DBε

l,mA (1.24)

conεl,mA := H l,m

A − 1

2r2DAG

l,m. (1.25)

Invece di cercare combinazioni gauge invarianti possiamo lavorarefissando uno specifico gauge. Possiamo scegliereEl,m eEl,mA tali cheGl,m =

H l,mA = 0. Questo e conosciuto come Regge-Wheeler gauge, e con questa scelta

abbiamo K l,m = K l,m, H lmAB .

In termini delle perturbazioni gauge invarianti K l,m e H l,mAB e possibile

definire quella che viene chiamata Zerilli-Moncrief master function

Ψl,mpari :=

2r

l(l + 1)

[K l,m +

2

Λ

(rArBH l,m

AB − rrADAK

l,m)], (1.26)

5


con Λ := (l − 1)(l + 2) + 6M/r. L’importanza di questa funzione derivadal fatto che usando le equazioni di campo di Einstein perturbate si puomostrare che la Ψl,m

pari obbedisce all’equazioni di Zerilli o piu generalmenteeven-parity master equation

∂2t Ψl,m

even − ∂2r∗Ψ

l,meven + V l,m

evenΨl,meven = Sl,meven, (1.27)

con r∗ coordinata radiale della tartaruga definita come r∗ := 1+2Mln(r/2M−1), e dove V l,m

even e Sl,meven sono termini di potenziale e sorgente che dipendonoda l e m. L’espressione esplicita la si puo trovare in [7].

1.2.3 Perturbazioni Dispari

Come per le perturbazioni pari possiamo costruire quantita gauge inva-rianti anche per le perturbazioni dispari. Si procede come prima prendendopero una trasformazione dispari ξµ = (ξA, ξa), che puo essere espansa intermini di scalari e vettori armonici dispari

ξA = 0, ξa =∑l,m

El,mX l,ma . (1.28)

Calcolando le nuove derivate covarianti si perviene alle seguenti trasfor-mazioni

hl,mA → hl,mA −DAEl,m +

2rArEl,m, (1.29)

hl,m → hl,m − 2El,m. (1.30)

Come prima possiamo trovare quantita gauge invarianti

hl,mA := hl,mA − 1

2DAh

l,m +rArhl,m. (1.31)

E’ chiaro che possiamo sempre scegliere El,m tale che hl,m = 0 e quindihl,mA = hl,mA , corrispondentemente al Regge-Wheeler gauge.

In termini di hl,mA definiamo la Cunningham-Price-Moncrief master func-tion

Ψl,modd :=

2r

(l − 1)(l + 2)εAB

[DAh

l,mB − 2rA

rhl,mB

](1.32)

:=2r

(l − 1)(l + 2)εAB

[DAh

l,mB − 2rA

rhl,mB

]. (1.33)

Come prima e possibile mostrare che Ψl,modd ubbidisce a

∂2t Ψl,m

even − ∂2r∗Ψ

l,meven + V l,m

evenΨl,meven = Sl,mevenΨl,m

even (1.34)

dove V l,meven e Sl,meven sono termini di potenziale e sorgente che dipendo-

no da l e m. L’equazione sopra e chiamata Regge-Wheeler equation o piugeneralmente odd-parity master equation.

6


1.2.4 TT Gauge

Le master function descritte sopra sono importanti perche non solo soddi-sfano ad una equazione tipo quella delle onde, ma anche perche si possonolegare con le ampiezze di polarizzazione “piu”e “croce”. ovvero h+ e h× del-le onde gravitazionali nel transverse and traceless (TT) gauge. Seguo il ragiona-mento di [8] per trovare la relazione sopra enunciata. Per prima cosa ci ser-viamo di una tetrade le cui componenti sono eµµ ≡ e

b, e−b, r−1, (rsinθ)−1,ovvero decomponiamo la perturbazione su di una tedrade e diagonale diun osservatore stazionario, ricordando che vale hµµ = eµµeννhµν . Secondacosa, tutte le quantita saranno calcolate lontane dalla sorgente, ovvero nellazona chiamata radiation gauge. Non considero nel seguito i coefficienti H l,m

AB

eH l,mA essendo possibile mostrare che decadono piu rapidamente degli altri.

Nella nostra base vale

h+ =1

2(hθθ − hϕϕ) =

1

2r2

(hθθ −

hϕϕsin2θ

)(1.35)

h× = hθϕ =hθϕ

r2sinθ(1.36)

A questo punto devono essere usate le espansioni in multipoli, in particola-re per il settore pari e possibile scrivere

(h+)l,meven =Gl,m

2

(Z l,mθθ −

Z l,mϕϕsin2θ

)(1.37)

(h×)l,meven = Gl,m(Z l,mθϕsinθ

)(1.38)

Se consideriamo pero la condizione di annullamento della traccia abbiamo

Ωab(K l,mΩabY

l,m +Gl,mZ l,mab)

= 0. (1.39)

Come detto precedentemente Z l,mab e a traccia nulla, il che implica K l,m = 0.Lafunzione di Zerilli-Moncrief definita in (1.26) allora si semplifica

Ψl,meven = rGl,m. (1.40)

Possiamo quindi riscrivere (1.37) come

(h+)l,meven =Ψl,meven

2r

(Z l,mθθ −

Z l,mϕϕsin2θ

)=

Ψl,meven

r

[∂2

∂θ+

1

2l(l + 1)

]Y l,m

(1.41)

(h×)l,meven =Ψl,meven

r

(Z l,mθϕsinθ

)=

Ψl,meven

r

( im

sinθ

)[ ∂∂θ− cotθ

]Y l,m,

(1.42)

usando la relazione ∂ϕY l,m = imY l,m.

7

1.3. Gli scalari di Weyl

Per quanto riguarda il settore dispari della perturbazione

(h+)l,modd =hl,m

2

(X l,mθθ −

X l,mϕϕ

sin2θ

)(1.43)

(h×)l,modd = hl,m(X l,m

θϕ

sinθ

)(1.44)

Dobbiamo legare hl,m a Ψl,modd. Si procede analogamente a prima, ma e

necessario usare maggior cautela dato che Ψl,modd dipende solo da hl,mA e non

da hl,m. Tuttavia queste quantita sono legate nel TT gauge e sfruttando lacondizione di trasversalita

∇µhµa = 0, (1.45)

in cui i simboli di Christoffel devono essere riscritti utilizzando l’espansionein multipoli, sostituendo la divergenza di X l,m

ab e la condizione di spazialitahµ0 = 0 troviamo

hl,m ∼ −rΨl,modd. (1.46)

Possiamo riscrivere la perturbazione dispari nel modo seguente

(h+)l,modd = −Ψl,modd

2r

(X l,mθθ −

X l,mϕϕ

sin2θ

)= −

Ψl,modd

r

( im

sinθ

)[ ∂∂θ− cotθ

]Y l,m

(1.47)

(h×)l,modd =Ψl,modd

r

(X l,mθϕ

sin2θ

)=

Ψl,meven

r

[∂2

∂θ+

1

2l(l + 1)

]Y l,m,

(1.48)

In definitiva unendo (1.41), (1.42), (1.47) e (1.48)si puo scrivere

h+ =1

2r

∑l,m

[Ψl,meven

(Z l,mθθ −

Z l,mϕϕsin2θ

)−Ψl,m

odd

(X l,mθθ −

X l,mϕϕ

sin2θ

)](1.49)

h× =1

r

∑l,m

[Ψl,meven

(Z l,mθϕsinθ

)−Ψl,m

odd

(X l,mθϕ

sin2θ

)]. (1.50)

1.3 Gli scalari di Weyl

Un secondo approccio per l’estrazione di informazione delle ondegravitazionali e basato sulla definizione del tensore di Weyl

Cαβµν := Rαβµν −2

n− 2

[gα[µRν]β − gβ[µRν]α

]+

2

(n− 1)(n− 2)gα[µgν]βR.

(1.51)

La definizione di questo tensore deriva dal fatto che si e voluto cercareun modo per riscrivere il tensore di Riemann Rαβµν in modo da separare

8

1.3. Gli scalari di Weyl

l’informazione portata dal tensore e dallo scalare di Ricci. Il tensore di Weylha le stesse simmetrie del tensore di Riemann, ma ha l’ulteriore proprietadi avere traccia nulla

Cαµαν = 0. (1.52)

In n dimensioni il tensore di Weyl possiede n(n + 1)(n + 2)(n − 3)/12

componenti indipendenti, ovvero 10 componenti in quattro dimensioni e siannulla per n‘3.

Un’altra importante proprieta del tensore di Weyl e legata al suo com-portamento sotto trasformazioni conformi della metrica del tipo

gij = Ωgij . (1.53)

Le due metriche gij e gij hanno in generale due tensori di curvatura diRiemann, ma hanno lo stesso tensore di Weyl

Cαβµν = Cαβµν . (1.54)

Serve tuttavia un ulteriore ingrediente, la tedrade nulla di Newman-Penrose.

1.3.1 La tetrade nulla di Newman-Penrose

L’idea Newman e Penrose e di introdurre una tedrade di vettori nulli.Il primo passo da compiere e trovare una tedrade ortonormale ~e(a) inmodo da poter riscrivere la metrica nel seguente modo

gµν = −e(0)µe(0)ν + e(1)µe(1)ν + e(2)µe(2)ν + e(3)µe(3)ν (1.55)

e~e(a) · ~e(b) = η(a)(b) (1.56)

con η(a)(b) corrispondente alla matrice di Minkowski.Solitamente si sceglie eµ(0) come il vettore normale all’ipersuperfice

spaziale eµ(0) = nµ, eµ(1) come il vettore radiale in coordinate sferiche eµ(1) =

eµ(r). Dato che in generale eµ(θ) e eµ(ϕ) non sono vettori unitari e ortogonali,per cui e necessario applicare la procedura di ortonormalizzazione di Gram-Schmidt per costruire eµ(2) ed eµ(3).

E’ ora possibile iniziare a costruire due vettori nulli

lµ :=1√2

(eµ(0) + eµ(1)

)kµ :=

1√2

(eµ(0) − e

µ(1)

)(1.57)

Per costruire altri due vettori nulli con cui formare una tedrade e necessariopassare al campo complesso:

mµ :=1√2

(eµ(2) + ieµ(3)

) −mµ

:=1√2

(eµ(2) − ie

µ(3)

). (1.58)

. I vettori sopra definiti formano una tedrade nulla

lµlmu = kµkmu = mµmmu =−mµ −mmu = 0. (1.59)

9

1.4. Energia e momento delle onde gravitazionali

Come precedentemente detto lo scalare di Weyl ha 10 componenti indi-pendenti e nel formalismo di Newman-Penrose appena introdotto queste10 componenti possono essere rappresentate da cinque scalari complessi,gli scalari di Weyl cosi definiti

Ψ0 := C(1)(3)(1)(3) = Cαβµν lαmβlµmν (1.60)

Ψ1 := C(1)(2)(1)(3) = Cαβµν lαkβlµmν (1.61)

Ψ2 := C(1)(3)(4)(2) = Cαβµν lαmβmµkν (1.62)

Ψ3 := C(1)(2)(4)(2) = Cαβµν lαkβmµkν (1.63)

Ψ4 := C(2)(4)(2)(4) = Cαβµνkα −m

βkµmν . (1.64)

E’ da notare che questi sono scalari rispetto a trasformazioni di coordinate,ma dipendono dalla tetrade nulla scelta. Nella prossima sezione verraillustrato come queste entita si legano alle onde gravitazionali.

1.4 Energia e momento delle onde gravitazionali

In relativita generale non c’e la nozione di densita di energia del campogravitazionale stesso. Questo puo essere facilmente capito ricordando chenella teoria Newtoniana la densita di energia del campo gravitazionalee proporzionale a∇iϕ∇iϕ, con ϕ potenziale gravitazionale. Dato che nellimite Newtoniano ϕ gioca il ruolo della metrica, ci aspetteremmo cheanche in relativita generale l’energia della gravita si potesse scrivere co-me combinazione di termini quadratici della derivata prima della metrica.Considerando che le derivate covarianti sono combinazioni di simboli diChristoffel che localmente si annullano in coordinate piatte, un termine co-me quello considerato non puo essere corretto. E’ possibile tuttavia, definireun flusso di energia e momento portato via dalle onde gravitazionali.

Ricordiamo anzitutto che le onde gravitazionali sorgono al primo ordinedi perturbazione della metrica di background. Ricordiamo anche che iltensore energia impulso associato a campi scalari o elettromagnetici e delsecondo ordine nei campi. Questo suggerisce di sondare al secondo ordineperturbativo la metrica di background se vogliamo trovare il nostro tensoreenergia-impulso:

gµν = g(0)µν + εh(1)

µν + ε2h(2)µν (1.65)

dove g(0)µν e la metrica di background (che assumiamo essere una solu-

zione delle equazioni di Einstein nel vuoto), e dove abbiamo introdottoesplicitamente il parametro ε 1 per tracciare l’oridine dei vari termini.

Il tensore di Ricci si puo espandere nella forma

Rµν = εR(1)µν + ε2R(2)

µν (1.66)

10


dove ho posto gia R(0)µν = 0, vero per costruzione. Un calcolo un po’ lungo

porta a riscrivere i differenti pezzi del tensore di Ricci nel seguente modo

R(1)µν = F (1)

µν (h(1)), (1.67)

R(2)µν = F (1)

µν (h(2)) + F (2)µν (h(1)), (1.68)

dove F (1)µν e F (2)

µν sono definiti come

F (1)µν :=

1

2(−h|µν − hµν|αα + hαµ|ν

α + hαν|µα), (1.69)

F (2)µν :=

1

4hαβ|µh

αβ|ν +

1

2hαβ(hαβ|µν + hµν|αβ − 2hα(µ|ν)β)

+ hα|βν hµ[α|β] +1

2

(hαβ|β −

1

2h|α)

(hµν|α − 2hα(µ|ν))

(1.70)

Consideriamo ora le equazioni di Einstein nel vuoto per ogni potenza di ε:

R(0)µν = 0 (1.71)

R(1)µν = F (1)

µν (h(1)) = 0 (1.72)

R(2)µν = F (1)

µν (h(2)) + F (2)µν (h(1)) = 0. (1.73)

La prima di queste equazioni e automaticamente soddisfatta, mentre laseconda descrive la dinamica dell’onda gravitazionale. La terza equazionepuo essere riscritta in termini del tensore di Einstein nel modo seguente

F (1)µν (h(2))− 1

2g0µνF

(1)(h(2)) = 8πtµν , (1.74)

con

tµν := − 1

8π

(F (2)µν (h(1))− 1

2g(0)µν F

(2)(h(1))

). (1.75)

Le ultime relazioni suggeriscono quindi di interpretare tµν come il tenso-re energia-impulso dell’onda gravitazionale descritta da h(1)

µν . Sfruttandol’identita di Bianchi all’ordine ε2 si puo vedere che

tµν|µ = 0, (1.76)

cosi che tµν risulta essere una quantita conservata.

1.4.1 La media di Isaacson

La quantita tµν non e gauge invariante, per questo motivo e chiamatapseudo-tensore energia-impulso. Questo non vuol dire che debba perderedi significato. Per quanto sia possibile annullare in un punto tµν , nonsara in generale possibile annullarlo in una regione finita dello spazio-tempo. Dobbiamo cercare di mediare questo tensore su di una regione finitadella nostra varieta. E’ necessario introdurre la cosiddetta short wavelengthapproximation, ovvero consideriamo una media di tµν su di una regioneche comprende diverse lunghezze d’onda, ma comunque piccola rispettoalla lunghezza d’onda della metrica di background. Isaacson mostra comecostruire questa media in [10].

11


Il risultato di un integrazione di un campo tensoriale su di uno spaziocurvo potrebbe non essere piu un tensore, dato che ad ogni punto questocampo ha differenti proprieta di trasformazione. Dato che e permessosommare tensori ad uno stesso punto e necessario costruire una media conqualche cosa che porti tutti i tensori in un unico punto comune. Per farecio in modo unico deve essere introdotto il bivettore di trasporto parallelosu geodetiche denotato con gβ

′α (x, x′). Questo trasforma come un vettore sia

per trasformazioni di x che x′ e assumendo che questi due punti sianosufficientemente vicini possiamo essere sicuri dell’esistenza di un’unicageodetica, della metrica che indicheremo con γαβ , che li collega. Dato ilvettore Aβ′ in x′, allora Aα = gβ

′α Aβ′ e l’unico vettore in x che puo essere

ottenuto trasportando parallelamente Aβ′ in x′ a x lungo la geodetica. Datoun tensore Tµν oscillante con frequenza λ e una metrica di background confrequenza di oscillazione L λ, definiamo la media del tensore Tµν come

< Tµν(x) >:=

∫tutto lo spazio

gα′

µ (x, x′)gβ′

ν Tα′β′(x′)f(x, x′)

√−γ(0)(x′)d4x′

(1.77)dove f(x, x′) e una funzione peso che decresce a zero in modo smoothquando x e x′ differiscono di una distanza d tale che λ d L e dove∫

tutto lo spazio

f(x, x′)√−γ(0)(x′)d4x′ = 1 (1.78)

Grazie a f(x, x′) possiamo considerare la regione dove avviene effettiva-mente la media, a metrica piatta e non dobbiamo aggiungere ulterioricondizioni per ottenere l’unicita della geodetica.

Questa media possiede interessanti proprieta:

1. le derivate covarianti commutano; es. < hhµν|αβ >=< hhµν|βα >.L’errore commesso va come O( λR)2 con R raggio di curvatura delbackground.

2. il gradiente media a zero; es. < (h|αhµν)|β >= 0. L’errore commessova come O( λR) con R raggio di curvatura del background.

3. integrando per parti si trova il salto delle derivate; es. < hh|µν >=<

h|νh|µ >

Con questo nuovo mezzo possiamo calcolare il valor medio del tensoreenergia impulso

Tµν :=< tµν >=1

32π

⟨hαβ|µh

αβ|ν −

1

2hµhν − 2hαβ|β hα(µ|ν)

⟩, (1.79)

dove ho tralasiato l’apice (1) in h(1) e hµν e definito come

hµν := hµν −1

2g(0µν)h. (1.80)

12


Il tensore Tµν sopra definito e chiamato Isaacson stress energy tensor. Tµνe gauge invariante e nel particolare caso del (TT) gauge in cui h = hαβ|β = 0,l’Isaacson stress energy tensor si riduce a

Tµν =1

32π

⟨hαβ|µh

αβ|ν⟩. (1.81)

1.4.2 Energia e momento irradiati

Assumiamo di essere lontani dalla sorgente, cosi da poter descrivere leonde gravitazionali come onde sferiche che si propagano in un backgroundpiatto. In coordinate locali Cartesiane, usando il TT gauge, la perturbazionedella metrica puo essere riscritta nel seguente modo

hµν = h+A+µν + h×A×µν , (1.82)

dove h+,× sono le ampiezze delle tue indipendenti polarizzazzioni dell’on-da e A+,×

µν sono tensori costanti simmetrici tali che

Aµν lν = 0, Aµmu = 0, Aµνu

ν = 0, (1.83)

con lµ vettore nullo (vettore d’onda) e uµ arbitrario vettore unitario di tipotempo.

Voglio riscrivere il tensore energia impulso in un modo piu convenien-te e per fare cio e necessario considerare una base spaziale ortonormale(er, eθ, eϕ), indotta dalle coordinate sferiche (r, θ, ϕ), un vettore di tipo tem-po uµ e lµ vettore radiale nullo uscente (se la sorgente e isolata avro soloonde uscenti). Cio che si trova e che le uniche componenti non nulle deltensore di polarizzazione dell’onda sono

A+

θθ= A+

ϕϕ = 1, (1.84)

A×θϕ

= A×ϕθ

= 1. (1.85)

Il tensore energia impulso puo quindi essere riscritto

Tµν =1

16π〈∂µh+∂νh

+ + ∂µh×∂νh

×〉, (1.86)

o in modo equivalente

Tµν =1

16πRe〈∂µH∂ν

−H〉, (1.87)

con H = h+ − ih× e dove Re(z) denota la parte reale di z.La cosa interessante e che la quantita complessa H si puo esprimere in

termini degli scalari di Weyl. Per vedere cio e necessario osservare che iltensore di Weyl e quello di Riemann coincidonolontani dalla sorgente delleonde gravitazionali. In particolare per un onda descritta nel TT gauge chesi propaga nella direzione radiale gli scalari di weyl diventano

Ψ1 = Ψ2 = Ψ3 = 0 , (1.88)

13


Ψ0 =− 1

4(∂2t h

+ + 2∂t∂rh+ + ∂2

rh×)

− i

4(∂2t h× + 2∂t∂rh

× + ∂2rh×) ,

(1.89)

Ψ4 =− 1

4(∂2t h

+ − 2∂t∂rh+ + ∂2

rh×)

+i

4(∂2t h× − 2∂t∂rh

× + ∂2rh×) .

(1.90)

Per un’onda uscente abbiamo h = h(r − t), cosi che ∂rh = −∂th. Gli scalaridi Weyl diventano

Ψ0 = Ψ1 = Ψ2 = Ψ3 = 0 , (1.91)

Ψ4 = −h+ + ih× = −H. (1.92)

Tutto cio implica che

H = −∫ t

−∞

∫ t′

−∞Ψ4dt

′′dt′. (1.93)

Adesso ci sono tutti gli elementi che servono per calcolare flussi dienergia e momento. Consideriamo il flusso di energia lungo la direzionedelle i, la quale e data in generale da T 0i. Il flusso lungo la direzione radialesara quindi dato da

dE

dtdA= T 0r =

1

16πRe〈∂0H∂rH〉 =

1

16πRe〈∂tH∂rH〉 , (1.94)

con dA elemento di area ortogonale alla direzione radiale. Per onde uscentipossiamo riscrivere l’ultima relazione nel modo seguente

dE

dtdA=

1

16πRe〈H ˙H〉 =

1

16πRe〈|H|2〉 . (1.95)

Per il flusso di momento si procede in modo analogo considerando perola componente spaziale del tensore energia impulso Tij

dPidtdA

= Tir =1

16πRe〈∂iH∂rH〉 . (1.96)

Se usiamo il fatto che distanti dalla sorgente la componente angolare puoessere trascurata e che ∂iH ≈ (xi/r)∂rH troviamo che

dPidtdA

≈1

16πli〈|

.H|2〉 , (1.97)

con ~l vettore radiale unitario nello spazio piatto.

14

CAPITOLO 2

LA SORGENTE DELLE ONDE

Indice2.1 Stella di neutroni . . . . . . . . . . . . . . . . . . . . . . 152.2 Soluzione di TOV . . . . . . . . . . . . . . . . . . . . . 162.3 EOS . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 172.4 Stabilita . . . . . . . . . . . . . . . . . . . . . . . . . . . 18

All’alta fantasia qui manco possa; — ma gia volgeva ilmıo disio e ’l velle — sı come rota ch’igualmente emossa — l’amor che muove il sole e l’altre stelle.

DANTE ALIGHIERI

2.1 Stella di neutroni

Nella simulazione eseguita e coinvolta una stella di neutroni, un oggettocon massa dell’ordine pari a quella solare, che si trova ad avere densitanucleare ed e globalmente neutro. Il raggio di un tale oggetto si aggiraintorno alla decina di chilometri. Si pensa essere lo stadio finale di unastella con massa compresa tra 1, 4 e 3 masse solari. A causa dell’altissimadensita e delle piccole dimensioni, una stella di neutroni possiede un campogravitazionale superficiale cento miliardi di volte piu intenso di quelloterrestre.

E’ proprio questo oggetto che sfruttiamo come sorgente di onde gravi-tazionali.

Durante il collasso gravitazionale, dovuto al fatto che la stella ha finitoil suo carburante per le reazioni nucleari, la stella puo dare origine a treoggetti compatti: nane bianche, stella di neutroni, buchi neri. Gli elettroninon potranno cadere, per il principio di esclusione di Pauli, tutti nello statofondamentale. Questo fenomeno genera una pressione diversa da zeroche sostiene la pressione dovuta al collasso gravitazionale e vengono cosi

15

2.2. Soluzione di TOV

generate le nane bianche. Nel casso la pressione gravitazionale fosse piuforte, gli elettroni andranno ad annichilarsi con i protoni dei nuclei atomicie rimarranno neutroni, in questo caso la sorgente di pressione maggiore edovuta alla forza repulsiva tra i nucleoni. Non si sa con esattezza quelloche accade pero almeno il guscio superficiale dovrebbe essere formato daneutroni (da qui il nome stella di neutroni), l’interno della stella e ancoraargomento di discussione. Quello che si sa e che la densita e pari a quelladel nucleo e l’oggetto e globalmente neutro. Se la massa della stella loconsente, vengono generati buchi neri.

La stella di neutroni utilizzata come sorgente e una stella non rotantee non sono stati presi in considerazione neanche i forti campi magneticipresenti in oggetti di questo tipo. Lo scopo di questo lavoro era semplice-mente testare i metodi di estrazione delle onde. Per simulare l’interno dellastella abbiamo adottato un’equazione di stato di tipo politropico.

2.2 Soluzione di TOV

E’ evidente che le equazioni di Einstein nel vuoto, che caratterizza-no lo spazio attorno alla stella, dovranno essere affiancate da equazioniche coinvolgono materia. La soluzione delle equazioni di Einstein perun corpo a simmetria sferica e statico e nota ed e chiamata soluzione diTolman-Oppenheimer-Volkoff, in breve soluzione di TOV.

Schematicamente, si trova che la geometria dello spazio tempo all’inter-no di una stella sferica, statica, fluida e data da

ds2 = −e2φdt2 +(

1− 2m(r)

r

)−1dr2 + r2dΩ2, (2.1)

dovem(r) = 4π

∫ r

0ρ(r

′)r′2dr

′(2.2)

e φ e determinato dadφ

dr=m(r) + 4πr3p

r[r − 2m(r)]. (2.3)

La condizione necessaria e sufficiente per l’equilibrio e che la seguenteequazione, chiamata equazione di Tolman-Oppenheimer-Volkoff per l’equilibrioidrostatico

dp

dr= (p+ ρ)

m(r) + 4πr3p

r[r − 2m(r)], (2.4)

sia soddisfatta.Ricavare le precedenti equazioni non e troppo complicato. Per quanto

detto nel primo capitolo, la metrica dovra avere la seguente forma

ds2 = f(r)dt2 + h(r)dr2 + r2(dθ2 + sin2θdφ2). (2.5)

Stiamo considerando che la stella sia composta da un fluido perfetto il cuitensore energia impulso e dato da

Tµν = ρuµν + p(gµν + uµunu). (2.6)

16

2.3. EOS

Per compatibilita con le nostre condizioni di simmetria dobbiamo porre

uµ = −(e0)µ = −f12 (dt)µ. (2.7)

Grazie alle precedenti condizioni e possibile calcolare il tensore e loscalare di Ricci e scrivere le opportune equazioni di Einstein in cui compa-rira anche il tensore energia impulso della materia fluida. Dall’equazioneG00 = 8πT00 non comparendo f si ricava h in funzione di m(r). Il motivoper cui f non compare e semplice: l’equazione considerata rappresenta uninitial value constraint. Risolvendo l’intero sistema si trovano la metrica ele condizioni scritte all’inizio del paragrafo.

Da notare pero che m(r) non rappresenta la massa vera e propria, datoche nell’integrale non si e usato l’elemento di volume

√gdx3. Possiamo

definire la massa propria Mp come

Mp = 4π

∫ R

0ρ(r)r2

[1− 2m(r)

r

]− 12dr (2.8)

La differenza tra M e Mp puo essere interpretata come energia di legamegravitazionale, EB , della configurazione,

EB = Mp −M (2.9)

che e sempre positiva essendo Mp > M .Vale la pena ricordare che l’equazione di Tolman-Oppenheimer-Volkoff

integrate porta in modo naturale a definire una massa limite per la stabilitadella stella.

Quello che pero e realmente importante ai fini di questo lavoro e capi-re come queste equazioni devono essere utilizzate nella simulazione. Permateria fluida con una determinata equazione di stato, p = p(ρ), la confi-gurazione di equilibrio puo essere determinata come segue: prescriviamoarbitrariamente una densita centrale ρc, e quindi una pressione centralep(ρc). Successivamente devono essere integrate le equazioni (2.2) e (2.4)fino a quando non si raggiunge la superficie della stella, dove p = ρ = 0 esi deve raccordare la soluzione con quella del vuoto di Schwarzschild ( 1.2).Risolviamo l’equazione per φ (2.3) richiedendo sempre che la metrica siinterlacci con quella nel vuoto sulla superficie della stella.

2.3 EOS

Le equazioni di evoluzione dell’idrodinamica relativistica, cosi comequelle di TOV non sono sufficienti da sole a descrivere il nostro sistema.Abbiamo bisogno di un ulteriore condizione che puo essere data dall’equa-zione di stato (EOS). Teoricamente quest’ultima dovrebbe essere ricavatada considerazioni microscopiche, ma ancora un po’ di foschia aleggia sullaconoscenza esatta di cio’ che avviene internamente ad una stella di neutroni.La conoscenza dell’interazione di materia nucleare ad elevate densita none ancora ben compresa. Esistono pero semplici equazioni di stato usatenelle simulazioni numeriche che risultano essere molto utili in particolare

17

2.4. Stabilita

l’equazione politropica. Questa e strettamente legata al gas ideale ed ha laseguente forma:

p = KρΓ0 ≡ Kρ

1+ 1N

0 , (2.10)

dove K e una costante e N e chiamato indice politropico. L’origine di questaeos e da ricercarsi in un processo adiabatico, ovvero senza scambio dicalore, per il quale si ricava, dalla prima legge della termodinamica

p = Kργ0 (2.11)

con γ =cpcV

rapporto tra calore specifico a pressione costante e a volumecostante. Confrontando questa equazione con (2.10) notiamo che Γ giocail ruolo di γ. E’ da sottolineare pero il fatto che l’equazione politropicae promossa ad eos anche per sistemi per cui avviene, o comunque puoavvenire, scambio di calore. La costante K e invece legata all’entropia delfluido. Piu precisamente, se σ e l’entropia per unita di volume possiamoscrivere

K = eσ/cV . (2.12)

Per maggiori dettagli consultare [2].Nel caso di oggetti compatti freddi, come una stella di neutroni, che si

possono formare grazie al principio di esclusione di Pauli, non potremmousare un equazione di stato di un gas ideale come descritto sopra, datoche applicandosi ad un caso classico non puo descrivere accuratamente ilivelli degeneri di un gas di Fermi. Tuttavia si riesce a trovare un equazionedi stato per un gas di Fermi a temperatura zero che si riduce proprioall’equazione politropica. Questo e vero sia nel caso non relativistico, chenel limite fortemente relativistico, con γ = 5/3 e γ = 4/3 rispettivamente.Piu generalmente si utilizzano equazioni politropiche con 1 < γ < 3 persemplici modelli di stelle di neutroni. Per i nostri scopi abbiamo utilizzatoK = e Γ =.

2.4 Stabilita

Non e sufficiente per una stella essere in equilibrio per poter esisterein natura. Deve esserci un equilibrio stabile. Le oscillazioni di una stellapossono essere scomposti in modi normali con frequenze definite. Se ilquadrato di tali frequenze e negativo, la stella sara instabile e non potrasopravvivere in natura. Per trovare queste oscillazioni e necessario risol-vere le equazioni di Einstein per piccole oscillazioni vicino all’equilibrio.Il calcolo non e molto complesso e da molte piu informazioni di quellenecessarie. Qui interessano solo un paio di fatti.

Si puo dimostrare che i modi radiali sono quelli che regolano la stabi-lita di una stella sferica. Un oscillazione sferica puo essere descritta dallospostamento radiale ξ(r, t) di un elemento fluido localizzato al raggio r altempo t. Gli spostamenti possibili possono essere analizzati tramite unospettro discreto di modi con frequenze definite ωnn = 1, 2, 3 .... Per ognimodo lo spostamento vale

ξn(r)e±iωnt. (2.13)

18

2.4. Stabilita

I modi stabili hanno ω2n > 0, hanno frequenze reali e oscillano. Modi

instabili hanno ω2n < 0, frequenze immaginarie e possono crescere esponen-

zialmente veloce. I modi devono annullarsi al centro della stella, mentresulla superficie devono preservare la condizione che la pressione svanisce.Come descritto in [19] se compissimo il calcolo per le gfrequenze troverem-mo un risultato in funzione della densita centrale dell’oggetto considerato,ω2n(ρc), n = 0, 1, 2, .... Schematicamente otterremmo qualche cosa simile

alla figura (2.1)

Figura 2.1: Rappresentazione schematica del quadrato delle frequenze infunzione della densita centrale di una stella, per i modi 0, 1, 2.

Al variare di ρc le frequenze dei vari modi cambiano facendo diventarela stella stabile o instabile. Un modo ha frequenza nulla quando cambiastabilita. Una frequenzapari a zero ha la proprieta di ottenero uno sposta-mento indipendente dal tempo, ha energia nulla e quindi deve conservarel’energia della massa totale. A basse densita tutti i modi sono stabili. AlladensitaaA il modo piu basso diventa instabile. Al puntoB anche il secondomodo diventa instabile. Queste frequenze riportano la stabilita nei puntiC e D. Le stelle sono stabili solo quando tutti i suoi modi hanno frequen-za positiva, quindi possiamo individuare delle configurazioni di densitaparticolari: 0 < ρc < A corrisponde a configurazioni che possono avere lenane bianche mentre tra C < ρc < D ci sono le configurazioni delle stelledi neutroni.

Dovremo vedere quindi anche la nostra stella simulata oscillare, chia-ramente le oscillazioni saranno differenti in funzione del raffinamento delreticolo che usato.

19

CAPITOLO 3

FORMALISMO 3+1 E FORMULAZIONE BSSNOK

Indice3.1 Il formalismo 3+1 . . . . . . . . . . . . . . . . . . . . . . . 213.2 Equazioni ADM . . . . . . . . . . . . . . . . . . . . . . 243.3 La formulazione BSSNOK . . . . . . . . . . . . . . . . 253.4 Condizioni di Gauge . . . . . . . . . . . . . . . . . . . . 27

3.4.1 Maximal slicing . . . . . . . . . . . . . . . . . . . 283.4.2 Condizione di slicing iperbolico . . . . . . . . . 283.4.3 Elliptic shift condition . . . . . . . . . . . . . . . 303.4.4 Driver Conditions . . . . . . . . . . . . . . . . . 32

Grazie al modello matematico posso dirvi Come e natol’universo: non chiedetemi il Perche.

STEPHEN HAWKING

3.1 Il formalismo 3+1

Le equazioni di Einstein sono scritte in modo covariante senza una nettadistinzione tra spazio e tempo. Il formalismo 3+1 permette una rappresenta-zione piu intuitiva che consente di pensare ad una evoluzione dinamicadel campo gravitazionale nel tempo.Suddividiamo lo spaziotempo in unospazio tridimensionale e uno monodimensionale che rappresenta il tempo.Questo e il formalismo piu utilizzato in relativita numerica.

La prima cosa da fare nello studio di un sistema fisico e porre in modoadeguato un problema di Cauchy: date le condizioni iniziali e al contorno, leequazioni fondamentali del sistema permettono di prevedere l’evoluzionefutura o passata del sistema. E’ necessario quindi riscrivere il le equazionidi Einstein come un problema di Cauchy ben posto cercando di suddividerebene i ruoli del tempo e dello spazio.

21

3.1. Il formalismo 3+1

Dobbiamo quindi considerare uno spazio globalmente iperbolico, ov-vero uno spazio che sia possibile fogliare con ipersuperfici di tipo spazio.Possiamo considerare queste ipersuperfici come livelli corrispondenti adun determinato parametro t, che puo essere considerato funzione univer-sale del tempo ( t non necessariamente corrispondera al tempo proprio).

Se consideriamo due ipersuperfici adiacenti Σt e Σt+dt, la geometriadello spazio tempo contenuta tra di esse puo essere descritta considerandoi seguenti tre ingredienti:

Figura 3.1: Due adiacenti superfici spacelike.

• La metrica tri-dimensionale γij (i, j = 1, 2, 3) che permetta la misura-zione della distanza propria dentro le ipersuperfici stesse:

dl2 = γijdxidxj . (3.1)

• La dilatazione del tempo proprio “lapse”, tra le ipersuperfici, mi-surato da un’osservatore che si mupve della direzione normale alleipersuperfici, ovvero da un osservatore Euleriano:

dτ = α(t, xi)dt, (3.2)

con α lapse function.

• La velocita relativa tra l’osservatore Euleriano e le linee di coordinataspaziale costante:

xit+dt = xit − βi(t, xj)dt, (3.3)

Il vettore βi e conoscito come shift vector.

La lapse function α e lo shift vector βi non sono uniche, sono funzionidi gauge che portano informazione sul sistema di coordinate scelto.

22

3.1. Il formalismo 3+1

La metrica dello spazio tempo si puo mostrare assumere la seguenteforma:

ds2 = (−α2 + βiβi)dt2 + 2βidtdx

i + γijdxidxj , (3.4)

con βi = γijβj .

Per riscrivere le equazioni di Einstein dobbiamo anche introdurre ilconcetto di curvatura estrinseca, ovvero la curvatura di una varieta chesi percepisce a causa della sua immersione in uno spazio con dimensionemaggiore. Dobbiamo distinguere l’intrinseca curvatura delle ipersuperficicon cui abbiamo foliato lo spazio con la curvatura estrinseca che questeassumono immergendole in uno spazio quadri-dimensionale. La curvaturaintrinseca e data dal tensore di Riemann tri-dimensionale definito in terminidella tre-metrica γi,j . La curvatura estrinseca si puo esprimere in termini dicio che succede al vettore normale alla ipersuperfice quando viene traspor-tato parallelamente da un punto ad un altro. La curvatura estrinseca Kαβ euna misura di quanto varia il vettore normale sotto trasporto parallelo.

Se l’operatore di proiezione sulle ipersuperfici spaziali e

Pαβ := δαβ + nαnβ, (3.5)

con nα vettore normale alle iper-superfici, la curvatura estrinseca e cosidefinita

Kµν := −Pαµ∇αnν . (3.6)

Si puo verificare che la curvatura estrinseca non e altro che la “velocita”della metrica spaziale vista da un osservatore Euleriano, infatti si puoriscrivere come

Kµν = −1

2£~nγµν . (3.7)

Essendo ~n normale alle ipersuperfici, per ogni scalare φ si puo scrivere

£~nαγµν =1

φ£φ~nγµν . (3.8)

Prendendo come funzione scalare il lapse α possiamo scrivere

Kµν = − 1

2α£α~nγµν = − 1

2α

(£~t −£~β

)γµν , (3.9)

dove si e sfruttato il fatto che il vettore tangente alle curve temporali siscrive tµ = αnµ + βµ. In definitiva possiamo scrivere, scegliendo l’adattosistema di coordinate

∂tγij = −2αKij +Diβj +Djβi, (3.10)

con Di rappresentazione tri-dimensionale della derivata covariante. Questae un equazione di evoluzione per la metrica spaziale ottenuta da considerazionicinematiche.

La vera dinamica e contenuta nell’evoluzione di Kij . Per riscrivere leequazioni relativistiche dobbiamo considerare contrazioni di queste ultimecon vettore normale alle ipersuperfici e con l’operatore di proiezione.

23

3.2. Equazioni ADM

E’ necessario esprimere il tensore di Riemann quadri-dimensionaleRαβµν in termini dell’intrinseco tensore di Riemann tri-dimensionale (3)Rαβµνe della curvatura estrinseca Kµν . In questo modo potremmo arrivare ariscrivere il tensore di Einstein. L’equazione di Gauss-Codazzi ci fornisce laproiezione desiderata:

P δαPκβ P

λµP

σν Rδκλσ = (3)Rαβµν +KαµKβν −KανKβµ, (3.11)

analogamente la proiezione sulle ipersuperfici e la contrazione con unvettore normale ci fornisce la relazione di Codazzi-Mainardi

P δαPκβ P

λµn

νRδκλσ = DβKαµ −DαKβµ. (3.12)

Sfruttando le due equazioni sopra e notando che

PαµP βνRαβµν = 2nµnνGµν (3.13)

si perviene a delle equazioni che non dipendono dal tempo e dalle funzionidi gauge, ovvero, ottengo i seguenti constraints che devono essere verificati atutti i tempi e si riferiscono a una data superficie:

(3)R+K2 −KijKij = 16πρ ENERGY CONSTRAINT, (3.14)

Dj(Kij − γijK) = 8πji MOMENTUM CONSTRAINT, (3.15)

con ρ := nµnνTµν densita di energia locale misurata da un osservatoreEuleriano e ji := −P iµnνTµν densita di momento sempre osservata da unosservatore Euleriano.

L’esistenza dei constraints mostra come non sia possibile specificarearbitrariamente le 12 quantita dinamiche (ricordo che Kµν e simmetrico)come condizioni iniziali.

3.2 Equazioni ADM

I constraints visti sopra forniscono 4 delle 10 equazioni di Einstein, enon corrispondono ad equazioni dinamiche. E’ chiaro che dovranno esserele altre sei equazioni a determinare l’evoluzione del campo gravitazionale.Per trovarle proiettiamo ancora il tensore di Riemann sulle ipersuperfici,ma questa volta contraiamolo con due vettori normali in modo da ottenere

P δµPκν n

λnσRδλκσ = £nKµν +KµλKλν +

1

αDµDνα. (3.16)

In questa equazione si puo notare la comparsa di α e di una derivatalungo la direzione normale alle ipersuperfici; stiamo considerando quindiun’evoluzione nel tempo. Se adesso sfruttiamo anche la (3.11) e l’equazionedi Einstein scritta in termini di Rµν ovvero

Rµν = 8π(Tµν −

1

2gµνT

), (3.17)

24

3.3. La formulazione BSSNOK

otteniamo

£tKµν −£βKµν = −DµDνα+ α[(3)

Rµν +KKµν − 2KµλKλν

]+ 4πα[γµν(S − ρ)− 2Sµν ]

(3.18)

se come prima ci occupiamo della parte spaziale e consideriamo che £~t = ∂t

∂~tKij −£~βKij = −DiDjα+ α

[(3)Rij +KKij − 2KikK

kj

]+ 4πα[γij(S − ρ)− 2Sij ] .

(3.19)

Queste equazioni ci forniscono l’evoluzione dinamica delle sei compo-nenti indipendenti della curvatura estrinseca Kij . Le equazioni (3.19) sononote come equazioni di Arnowitt-Deser-Misner (ADM), anche se differisco-no dalla loro forma originale [14], ma sono una riformulazione dovuta aYork [15]. La differenza nasce dal fatto che il primo gruppo di equazioni ederivato passando per il tensore Gµν mentre il secondo per Rµν . Entrambii set di equazioni sono equivalenti anche se differiscono per un termineadditivo proporzionale all’energy constraint, il quale deve annullarsi peruna qualsiasi soluzione fisica. Tuttavia la differenza matematica si puo farsentire in particolare alla risposta di piccole violazioni del constraint, cherisordo contiene derivate spaziali della metrica al secondo ordine e quindicambia la natura dell’equazione differenziale. La formulazione qui usata eben posta, altro notevole vantaggio.

3.3 La formulazione BSSNOK

Data la non unicita delle equazioni di evoluzione e logico aspettarsialtri metodi di risoluzione. Una riformulazione particolarmente robusta perl’evoluzione numerica, in presenza e non, di materia, e quella conosciutacon il nome di: formulazione BSSNOK. Questa e una riformulazione delleequazioni ADM basate su di una trasformazione conmforme. Il motivodella maggiore stabilita e da ricercarsi nell’iperbolicita delle equazioniottenute. Nel nome di questa formulazione sono presenti le iniziali deinomi delle persone che ci hanno lavorato (Nakamura, Oohara per primi epoi Baumgarte, Shibata, Shapiro ). Per maggiori dettagli vedere [11] e [12].

Nella formulazione BSSNOK consideriamo un rescaling della metricaspaziale della forma

γij = ψ−4γij =⇒ ψ = γ1/12 (3.20)

con γ determinante di γij . Si chiede anche che questa relazione sia mante-nuta durante l’evoluzione. Da (3.10) si ricava

∂tγ = γ(−2αK + 2Diβi) = −2γ(αK − ∂iβi) + βi∂iγ, (3.21)

da cui∂tψ = −1

6ψ(αK − ∂iβi) + βi∂iψ. (3.22)

25

3.3. La formulazione BSSNOK

Solitamente si usa φ = lnψ = 112 lnγ cosi che γij = e−4φγij e

∂tφ = −1

6(αK − ∂iβi) + βi∂iφ. (3.23)

La seguente formulazione inoltre, separala curvatura estrinseca nellasua traccia e una parte senza traccia

Aij = kij −1

3γijK. (3.24)

si fa anche un rescaling conforme della curvatura estrinseca senza traccia

Aij = ψ−4Aij = e−4φAij . (3.25)

Il punto cruciale e che vengono introdotte tre nuove variabili ausiliarieconosciute con il nome di funzioni di connessione conformi e definite come

Γi := γjkΓijk = −∂j γij , (3.26)

dove Γijk sono i simboli di Christoffel con la metrica conforme. Ora invecedi avere 12 variabili γij e Kij usiamo 17 variabili φ, K, γij , Aij e Γi. Da-to pero che Aij e senza traccia e γij ha determinante unitario dobbiamogestire, in realta, 15 variabili dinamiche. L’evoluzione di φ e scritta sopra,mentre l’evoluzione di Aij , K e γij si ottengono direttamente dal sistemadi equazioni di evoluzione

d

dtγij = −2αAij , (3.27)

d

dtφ = −1

6αK, (3.28)

d

dtAij = e−4φ

−DiDjα+ αRij + 4πα[γij(S − ρ)− 2Sij ]

TF+ α(KAij − 2AikAkj),

(3.29)

d

dtK = −DiD

iα+ α(AijA

ij +1

3K2)

+ 4πα(ρ+ S), (3.30)

con d/dt := ∂t − £~βe dove TF denota la parte libera dalla traccia del

termine tra parentesi. Nell’evoluzione di evoluzione di K e stato usatol’energy constraint per eliminare lo scalare di Ricci.

R = KijKij −K2 + 16πρ = AijA

ij − 2

3K2 + 16πρ. (3.31)

E’ importante notare che nelle equazioni di evoluzione per Aij e K appaio-no derivate covarianti rispetto alla metrica fisica. Queste possono esserecalcolate considerando il legame tra i simboli di Christoffel fisici e conformi.L’equazione di evoluzione per Γi si ricava direttamanete da (3.10) e (3.26)

∂tΓi = −∂j

(£~βγij

)− 2(α∂jA

ij + Aij∂jα). (3.32)

Usando il momentum constraint

∂jAij = −ΓijkA

jk − 6Aij∂jφ+2

3γij∂jK + 8πji (3.33)

26

3.4. Condizioni di Gauge

con ji := e4φji, otteniamo

∂tΓi = γjk∂j∂kβ

i +1

3γij∂j∂kβ

k − 2Aij∂jα

+ 2α(

ΓijkAjk + 6Aij∂jφ−

2

3γij∂jK − 8πji

).

(3.34)

Possiamo evolvere ora liberamente le componenti di Aij e γij , conl’unica condizine A = 0. Tutto cio aumenta notevolmente la stabilita dellasimulazione.

3.4 Condizioni di Gauge

Nel formalismo 3+1 la scelta del sistema di coordinate e dato in terminidelle variabili di gauge: α lapse function e βi shift vector. Queste funzioniappaiono nelle equazioni di evoluzione (3.10) e (3.19) per la metrica ela curvatura estrinseca. Le equazioni di Einstein non ci dicono nulla ariguardo e questo e chiaro perche il sistema di coordinate puo essere sceltoliberamente. Tutta questa liberta puo essere molto vantaggiosa, possiamotrovare un sistema di coordinate che meglio si adatta al nostro problemae soluzioni con un comportamento migliore. Sorge spontanea la seguentedomanda: qual’e la scelta migliore per queste funzioni?

In generale, si scelgono α e β in modo che si adattino alle simmetriedel problema, in modo che abbiano un buon comportamento matematico euna facile implementazione numerica.

Per prima cosa e necessario legare a quantita geometriche queste fun-zioni di gauge.

Per specificare una fogliazione dello spaziotempo in ipersuperfici ditipo spazio, abbiamo bisogno di una ricetta che ci permetta di calcolare α.Tipicamente il lapse si specifica a priori come funzione dello spaziotempo,per esempio α = 1 (geodesic slicing) o con condizioni algebriche di slicingdove il lapse e definito come funzione algebrica di variabili geometriche.

Per discutere la condizione di slicing utilizzata ee conveniente consi-derare il moto di un osservatore Euleriano. In generale questo osservatorenonsara in caduta libera e quindi avra una propria accelerazione. In ter-mini della derivata direzionale delle 4-velocita dell’osservatore Euleriano,ovvero il vettore normale nµ lungo se stesso

aµ = nµ∇νnµ. (3.35)

Osservare come il vettore aµ sia ortogonale a nν , l’accelerazione e puramen-te spaziale. Possiamo legarla alle funzioni di gauge in modo che

a0 = βm∂mlnα (3.36)

ai = ∂ilnα (3.37)

Un’altra relazione importanete deriva dall’evoluzione del volume associatoad un osservatore Euleriano. Le variazioni nel tempo di questo volume

27


sono date dalla divergenza dell 4-velocita ∇µnµ che si puo legare allacurvatura estrinseca

∇µnµ = −K (3.38)

ovvero la rapidita di variazione dell’elemento di volume e semplicementel’opposto della traccia della curvatura estrinseca.

3.4.1 Maximal slicing

Come suggerito nel precedente paragrafo potremmo impostare il valoredel lapse ad una costante, ottenendo il cosiddetto geodesic slicing. Il nomederiva dal fatto che l’accelerazione propria di un osservatore Eulerianoscompare. Che implicazioni ci sono ? L’osservatore segue geodetiche ditipo tempo. Sorge un grosso problema pero: osservatori in caduta liberain un campo gravitazionale non uniforme potrebbero collidere. Otteniamouna singolarita. Da (3.19) si ottiene

∂tK − βi∂iK = −D2α+ α[KijK

ij + 4π(ρ+ S)], (3.39)

e per geodesic slicing

∂tK − βi∂iK = α[KijK

ij + 4π(ρ+ S)]. (3.40)

Il secondo membro e sempre positivo, la traccia della curvatura estrinsecalungo la direzione normale continua a crescere il che implica (3.38) che ilvolume dell’osservatore Euleriano collassa a zero. Per questo motivo lascelta di questo gauge non viene mai untilizzata se non in test numerici.

Per risolvere questo problema dobbiamo costruire una miglior condi-zione di slicing, per esempio imponendo che il volume dell’osservatoreEuleriano rimanga costante, da (3.38)

K = ∂tK = 0. (3.41)

E’ chiaro che questo deve essere richiesto sia all’inizio, che durante lasimulazione. Inoltre da (3.19) deve essere soddisfatta

D2α = α[KijK

ij + 4π(ρ+ S)], (3.42)

equazione ellittica. Questa condizione e chiamata maximal slising perche sipuo dimostrare che quando K = 0 il volume delle ipersuperfici spazialie massimale rispetto a piccole variazioni nella superficie stessa. Questogauge si comporta meglio vicino alle singolarita. Nelle vicinanze di un blackhole per esempio, il lapse collassa esponenzialmente vicino alla singolarita,tuttavia e come se l’interno del black hole fosse congelato mentre l’esternocontinua ad evolversi. Il problema e la lentezza nella risoluzione di unequazione ellittica.

3.4.2 Condizione di slicing iperbolico

Dato il grande dispendio di tempo impiegato in un equazione di tipoellittico, sono state tentate altre strade per determinare il valore di α. Sto-ricamente ci sono state due strade per la risoluzione di questo problema,

28


che hanno sfruttato equazioni di evoluzione per la funzione di lapse neltrovare condizioni di slicing.

La prima via sfrutta il noto gauge armonico delle coordinate

xα = gµν∇µ∇νxα = 0, (3.43)

che puo essere riscritto nel seguente modo espandendo l’espressione sopra

Γα := gµνΓαµν = 0. (3.44)

Usando l’espressione dei simboli di Christoffel nel formalismo 3 + 1, l’e-spressione (3.44) diventa

d

dtα ≡

(∂t −£~β

)α = −α2K, (3.45)

dove ho usato il sistema tale che x0 = t. Questa e conosciuta come condizionedi slicing armonico. Se sfruttiamo le equazioni ADM questa condizione sitraduce in

d

dtα = 0, (3.46)

con α := α/√γ. Nel caso lo shift sia nullo possiamo integrare l’ultima

espressione ed ottenere α = h(xi)√γ, con h(xi) funzione positiva arbitraria

indipendente dal tempo. La relazione persiste solo se integriamo lungo larelazione normale alle ipersuperfici, lega quindi il lapse con il l’elementodi volume che si muove con l’osservatore Euleriano.

Un altra via per utilizzare una condizione di slicing dinamica, partedal considerare la forma integrata dello slicing armonico discusso prece-dentemente α =

√γ. Tuttavia questo tipo di slicing non si comporta bene

vicino a delle singolarita. La ricerca empirica di una migliore condizionedi slicing algebrica ha condotto alla seguente forma per la funzione dilapse: α = 1 + logγ chiamata 1 + log slicing. Nella pratica si e dimostratoessere molto robusto e imita le proprieta del maximal slicing vicino allesingolarita.

Le due vie discusse sono state riunite in un lavoro di Bona e Massoche riunirono diverse condizione di slicing in un’unica famiglia [18], chenon e altro che una generalizzazione dello slicing armonico, per la quale lafunzione di lapse e scelta in modo da soddisfare la seguente equazione

d

dtα = −α2f(α)K, (3.47)

con f(α) positiva. Se f = 1 si riottiene lo slicing armonico, mentre seponiamo f = N/α, con N costante, otteniamo α = h(xi) + logγN/2. Questoe il tipo di condizione di slicing che abbiamo utilizzato nelle simulazionisvolte.

Se prendiamo un ulteriore derivata temporale nell’equazione (3.47)otteniamo

d2

dt2α = −α2f

[ ddtK − α(2f + αf

′)K2

], (3.48)

29


con f′

:= df/dα. Sfruttando anche l’equazione di evoluzione di K, (3.39),troviamo (nel vuoto) un’ equazione delle onde con un termine di sorgentequadratico in kij

d2

dt2α− α2fD2α = −α3f

[KijK

ij − (2f + αf′K2)

]. (3.49)

Questo e il motivo per cui viene chiamata condizione di slicing iperbolico.La velocita di gauge associata all’onda in una direzione specifica xi e

vg = α√fγii. (3.50)

Notare come questa velocita sia reale solo se f(α) ≥ 0, motivo dellarichiesta di f positiva.

3.4.3 Elliptic shift condition

Molto meno e conosciuto sulle condizioni di shift. La classica condizionedi shift ellittico fu proposta da Smarr e York in [16, 17]. COnsideriamodelle curve di tipo tempo congruenti con vettore tangente Zµ. L’operatoredi proiezione normale a Zµ e

Pµν = gµν + ZµZν . (3.51)

In termini di Pµν , la torsione ωµν , la deformazione θµν e l’accelerazione ζµ dellecurve congruenti sono definiti come

ωµν := Pαµ Pβν ∇[αZβ] (3.52)

θµν := Pαµ Pβν ∇(αZβ) (3.53)

ζµ := Zν∇νZµ. (3.54)

Il tensore deformazione θµν puo essere ulteriormente decomposto nellaparte senza traccia θ = ∇µZµ, conosciuto come espansione e nella partelibera σµν = θµν − (θ)/3Pµν , nota come shear. In termini di queste quantitala derivata covariante puo essere riscritta come

∇µZν = ωµν + θµν − ζµZν . (3.55)

Assumiamo che le nostre congruenze corrispondano alle linee di uni-verso degli osservatori Euleriani. Abbiamo percio Zµ = nµ, con nµ vettorenormale alle ipersuperfici di tipo spazio. Chiaramente ωµν = 0 dato che lecongruenze sono ortogonali alle ipersuperfici. Abbiamo inoltre, ζµ = aµ,con aµ accelerazione dell’osservatore Euleriano e θµν = −Kµν , con Kµν

curvatura estrinseca. Il tensore deformazione θµν = −Kµν = 1/2£~nγµνcorrisponde al moto nella direzione normale. In modo analogo potremmodefinire un tensore deformazione lungo la linea temporale

Θij :=1

2£~tγij = −αKij +

1

2£~βγij , (3.56)

in cui abbiamo considerato solo le componenti spaziali.

30


Questo nuovo tensore ci fornisce indicazioni sulla deformazione del-la forma e nel volume del volume lungo le linee temporali. Potremmoquindi provare ad usare il vettore di shift per minimizzare qualche misu-ra ottenuta con Θij , in modo da ottenere una deformazione minore nellametrica spaziale. Smarr e York proposero di minimizzare la radice nonnegativa di ΘijTheta

ij considerata globalmente sulla ipersuperfice. Se siintegra ΘijTheta

ij su di uno slice, la variazione rispetto al vettore di shiftsi traduce nella seguente equazione

D2βi +DiDjβj +Rijβ

j = 2Dj(αKij), (3.57)

con il tensore di Ricci che compare quando vengono commutate le derivatecovarianti. Questa equazione e conosciuta come shift di minima deformazionee ci fornisce tre equazioni ellittiche per trovare le componenti dello shift.Dato che il cambio di volume e legato alla traccia della curvatura estrin-seca, potremmo chiederci se non sia meglio usare il vettore di shift perminimizzare il cambio nella forma del volume senza occuparci della taglia.Definiamo lo shear associato alle linee temporali come la parte a traccialibera di Ωij :

Σij := Θij −1

3γijΘ. (3.58)

Smarr e York chiamano Σij tensore di distorsione. Dalla definizione sopra sitrova

Σij = −αAij +1

2(Lβ)ij , (3.59)

doveAij = Kij−(γij/3)K e la parte senza traccia della curvatura estrinsecae (Lβ)ij e la forma conforme di killing associata con il vettore di shift [2]:

(Lβ)ij := Diβj +Djβi −2

3γijDkβ

k. (3.60)

Si puo anche riscrivere il tensore di distorsione come

Σij =1

2γ1/3£~tγij , (3.61)

con γ = γ−1/3γij metrica conforme.Se ora minimizzassimo l’integrale di ΣijΣij sulla ipersuperficie rispetto

al vettore di shift troveremmo la seguente condizione

DjΣij = 0, (3.62)

cje implica∆Lβ

i = 2Dj(αAij), (3.63)

dove l’operatore ∆L e definito come

∆Lβi := Dj(Lβ)ij = D2βi +DjD

iβj − 2

3DiDjβ

j

= D2βi +1

3DiDjβ

j +Rijβj .

(3.64)

31


L’equazione (3.61) puo essere riscritta

Dj(∂tγij) = 0. (3.65)

Questa e molto simile alla condizione proposta da Dirac in un suo lavoroper trovare un analogo della radiazione di gauge nella relativita generale[20, 17]

∂t∂j γij = ∂j∂tγij = 0. (3.66)

Dato che le funzioni di connessione conforme sono date da Γi = −∂j γij , lacondizione di gauge sopra si traduce in

∂tΓi = 0. (3.67)

Questa condizione e conosciuta come gamma freezing dato che congela tregradi di liberta indipendenti. E’ spesso usata nel contesto della formula-zione BSSNOK. Si puo dimostrare la condizione di minima distorsione eil gamma freezing sono equivalenti a meno di termini che coinvolgonoderivate prime spaziali della metrica conforme e del fattore conforme.

3.4.4 Driver Conditions

Come gia sottolineato in precedenza risolvere equazioni ellittiche edispendioso in termini di tempo computazionale. Un’idea per aggirarequesto problema e cercare di sostituire le nostre equazioni per la funzionedi shift con altre di piu facile risoluzione, per esempio di tipo parabolico,risolte in un tempo fittizio ed evolvendo queste ultime equazioni finche lasoluzione non rilassa a quella cercata. Consideriamo per esempio l’equa-zione di Laplace con date condizioni al contorno (Dirichlet o Newmannnon importa)

D2φ = 0. (3.68)

Potremmo rimpiazzarla con un’ altra equazione

∂tφ = εD2φ, (3.69)

dove τ e il tempo fittizio, mentre ε > 0 e una qualche costante che fissa lascala temporale del rilassamento. Potremmo evolvere φ con dati inizialiarbitrari fino a quando questo campo non rilassa ad una soluzione stazio-naria che soddisfa le condizioni iniziali. Naturalmente dovremo porre ,se cerchiamo una situazione stazionaria ∂tφ = 0. Purtroppo le condizionidi stabilita numerica per un equazione parabolica, come quella del calore,richiedono 4τ . ε(4x)2, cosi che se si raffina troppo il reticolo il passotemporale diventa proibitivamente piccolo. Potremmo pensare di usare unequazione tipo onde

∂2τφ− v2D2φ = 0, (3.70)

per la quale le condizioni di stabilita numerica sono meno restrittive,v4τ . 4x, con v velocita caratteristica del sistema. E’ necessaria cautelapero, perche quest’ultimo approccio non rilassa con arbitrarie condizioni alcontorno, l’onda si riflettera sui bordi.

32


E’ fondamentale ricordare pero, che noi stiamo risolvendo un’equazioneper una funzione di gauge, arbitraria, quindi non e importante trovareuna soluzione esatta, anche una approssimata puo funzionare altrettantobene. Possiamo anche usare un tempo reale invece di uno fittizio. Questiragionamenti conducono a quella che viene detta driver conditions, dellequali la prima e la versione parabolica della distorsione minimale “guidata”

∂tβi = ε

[D2βi +

1

3DiDjβ

j +Rijβj − 2Dj(αA

ij)], (3.71)

dove ε > 0 e scelto piccolo a sufficenza da garantire la stabilita numerica.Da osservare che la versione parabolica non permette una risposta

celere dello shift ai cambi di geometria, a meno di utilizzare ε grande.Per questo motivo si preferisce la versione iperbolica, hyperbolic driver

∂2t β

i = α2ξ[D2βi +

1

3DiDjβ

j +Rijβj − 2Dj(αA

ij)], (3.72)

dove ξ > 0 e un’ arbitraria funzione positiva della posizione e/o di α,che controlla la velocita dell’onda nelle differenti direzioni. Solitamentesi utilizzano insieme a questa equazione condizioni al contorno per ondeuscenti, che riducono al minimo la riflessione, ovvero permettono l’uscitadell’onda dal dominio computazionale.

33

CAPITOLO 4

IL SOFTWARE

Indice4.1 CACTUS e Einstein Toolkit . . . . . . . . . . . . . . . . 354.2 Runge-Kutta e MOL . . . . . . . . . . . . . . . . . . . . . 51

Etiam et etiam iterare oportet.

4.1 CACTUS e Einstein Toolkit

Cactus, [27], e un ambiente open source dedicato al problem solving,sviluppato per scienziati ed ingegneri. La sua struttura modulare si adattafacilmente al calcolo parallelo, e portabile e sostenuta da differenti gruppi.

Il nome Cactus deriva dalla struttura di questo ambiente, un nucleocentrale al quale si possono connettere le spine (“Thorns”), le quali rap-presentano estensioni dell’ambiente. I Thorns possono implementare ap-plicazioni scientifiche personalizzate, come per esempio fluidodinamicacomputazionale.

L’Einstein Toolkit [29] consiste di un set di piu di 100 Thorns di Cactusper la relativita numerica, con tools per gestire la simulazione e la visualiz-zazione. Questo toolkit comprende codice per la risoluzione delle equazionidi Einstein nel vuoto, o con materia. La risoluzione dell’idrodinamica relati-vistica avviente tramite il Whisky code [30] . Tramite questo tool sono staterisolte le equazioni di Einstein nella formulazione di BSSNOK ed estrat-te, tramite altri thorns, le informazioni sulle onde gravitazionali. Il thornWaveExtract si occupa dell’estrazione dell’informazione delle onde gravi-tazionali tramite il metodo dello sviluppo in multipoli, mentre Psikadeliasfrutta gli scalari di Weyl. Altri thorns si occupano della lettura e scritturain parallelo dell’input e dell’output sia in formato ASCII che in HDF5. Non

35

4.1. CACTUS e Einstein Toolkit

meno importante e CARPET, che permette di generare la griglia, ed i suoiraffinamenti, per la risoluzione delle equazioni alle derivate parziali.

Figura 4.1: Schema di funzionamento dell’Einstein Toolkit

Nel seguito discutero brevemente i thorns sfruttati prendendo in consi-derazione un file di parametri utilizzato per eseguire la simulazione.

La prima cosa da fare e indicare quali thorns devono essere attivatitramite l’opzione ActiveThorns. Nel listato e mostrato l’elenco dei thornsche abbiamo sfruttato e che descrivero brevemente. Per ogni thorn attivatoe possibile settarne le funzionalita opzionali.

CartGrid3D permette di settare le coordinate su di una griglia carte-siana a tre dimensioni in modo flessibile. E’ possibile la scelta di dominidifferenti da utilizzare sulla griglia, per esempio si possono sfruttare gliottanti, in modo da poter usufruire delle simmetrie differenti del problemaoggetto di studio e facilitare l’esecuzione della simulazione. Questo thornsprevede anche routin per la registrazione di simmetrie delle funzioni defini-te sulla griglia e di applicare condizioni di simmetria attraverso gli assi dellagriglia stessa. CartGrid3D::avoid_origin="no" non permette di ave-re un punto sulla griglia esattamente uguale a zero. La griglia viene cioetraslata di mezzo passo dall’origine. CartGrid3D::type="coordbase"specifica come inizializzare la griglia. In questo caso diciamo al program-ma che per fissare la dimensione della griglia, la sua risoluzione e ne-cessario fare riferimento ad un altro thorn, nello specifico CoordBase.Questa applicazione fornisce meccanismi per la registrazione del sistemadi coordinate, mantenendo un database di sistemi di coordinate usandotabelle chiave-valore. Come sopra citato, permetta la scelta delle dimen-sioni della griglia, tramite CoordBase::xmax...xmin e le opzioni ana-loghe, della sua discretizzazione tramite CoordBase::dx..., che comesi vedra nel seguito potra essere ulteriormente raffinata. Con queste im-postazioni abbiamo inizializzato un ottante. Le altre opzioni di questothorn hanno la funzionalita di estendere superiormente ed inferiormentela griglia discretizzata in modo da poter sfruttare l’algoritmo di calcolodi quantita sul reticolo, per le quali serve conoscere non solo il punto sulquale vengono calcolate, ma anche i primi vicini. Il dominio fisico quin-di non coincide necessariamente con quello computazionale che usual-

36


Listing 4.1: File dei parametri, attivazione dei thorns1 #=========================================================

#3 # Active Thorns

#5 #=========================================================

7

ActiveThorns = ”Time MoL Coordbase CartGrid3d9 Boundary StaticConformal

SymBase ADMBase TmunuBase11 HydroBase InitBase ADMCoupling

ADMMacros IOUtil Formaline13 SpaceMask CoordGauge Constants

LocalReduce aeilocalinterp LoopControl15 Carpet CarpetLib CarpetReduce

CarpetRegrid2 CarpetInterp17 CarpetIOASCII CarpetIOScalar

CarpetIOHDF5 CarpetIOBasic19 ReflectionSymmetry ADMConstraints

NaNChecker EOS Base EOS IdealFluid21 EOS Polytrope SphericalSurface

GenericFD NewRad SummationByParts23 ML BSSN ML BSSN Helper

Whisky25 Whisky TOVSolverC

WeylScal427 PsiKadelia

ManualTermination29 WaveExtract”

mente e allargato. CoordBase::boundary_shiftout_x_lower=1 eCoordBase::boundary_size_x_lower=3 nel listato, dicono appun-to che devo considerare tre punti aggiuntivi al reticolo nella parte inferioree questi boundary devono essere pensati fuori dal dominio fisico. Per l’as-segnazione dei valori a questi punti si fara riferimento alle simmetrie diriflessione attorno agli assi, le quali sono attivate in un thorn che vedremonel seguito.

37


Listing 4.2: File dei parametri, CartGrid e CoordBase. Come creare la griglia30 #==================================

# CartGrid3D32 #==================================

34 CartGrid3D::type = ”coordbase”CartGrid3D::avoid origin = ”no”

36

#==================================38 # CoordBase

#==================================40

42 CoordBase::boundary size x upper = 3CoordBase::boundary size y upper = 3

44 CoordBase::boundary size z upper = 3CoordBase::boundary shiftout x upper = 0

46 CoordBase::boundary shiftout y upper = 0CoordBase::boundary shiftout z upper = 0

48 CoordBase::boundary shiftout x lower = 1CoordBase::boundary shiftout y lower = 1

50 CoordBase::boundary shiftout z lower = 1CoordBase::boundary size x lower = 3

52 CoordBase::boundary size y lower = 3CoordBase::boundary size z lower = 3

54 CoordBase::domainsize = ”minmax”CoordBase::dx = 4

56 CoordBase::dy = 4CoordBase::dz = 4

58 CoordBase::xmax = 120CoordBase::xmin = 0

60 CoordBase::ymax = 120CoordBase::ymin = 0

62 CoordBase::zmax = 120CoordBase::zmin = 0

38


Listing 4.3: File dei parametri, Riflessioni, come sfruttare le simmetrie64 #==================================

# ReflectionSymmetry66 #==================================

68

ReflectionSymmetry::reflection x = ”yes”70 ReflectionSymmetry::reflection y = ”yes”

ReflectionSymmetry::reflection z = ”yes”72 ReflectionSymmetry::avoid origin x = ”no”

ReflectionSymmetry::avoid origin y = ”no”74 ReflectionSymmetry::avoid origin z = ”no”

Listing 4.4: File dei parametri, come cacolare Tµν#==================================

76 # TmunuBase#==================================

78

TmunuBase::support old CalcTmunu mechanism = ”no”80 TmunuBase::stress energy storage = yes

TmunuBase::stress energy at RHS = yes82 TmunuBase::timelevels = 1

TmunuBase::prolongation type = none

ReflectionSymmetry permette di attivare le simmetrie, in questocaso riflessioni come citato sopra.

Il thorn TmunuBase fornisce funzioni definite sulla griglia che per-mettono di calcolare il tensore energia impulso Tµν e memorizza quandoqueste funzioni sono state calcolate. Tutto cio permette a differenti thorns dicooperare senza esplicitamente dipendere gli uni dagli altri. TmunuBase eimportante per il calcolo del tensore Tµν quanto ADMBase e importante peril calcolo del tensore metrico. TmunuBase::prolongation_type=nonee l’opzione necessaria nel caso, come questo, siano indicati pochi livellitemporali. TmunuBase::stress_energy_at_RHS=yes specifica che iltensore in questione si trova al lato destro delle equazioni di Einstein.

Il thorn ADMMacros fornisce Macro che possono essere utilizzate per lacomputazione dei simboli di Christoffel o per le componenti del tensore diRiemann usando le variabili base del thorn ADMBase. Le macro possonoessere utilizzate sia nel codice Fortran sia nel codice C. Lavorano puntual-mente per il calcolo di quantita sulla griglia in un punto specifico (i,j,k).ADMMacros::spatial_order=4 definisce una differenziazione spazialedel quarto ordine.

ADMBase e una parte fondamentale del software, fornisce il cuore del-l’infrastruttura per i thorns che implementano la relativita generale su diuna griglia 3D nel formalismo 3+1. Grazie alle sue variabili fondamentali(3-metrica, curvatura estrinseca, lapse, vettore di shift...) per il formalismo3+1 e ad una serie di ulteriori parametri, permette di regolare l’evoluzionedelle equazioni del sistema. Queste variabili sono utilizzate per comuni-care tra i vari thorns che forniscono dati iniziali, metodi di evoluzione e

39


Listing 4.5: File dei parametri, ADMBase, ADMConstraints, ADMMacros84 #==================================

# ADMMacros86 #==================================

88 ADMMacros::spatial order = 4

90 #==================================# ADMBase

92 #==================================

94 ADMBase::metric type = ”physical”ADMBase::evolution method = ”ML BSSN”

96 ADMBase::lapse evolution method = ”ML BSSN”ADMBase::shift evolution method = ”ML BSSN”

98 ADMBase::dtlapse evolution method= ”ML BSSN”ADMBase::dtshift evolution method= ”ML BSSN”

100 ADMBase::initial data = ”tov”ADMBase::initial lapse = ”tov”

102 ADMBase::initial shift = ”tov”ADMBase::initial dtlapse = ”zero”

104 ADMBase::initial dtshift = ”zero”

106

#==================================108 # ADMConstraints

#==================================110

ADMConstraints::bound = ”static”112 ADMConstraints::constraints timelevels = 3

ADMConstraints::constraints persist = yes

routines di analisi. ADMBase::metric_type="physical" dichiara chela metrica e la curvatura estrinseca sono quelle fisiche che fornisce questothorn corrispondono a quelle fisiche. Si potrebbe richiedere anche quellaconforme. Come si puo notare dal listato che contiene i parametri ADM, ilmetodo di evoluzione utilizzato e ML BSSN, ovvero si sfrutta il formalismoBSSN con il metodo delle linee per il calcolo delle equazioni differenziali.

ADMConstraints permette il calcolo del constraints energia e impulso.ADMConstraints::bound="static" sostanzialmente dice di non farenulla di speciale al bordo e ADMConstraints::con...timelevels=3che il calcolo deve avvenire ogni tre livelli temporali.

40


Listing 4.6: File dei parametri, SpaceMask114 #==================================

# SpaceMask116 #==================================

118 SpaceMask::use mask = ”yes”

Listing 4.7: File dei parametri, Cactus#==================================

120 # Cactus#==================================

122

Cactus::cctk initial time = 0124 Cactus::terminate = ”time”

Cactus::cctk final time = 100000

SpaceMask fornisce funzioni sulla griglia che possono essere sfruttateper immagazzinare stati definiti dall’utente ad ogni punto della griglia.Un programmatore potrebbe voler sapere se un punto della griglia siasettato allo stato “interno”, “eccitato”, “bordo”. E’necessario includerlo perl’utilizzo di altri thorns.

Le variabili di Cactus fissano il tempo di partenza e di fine dellasimulazione. Come si vedra nel seguito e stata impostata anche una termi-nazione manuale, nel caso la simulazione si rivelasse troppo lunga si fer-merebbe dopo il numero di ore indicato nel thorn ManualTermination.Cactus::terminate="time" specifica che il controllo deve essere fattosul tempo e non sul numero di iterazioni, parametro che potrebbe essereimpostato al posto di questi per indicare il termine della simulazione.

Carpet e una libreria scritta in C++ che permette il raffinamento dellagriglia di calcolo e puo essere utilizzato al posto di un’altra libreria chiamataPUGH. Carpet supporta multi-livelli di raffinamento e multi-patches di gri-glia. Puo girare in parallelo anche se ancora non molto efficentemente. Nellasimulazione eseguita con questo file di parametri anche se il livello massimodi raffinamenti e fissato a dieci Carpet::max_refinement_levels=10abbiamo utilizzato cinque livelli CarpetRegrid2::num_levels_1=5,

Listing 4.8: File dei parametri126 #==================================

# NaNChecker128 #==================================

130 NaNChecker::check after = 6000NaNChecker::check every = 1

132 NaNChecker::action if found = ”just warn”NaNChecker::check vars = ”ADMBase::metric ADMBase::lapse

134 ADMBase::shift HydroBase::rhoHydroBase::eps HydroBase::press

136 HydroBase::vel admbase::gxx”

41


Listing 4.9: File dei parametri, Carpet#==================================

138 # Carpet#==================================

140

142 Carpet::convergence level = 0Carpet::domain from coordbase = ”yes”

144 Carpet::enable all storage = ”no”Carpet::ghost size = 3

146 Carpet::init each timelevel = ”no”Carpet::max refinement levels = 10

148 Carpet::num integrator substeps = 4Carpet::prolongation order space = 3

150 Carpet::prolongation order time = 2Carpet::regrid during initialisation = ”no”

152 Carpet::regrid during recovery = ”no”Carpet::regrid in level mode = ”yes”

154 Carpet::use buffer zones = ”yes”Carpet::poison new timelevels = ”yes”

156 Carpet::check for poison = ”no”Carpet::poison value = 113

158 Carpet::init 3 timelevels = noCarpet::init fill timelevels = ”yes”

160

CarpetLib::poison new memory = ”yes”162 CarpetLib::poison value = 114

164 CarpetRegrid2::regrid every = 0CarpetRegrid2::num centres = 1

166 CarpetRegrid2::num levels 1 = 5CarpetRegrid2::radius 1[1] = 60.0

168 CarpetRegrid2::radius 1[2] = 30.0CarpetRegrid2::radius 1[3] = 15.0

170 CarpetRegrid2::radius 1[4] = 8.0

che cominciano in punti prefissati con l’opzioneCarpetRegrid2::radius_...sono tutti centrati sullo stesso punto CarpetRegrid2::num_centres=1e calcolati solo una volta, durante l’inizializzazione dei dati inizialiCarpetRegrid2::regrid_every=0 . E’ inoltre specificato che stiamoutilizzando il thorn coordbase. Carpet::ghost_siz=3 indica quantodeve essere grande la zona fantasma, necessaria per la paralizzazione.Quando vengono caricate parti differenti della griglia in parallelo e necessa-rio conoscere il bordo di tale griglia, la cui taglia varia in base all’algoritmousato.

42


Listing 4.10: File dei parametri, Method of Line#==================================

172 # MoL#==================================

174 MoL::ODE Method = ”rk4”MoL::MoL Intermediate Steps = 4

176 MoL::MoL Num Scratch Levels = 1

Figura 4.2: Come appare una griglia con diversi livelli di raffinamento nella partecentrale.

Tramite il thorn MoL fissiamo i parametri da usare per risolvere le equa-zioni differenziali a derivate parziali con il metodo delle linee, nel caso,come il nostro, ci siano thorn che sfruttino questo metodo. Nel sottopa-ragrafo successivo spieghero la teoria che sta alla base di questo metodo.MoL::MoL_Intermediate_Steps=4 controlla l’accuratezza del metodose l’algoritmo usato e il Runge-Kutta MoL::ODE_Method="rk4". In que-sto caso si sfrutta il Runge-Kutta al quarto ordine.MoL::MoL_Num_Scratch_Levels=1controlla l’ammontare di spazio utilizzato per il calcolo, se insufficenteinvia un errore.

ML_BSSN e il metodo di evoluzione scelto nel thorn ADMBase. Fissia-mo l’evoluzione del lapse ML_BSSN::harmonicF=1.0 che significa averscelto il metodo di evoluzione 1+log e imostiamo un controllo sul limitedei valori dei parametri. Piu che altro sono utili questi ultimi parametriquando si vuole simulare un black hole.

43


Listing 4.11: File dei parametri, metodo di evoluzione#==================================

178 # ML BSSN#==================================

180

ML BSSN::timelevels = 3182 ML BSSN::harmonicN = 1 # 1+log

ML BSSN::LapseACoeff = 1.0184 ML BSSN::ShiftBCoeff = 1.0

ML BSSN::ShiftGammaCoeff = 0.0186 ML BSSN::AlphaDriver = 0.0

ML BSSN::BetaDriver = 0.0188 ML BSSN::LapseAdvectionCoeff = 0.0

ML BSSN::ShiftAdvectionCoeff = 0.0190 ML BSSN::MinimumLapse = 1.0e−8

ML BSSN::my initial boundary condition = ”extrapolate−gammas”192 ML BSSN::my rhs boundary condition = ”NewRad”

ML BSSN::ML log confac bound = ”none”194 ML BSSN::ML metric bound = ”none”

ML BSSN::ML Gamma bound = ”none”196 ML BSSN::ML trace curv bound = ”none”

ML BSSN::ML curv bound = ”none”198 ML BSSN::ML lapse bound = ”none”

ML BSSN::ML dtlapse bound = ”none”200 ML BSSN::ML shift bound = ”none”

ML BSSN::ML dtshift bound = ”none”

44



# InitBase204 #==================================

InitBase::initial data setup method = ”init all levels”

Listing 4.13: File dei parametri, IO206 #==================================

# IO208 #==================================

210 IOBasic::outInfo every = 1IOBasic::outInfo vars = ”HydroBase::rho ADMBase::lapse”

212

IO::out dir = data1 dx4/IOScalar/214 IOScalar::outScalar every = 32

IOScalar::one file per group = yes216 IOScalar::outScalar vars = ”

HydroBase::rho218 HydroBase::press

HydroBase::eps220 HydroBase::vel

ADMBase::lapse222 ADMBase::metric

ADMBase::curv224 ADMConstraints::ham

ADMConstraints::momentum226

”

Il thorn InitBase specifica come i dati devono essere inizializzati.Non imposta il valore dei dati, ma si comporta da conveniente repositoryche ricorda come i dati devono essere settati. Seleziona il meccanismo concui attribuire valori ai dati.

Nel seguito vengono presentati alcuni thorns di input output. Il fun-zionamento e analogo per tutti questi programmi, siano essi applicazioniche trattano dati ascii o hdf5. Si deve decidere ogni quanto far scrivere idati, quali titi di dato e dove. Nel gruppo che specifica l’IO per i dati HDF5e possibile chiedere quale dei livelli di raffinamento si desidera stampa-re. CarpetIOScalar immagazzina dati sommari, come il massimo ed ilminimo di una data variabile, per esempio della densita. Gli altri thornsdanno i valori delle variabili sulla griglia.

HydroBase estende il framework dell’EisteinToolkit fornendo un inter-faccia che permette l’implementazione della magnetoidrodinamica. Datoche la fluidodinamica relativistica che abbiamo utilizzato e quella imple-mentata in Whisky abbiamo dovuto specificareHydroBase::evolution_method="Whisky".

Il thorn Whisky fornisce il codice per gestire l’idrodinamica relativisticain tre dimensioni. Dobbiamo indicare l’equazione di stato da sfruttare ed ilmetodo per il calcolo del tensore di Riemann. Whisky_TOVSolverC si oc-

45


Listing 4.14: File dei parametri, IO228 IOASCII::out1D dir = data1 dx4/IO 1D/

IOASCII::out1D every = 20230 IOASCII::one file per group = yes

IOASCII::output symmetry points = no232 IOASCII::out3D ghosts = no

IOASCII::out3D outer ghosts = no234 IOASCII::out1D vars = ”

HydroBase::rho236 HydroBase::press

HydroBase::eps238 HydroBase::vel

ADMBase::lapse240 ADMBase::metric

ADMBase::curv242 ADMConstraints::ham

ADMConstraints::momentum244 PsiKadelia::WeylComponents

PsiKadelia::IJinvariants246 ”

Listing 4.15: File dei parametri, IO#==================================

248 # CarpetIOScalar#==================================

250

CarpetIOScalar::outScalar dir = ”data1 dx4/DataCARPET HydroAnalysis”252 CarpetIOScalar::outScalar reductions = ”minimum maximum

norm1 norm2254 norm inf”

Listing 4.16: File dei parametri, IO#==================================

256 # CarpetIOHDF5#==================================

258

#CarpetIOHDF5::checkpoint = ”yes”260 CarpetIOHDF5::open one input file at a time = ”yes”

CarpetIOHDF5::out dir = ”data1 dx4/DataCARPET H5”262 CarpetIOHDF5::out every = 512

CarpetIOHDF5::use grid structure from checkpoint = ”yes”264 iohdf5::out vars = ”

hydrobase::rho compression level=9 refinement levels=4266 hydrobase::presscompression level=9 refinement levels=4

hydrobase::eps compression level=9 refinement levels=4268 hydrobase::vel compression level=9 refinement levels=4

ADMBase::lapse compression level=9 refinement levels=4270 ADMBase::shift compression level=9 refinement levels=4

ADMBase::curv compression level=9 refinement levels=4272 ADMBase::metriccompression level=9 refinement levels=4 ”

46


Listing 4.17: File dei parametri#==================================

274 # HydroBase#==================================

276

HydroBase::timelevels = 3278 HydroBase::evolution method = ”Whisky”

Listing 4.18: File dei parametri, Whisky#==================================

280 # Whisky#==================================

282

Whisky::riemann solver = ”Marquina”284 Whisky::Whisky eos type = ”Polytype”

Whisky::Whisky eos table = ”2D Polytrope”286 Whisky::recon method = ”ppm”

Whisky::Whisky stencil = 3288 Whisky::bound = ”none”

Whisky::rho abs min = 1.e−10290

Whisky TOVSolverC::TOV Rho Central[0] = 1.28e−3292 Whisky TOVSolverC::TOV Gamma[0] = 2.0

Whisky TOVSolverC::TOV K[0] = 100.0

cupa della risoluzione delle equazioni di Tolman-Oppenheimer-Volkof perl’equilibrio idrostatico di una stella sferica non rotante. Abbiamo impostatocome suggerito precedentemente i valori Γ = 2 e k = 100 nell’equazione distato. Whisky::riemann_solver="Marquina".

Time calcola il time step per l’evoluzione secondo la regola seguente4t = dtfac ∗min(4xi).

Il thorn WaveExtract calcola le Ψeven e le Ψodd, informazioni che per-mettono di estrarre informazioni sulle onde gravitazionali. L’idea di fondoe fissare dei detectors, ovvero dei punti di osservazione sui quali si pre-sume che la metrica sia circa quella di Schwarzschild. Calcolare il raggiodi Schwarzschild per calcolare la corrispondente metrica ed utilizzarequest’ultima come quantita da sottrarre alla metrica totale calcolata dalprogramma. Otteniamo cio la perturbazione hµν . Tramite questi valori sicalcolano i vari coefficenti dello sviluppo delle armoniche sferiche, de-scritte nel primo capitolo, ed infine si calcolano le funzioni di Zerilli eMoncrief. I dati verrano presi ogni quattro iterazioni saranno immagazzi-nati in un’apposita directory da noi scelta. WaveExtract::l_mode=6 e


# SummationByParts296 #==================================

298 SummationByParts::order = 4

47


Listing 4.20: File dei parametri, Time#==================================

300 # Time#==================================

302

Time::dtfac = 0.25


# PsiKadelia306 #==================================

308 PsiKadelia::psif vec = ”standard−radial”PsiKadelia::psikadelia persists = ”yes”

310 PsiKadelia::ricci prolongation type = ”none”PsiKadelia::weyl timelevels = 3

WaveExtract::m_mode=6 limitano i modi che devonmo essere calcolati.

Gli ultimi due listati mostrano come sia possibile definire dei checkpointin modo da poter proseguire una simulazione, sia che sia stata interrotta, siache abbia finito le iterazioni impostate e come interrompere la simulazionetramite timer. Si po scegliere la directory in cui salvare i checkpoint per poipoterli recuperare IO::checkpoint_dir="data1_dx4/CHECKPOINT"e ogni quanto sono fissati IO::checkpoint_every=2048. Il file cheviene restituito e di tipo HDF5 IOHDF5::checkpoint="yes"

Grazie al thorn Termination fissiamo un orario limite prima dellaterminazione della simulazione. In questo caso abbiamo fatto durare lasimulazione dieci giorni.

48



# WaveExtract314 #==================================

316

WaveExtract::detector radius[0] = 30318 WaveExtract::detector radius[1] = 35









WaveExtract::interpolation operator = ”Lagrange polynomial interpolation”336 WaveExtract::interpolation order = 2

WaveExtract::l mode = 6338 WaveExtract::m mode = 6

WaveExtract::maximum detector number = 18340 WaveExtract::out dir = ”data1 dx4/WaveExtraction”

WaveExtract::out every = 4342 WaveExtract::rsch2 computation = ”average Schwarzschild metric”

WaveExtract::switch output format = 21344 WaveExtract::verbose = 1

Listing 4.23: File dei parametri, Checkpoint#==================================

346 # Checkpoint parameters#==================================

348

IO::checkpoint dir = ”data1 dx4/CHECKPOINT”350 IO::recover dir = ”data1 dx4/CHECKPOINT”

IO::checkpoint every = 2048352 IO::checkpoint keep = 2

IO::recover = ”autoprobe”354 IO::checkpoint on terminate = ”yes”

IO::recover file = ”checkpoint.chkpt”356 IOHDF5::checkpoint = ”yes”

IOHDF5::use reflevels from checkpoint = ”yes”

49


Listing 4.24: File dei parametri, Terminazione Manuale358 #==================================

# Termination360 #==================================

362 ManualTermination::on remaining walltime = 1 # minutes before terminationManualTermination::max walltime = 241 # hours

364 ManualTermination::termination from file = yesManualTermination::termination file = ”terminate job”

366 ManualTermination::check file every = 16 # evolution steps

50

4.2. Runge-Kutta e MOL

4.2 Runge-Kutta e MOL

La formula del metodo di Eulero per la risoluzione di equazioni diffe-renziali e

yn+1 = yn + hf(xn, yn), (4.1)

che permette l’ avanzamento di una soluzione da xn a xn+1 ≡ xn + h.Questa formula non e simmetrica, avanza la soluzione di un intervallo h,ma usa l’informazione della derivata solo in un verso. Questo significa chel’errore commesso e dell’ordine di h.

Questo metodo non e accurato e nemmeno molto stabile.Se invece consideriamo uno step intermedio, sia nel passo x che y,

otteniamo il seguente schema

k1 = hf(xn, yn) (4.2)

k2 = hf(xn +1

2h, yn +

1

2k1) (4.3)

yn+1 = yn + k2 +O(h3) (4.4)

il cui errore va come il cubo di h. Il metodo sopra e chiamato metodo diRunge-Kutta del secondo ordine. Ci sono altri modi per valutare il secondomembro f(x, y), di (4.1), tutti al primo ordine, ma che hanno coefficientiche portano un errore di ordine piu elevato. Uno schema molto preciso, diordine O(h5) e il cosiddetto schema di Runge-Kutta del quarto ordine

k1 = hf(xn, yn) (4.5)

k2 = hf(xn +1

2h, yn +

1

2k1) (4.6)

k3 = hf(xn +1

2h, yn +

1

2k2) (4.7)

k4 = hf(xn + h, yn + k3) (4.8)

yn+1 = yn +k1

6+k2

3+k3

3+k4

6+O(h5). (4.9)

In questo schema si valuta l’evoluzione passando per quattro step.MOL (method of line) e un particolare approccio per la risoluzione delle

equazioni differenziali. Si comincia prima discretizzando le dimensionispaziali lasciando la dimensione temporale continua. Il vantaggio di questometodo sta nella possibilita di disaccoppiare la scelta della differenziazionespaziale e temporale. Permette di accoppiare codici differenti sviluppatiper la risoluzione di differenti equazioni, motivo per cui viene utilizzatospesso in relativita numerica, si riescono ad implementare facilmente leequazioni della fluidodinamica con le equazioni di evoluzione 3 + 1. Ingenerale, se abbiamo un equazione del tipo

∂tu = S(u), (4.10)

con S un qualche operatore differenziale e utilizziamo le differenze finitestandard per risolvere la parte spaziale, ci riduciamo ad un sistema diequazioni del tipo

dudt

= Su, (4.11)

51

4.2. Runge-Kutta e MOL

con u vettore costruito con i valori di u sui differenti punti della grigliaspaziale, mentre S e una matrice tipicamente sparsa che accoppia i varipunti sulla griglia. A questo punto possiamo usare un altro schema per larisoluzione della parte temporale, per esempio il RK4.

52

CAPITOLO 5

SIMULAZIONI

Indice5.1 Introduzione . . . . . . . . . . . . . . . . . . . . . . . . 535.2 Il sistema fisico e la sua evoluzione . . . . . . . . . . . 545.3 CACTUS-CARPET-WHISKY:codice non-lineare 3D per

la relativita generale . . . . . . . . . . . . . . . . . . . . 545.4 Gauges . . . . . . . . . . . . . . . . . . . . . . . . . . . . 565.5 Evoluzione della materia . . . . . . . . . . . . . . . . . 565.6 Estrazione delle onde con lo scalare di Newmann-Penrose

Ψ4 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 575.7 Estrazione delle onde tramite multipoli . . . . . . . . . 585.8 I dati iniziali . . . . . . . . . . . . . . . . . . . . . . . . . 595.9 Simulazioni . . . . . . . . . . . . . . . . . . . . . . . . . 60

La cosa piu incomprensibile dell’Universo e che esso siacomprensibile.

ALBERT EINSTEIN

5.1 Introduzione

Lo studio e proseguito con l’esecuzione di simulazioni svolte utilizzan-do il sistema di calcolo Cactus/Carpet/Whisky per la computazione in3+1 dimensioni dell’equazioni di Einstein accoppiate ad un fluido perfetto.Dopo aver compreso il funzionamento di questo codice [27, 28, 29, 30]e stato possibile analizzare il comportamento di un fenomeno non fisico,che si evidenzia durante l’estrazione del segnale delle onde, dovuto a unadiscretizzazione dello spazio tempo, necessariamente presente in calcoli sureticolo. Per questo scopo abbiamo evoluto un modello di stella di neutronipolitropico, oscillante, non rotante, in tre dimensioni in relativita generale.

L’estrazione delle onde e avvenuta attraverso due differenti approcci:l’identificazione delle perturbazioni dello spazio tempo di Schwarzschild

53

5.2. Il sistema fisico e la sua evoluzione

con la decomposizione in multipoli, che e stato il metodo tradizionale,sviluppato originariamente da Regge e Weeler [4], Zerilli [5, 6], Moncrief[9] e l’estrazione dei segnali dell’informazione di onde gravitazionali intermini delle componenti del tensore di curvatura di Weyl rispetto aduna tetrade nulla, usando quello che viene conosciuto come formalismo diNewman-Penrose, molto in uso negli anni recenti.

5.2 Il sistema fisico e la sua evoluzione

Le simulazioni usano l’intero set delle equazioni di Einstein

Gµν = 8πTµν (5.1)

accoppiate con un fluido perfetto di materia, il cui tensore energia-impulsoe dato da

Tµν = ρ(

1 + ε+p

ρ

)uµuν + pgµν , (5.2)

dove uµ e la quadri-velocita del fluido, p la pressione del fluido, ε e l’energiaspecifica interna (energia per unita di massa) e ρ e la densita della massaa riposo. In questi termini e = ρ(1 + ε) e la densita di energianel sistemadi riferimento a riposo, rispetto al fluido, e h = ρ(1 + ε) + p e l’entalpiarelativistica specifica 1. Accanto alle equazioni di Einstein nello spazio-tempo dovranno quindi esserci anche quelle dell’idrodinamica relativisticaper l’evoluzione della materia, ovvero dovranno essere implementate lalegge di conservazione del tensore energia impulso ∇µTµν = 0 , la leggedi conservazione del numero barionico ∇µ(ρuµ) = 0 e l’equazione distato (EOS) del tipo p = p(ρ, ε). Per i nostri propositi bastava restringerel’attenzione su di un’equazione di stato politropica del tipo

p = KρΓ, ε =K

Γ− 1ρΓ−1, (5.3)

con K = 100 e Γ = 2.

5.3 CACTUS-CARPET-WHISKY:codice non-lineare 3Dper la relativita generale

Abbiamo usato la formulazine “3 + 1” conforme senza traccia delleequazioni di Einstein, nella quale lo spazio-tempo e decomposto in ipersu-perfici tri-dimesionali, descritte da una metrica γij , immerse nello spaziototale a quattro dimensioni, specificate dalla curvatura estrinseca Kij edalle funzioni di gauge α (lapse) e βi (shift) che determinano un sistema

1In termodinamica l’entalpia H e definita come la somma tra energia interna U e ilprodotto pressione per volume, H = U + pV . In relativita si aggiunge anche l’energia dellamassa a riposo M , cosi che l’entalpia diventa H =M +U +pV =M(1+ ε)+pV . L’entalpiaspecifica e quindi definita come entalpia per unita di massa: h = H/M = 1 + ε+ p

ρ.

54

5.3. CACTUS-CARPET-WHISKY:codice non-lineare 3D per la relativitagenerale

di riferimento. Trasformiamo le variabili standard ADM come segue. Latre-metrica γij e trasformata in modo conforme

γij = e−4φγij , φ =1

12ln detγij (5.4)

e il fattore conforme φ e evoluto come variabile indipendente, mentre γije soggetta al constraint detγij = 1. La curvatura estrinseca e soggetta allostessa trasformazione conforme e la sua traccia e evoluta come variabileindipendente. Invece di Kij abbiamo evoluto

K = trKij = γijKij , Aij = e−4φ(Kij −

1

3γijK

), (5.5)

con trAij = 0. Infine le nuove veriabili da evolvere definite in termini deisimboli di Christoffel e della tre-metrica

Γi = γjkΓijk. (5.6)

Le equazioni di Einstein specificano un set di equazioni ben noto, dato da

d

dtγij = −2αAij , (5.7)

d

dtφ = −1

6αK, (5.8)

d

dtAij = e−4φ

−DiDjα+ αRij + 4πα[γij(S − ρ)− 2Sij ]

TF+ α(KAij − 2AikAkj),

(5.9)

d

dtK = −DiD

iα+ α(AijA

ij +1

3K2)

+ 4πα(ρ+ S), (5.10)

con d/dt := ∂t − £~βe dove TF denota la parte libera dalla traccia del

termine tra parentesi. Nell’evoluzione di evoluzione di K e stato usatol’energy constraint per eliminare lo scalare di Ricci. Abbiamo anche

∂tΓi = γjk∂j∂kβ

i +1

3γij∂j∂kβ

k − 2Aij∂jα

+ 2α(

ΓijkAjk + 6Aij∂jφ−

2

3γij∂jK − 8πji

).

(5.11)

ρ, Sij , ji sono termini sorgenti per la materia definiti come

ρ ≡ nαnβTαβ, (5.12)

ji ≡ −γiαnβTαβ, (5.13)

Sij ≡ γiαγjβTαβ (5.14)

con nα = (−α, 0, 0, 0) vettore normale alle ipersuperfici. Le equazioni diEinstein portano anche a un set di constraints dentro ogni ipersuperfice

H ≡ R(3) +K2 −KijKij − 16πρ = 0 (5.15)

Mi ≡ Dj(Kij − γijK)− 8πji = 0, (5.16)

55

5.4. Gauges

Qui R(3) = Rijγij . La particolare scelta fatta introduce cinque ulteriori

constraints

detγij = 1, (5.17)

trAij = 0, (5.18)

Γi = γjkΓijk. (5.19)

Il codice utilizzato implemente i constraints algebrici (5.17) e (5.18). I con-straints rimanenti, H,Mi e (5.19) non sono stati implementati e possonoessere utilizzati per monitorare la bonta della simulazione.

5.4 Gauges

Abbiamo specificato il gauge in termini delle variabili standard diArnowitt-Deser-Misner (ADM), α funzione di lapse e βi vettore di shift.Abbiamo evoluto il lapse in accordo al “1 + log” slicing condition:

∂t − βi∂iα = −2αK. (5.20)

Il vettore di shift e stato evoluto tramite la condizione Γ-driver iperbolica:

∂tβi − βj∂jβi =

3

4αBi, (5.21)

∂tBi − βj∂jBi = ∂tΓ

i − βj∂jΓi − ηBi, (5.22)

dove η e un parametro che agisce come coefficente di dumping.Tutte le equazioni scritte sopra sono state risolte utilizzando un codice

tri-dimensionale alle differenze finite contenuto nel Cactus ComputationalToolkit.

5.5 Evoluzione della materia

Le equazioni per l’idrodinamica in relativita generale sono state ri-solte con il codice WHISKY. Una caratteristica importante del codice el’implementazione di una flux conservativ formulation delle equazioni del-l’idrodinamica, in cui il set di equazioni di conservazione per il tensoreenergia-impulso Tµν e la densita di corrente di materia Jµ = ρuµ, ovvero

∇µTµν = 0, ∇µJµ = 0, (5.23)

e scritto in forma iperbolica, al primo ordine e tipo conservazione di flusso,ovvero

∂tq + ∂ifi(q) = s(q), (5.24)

dove fi(q) e s(q) sono il vettore di flusso e la sorgente rispettivamente.Per riscrivere le equazioni di conservazione in questo modo, [21], e

necessario mappare la ρ densita nel sistema a riposo, la pressione delfluido p misurata nel sistema a risposo, la tri-velocita vi misurata da unosservatore locale con momento angolare nullo, l’energia specifica ε, e il

56

5.6. Estrazione delle onde con lo scalare di Newmann-Penrose Ψ4

fattore di Lorentz W , nelle cosi chiamate variabili conservate q ≡ (D,Si, τ)

tramite le seguenti relazioni

D ≡ √γWρ, (5.25)

Si ≡ √γρHW 2vi, (5.26)

τ ≡ √γ(ρHW 2 − p)−D. (5.27)

E’ necessario osservare che nel caso di un’equazione generale EOS del tipop = P (ρ, ε) solo cinque delle sette variabili primitive sono indipendenti.Se adottiamo pero un equazione del tipo p = p(ρ), allora pressione edenergia specifica sono completamente determinate dalla densita a riposo el’equazione (5.27) e ridondante e puo evitare di essere risolta.

5.6 Estrazione delle onde con lo scalare di Newmann-Penrose Ψ4

L’uso degli scalari di Weyl per l’estrazione delle onde e diventato di usocomune in relativita numerica negli ultimi anni.

Data una ipersuperfice con un vettore di tipo tempo normale nµ e datoun vettore unitario rµ che indica la direzione di propagazione delle onde,ricordo che la definizione standard dello Ψ4, dato il tensore di curvatura diWeyl, e

Ψ4 = −Cαµβνkµkνmαmβ (5.28)

dove kµ ≡ 1/√

2(nµ − rµ) e mµ e un vettore complesso nullo tale che(mµmµ = 1) che e ortogonale sia a nµ sia a rµ. Abbiamo definito una baseortonormale nello spazio a tre dimensioni (er, eθ, eφ), centrata sull’originedegli assi cartesiani e orientata con i poli lungo l’asse z. Il vettore normalealle ipersuperfici definisce il vettore tipo tempo, possiamo dunque costruirela tedrade nulla come

k =1√2

(et − er), (5.29)

l =1√2

(et + er), (5.30)

m =1√2

(eθ + ieφ), (5.31)

(5.32)

Quando calcoliamo Ψ4 tramite (5.28) in termini delle variabili ADM [22]su di uno slice

ψ4 = Cijmimj , (5.33)

doveCij ≡ Rij −KKij +Kk

i Kkj − iεkli ∇lKjk. (5.34)

L’ampiazza della polarizzazione delle onde gravitazionali h+ e h× sonolegate a ψ4 da

h+ − ih× = Ψ4. (5.35)

57

5.7. Estrazione delle onde tramite multipoli

Se usiamo l’espansione di Ψ4 in armoniche sferiche di spin pesate con pesos = −2

Ψ4(t, r, θ, φ) =∞∑l=2

l∑m=−l

Ψl,m4 (t, r)−2Y

l,m(θ, φ), (5.36)

in modo che Ψl,m4 e i multipoli della metrica hl,m sono legati da

hl,m(t, r) = Ψl,m4 (t, r), (5.37)

hl,m(t, r) e allora il doppio integrale indefinito di Ψl,m4 (t, r), che dopo aver

moltiplicato per r entrambi i membri di (5.37), si puo calcolare computazio-nalmente come

rhl,m(t, r) ≡∫ t

0dt′∫ t′

0dt′′

(5.38)

che risultarhl,m(t, r) = rhl,m(t, r) +Q0 + tQ1, (5.39)

dove sono state scritte esplicitamente le due costanti di integrazione Q0 eQ1. Possono essere determinate dai dati e il loro significato fisico e Q0 =

rhl,m(0, r) e Q1 = rhl,m(0, r).E’ necessario osservare tuttavia, che l’equazione discussa si riferisce ad

un segnale estratto ad una distanza finita, mentre cio che interessa calcolaree Ψ4 ad una distanza spaziale infinita. Si puo immaginare che ci sia unoffset da tenere presente ovvero dovremmo poter scrivere

rΨl,m4 (t) ≡ rΨl,m

4 (t, r) + 2Q2(r) (5.40)

con Q2(r) funzione di offset. L’integrazione temporale dovrebbe quindiessere

rhl,m(t, r) = rhl,m(t, r) +Q0 + tQ1 + t2Q2(r). (5.41)

5.7 Estrazione delle onde tramite multipoli

L’estrazione delle onde basata sulla teoria delle perturbazioni di unospaziotempo di Schwarzschild fu introdotta da Abrahams e Price [23].L’assunzione che viene fatta e che lontano da regioni di campo forte, lospazio-tempo numerico e ben approssimato dalla somma di una metricasfericamente simmetrica di Schwarzschild g(0)

µν e di una perturbazione hµνnon sferica. Come visto nel primo capitolo e possibile scrivere

gµν(r, t, θ, φ) = g(0)µν +

∞∑l=2

l∑m=−l

[(hl,mµν )(o) + (hl,mµν )e

]. (5.42)

I multipoli (hl,mµν )(o/e) e le loro derivate possono essere combinati insiemein due master functions gauge invarianti, a parita pari (Zerilli-Moncrief)Ψ

(e)l,m e a parita dispari (Regge-Wheeler) Ψ

(o)l,m. Queste due master functions

soddisfano due disaccoppiate equazioni delle onde con potenziale.

58

5.8. I dati iniziali

Tabella 5.1: Proprieta di equilibrio del modello A0. Da sinistra a destra le colonneriportano: densita centrale a riposo, densita totale di energia centrale,massa gravitazionale, raggio, compattezza.

ρc ec M R M/R

A0 1.28× 10−3 1.44× 10−3 1.40 9.57 0.15

5.8 I dati iniziali

Come modello rappresentativo di stella a neutroni, abbiamo scelto unmodello descritto da una EOS politropica con Γ = 2 e K = 100. Abbiamoconsiderato una densita centrale ρc = 1.28×10−3 cosi che la massa a riposovale approssimativamente M ' 1.4. In letteratura questo modello e statoampiamente utilizzato ed e conosciuto come modello A0 in [24]. Nellatabella (5.1) sono riepilogate alcune proprieta del modello.

E’inoltre necessario

1. risolvere le equazioni di Tolman-Oppenheimer-Volkov (TOV) pergenerare la configurazione di equilibrio

2. fissare una pressione perturbante assisimmetrica

3. risolvere i constraints linearizzati per la perturbazione della metrica.

La pressione perturbante e inizializzata come un multipolo assissimme-trico:

δp(r, θ) ≡ (p+ e)Hl0(r)Yl0(θ). (5.43)

Abbiamo limitato lo studio al caso l = 2 e m = 0. Dato il nostro interesseriguardava il confronto tra la forma delle onde e non la fisica dei processidi oscillazione della stella di neutroni, abbiamo considerato un sistemaoscillante esattamente ad una frequenza, ovvero Hl,0 corrisponde ad unaautofunzione della stella. Un autofunzione approssimata, come suggeritoda [25], per un dato modo di un fluido puo essere espressa da

Hl,0 = λsin[(n+ 1)πr

2R

], (5.44)

dove n e un intero che controlla il numero di nodi di Hl,0, λ e l’ampiezzadella perturbazione e R e il raggio della stella nelle coordinate di Sch-warzschild. Per semplicita ci siamo occupati solamente del caso n = 0 eλ = [0.1, 0.05, 0.01, 0.001].

59

5.9. Simulazioni

5.9 Simulazioni

Sono state inizialmente eseguite delle simulazioni di prova senza uti-lizzare la perturbazione, tabella (5.2). Questo ha consentito di prenderedimestichezza con il programma.

Tabella 5.2: Test preliminari.

dx rl A

Test 1 8 5 0Test 2 4 5 0

Come risulta evidente dalle figure in 5.1, il passo reticolare gioca unruolo fondamentale per la corretta riuscita di un test. Si puo notare comel’instabilita, dovuta allo scarso raffinamento reticolare, si manifesta nelgrafico che riporta la variazione della densita centrale per i due test equello che riporta la variazione temporale del constraint hamiltoniano. E’interessante notare anche il drift con oscillazioni per il test con dx = 4.Come descritto nei paragrafi precedenti ci si aspettava un oscillazione nelladensita, che come si puo vedere aumenta all’avvicinarsi dell’instabilitanumerica.

(a) Confronto densita massime per testtabella (5.2).

(b) Oscillazioni con dx = 4.

(c) Constraints hamiltoniano .

Figura 5.1: Alcuni test iniziali riferiti alla tabella 5.2

60

5.9. Simulazioni

Sono state eseguite altre otto simulazioni variando l’ampiezza dellaperturbazione e il raffinamento della griglia, tabelle 5.3 e 5.4. Su questi testsono stati messi alla prova i due differenti metodi di estrazione del segnaledi onde gravitazionali descritti.

Tabella 5.3: Test eseguiti con la massima risoluzione, l’intervallo sulla griglia piupiccolo e 0.250

dx rl A

Test 3 2 4 0.1Test 4 2 4 0.05Test 5 2 4 0.01Test 6 2 4 0.001

Tabella 5.4: Test eseguiti con risoluzione inferiore, l’intervallo sulla griglia piupiccolo e 0.375

dx rl A

Test 7 3 4 0.1Test 8 3 4 0.05Test 9 3 4 0.01Test 10 3 4 0.001

Per prima cosa e stato necessario inserire e attivare i thorns seguenti:Spherical Harmonics e Whisky_PerturbTOV. Grazie a questi pro-grammi possiamo generare la perturbazione della configurazione iniziale ecalcolare la Ψ4 con Psikadelia.

Con questi parametri modifichiamo l’ampiezza della perturbazione edimpostiamo i raggi di estrazione della Ψ4.

61

5.9. Simulazioni

Listing 5.1: File dei parametri, Spherical Harmonics1 #=========================================================

#3 # Spherical Harmonics

#5 #=========================================================

7

SphericalHarmonics::number of radii = 189 SphericalHarmonics::ex radii [0] = 30

SphericalHarmonics::ex radii [1] = 3511 SphericalHarmonics::ex radii [2] = 40






SphericalHarmonics::ex radii [13] = 9523 SphericalHarmonics::ex radii [14] =100

SphericalHarmonics::ex radii [15] =10525 SphericalHarmonics::ex radii [16] =110

SphericalHarmonics::ex radii [17] =11527

SphericalHarmonics::lmax = 629 SphericalHarmonics::interp integration order = 4

SphericalHarmonics::grid type = cart3d31 SphericalHarmonics::InterpPointsTheta = 49

SphericalHarmonics::InterpPointsPhi = 15633 SphericalHarmonics::number of vars = 2

SphericalHarmonics::vars [0] = ”PsiKadelia::psi4re”35 SphericalHarmonics::SH spin weight [0] = −2

SphericalHarmonics::vars [1] = ”PsiKadelia::psi4im”37 SphericalHarmonics::SH spin weight [1] = −2

62

5.9. Simulazioni

Listing 5.2: File dei parametri, PerturbTOV1 #=========================================================

#3 # Whisky PerturbTOV

#5 #=========================================================

7 Whisky PerturbTOV::PertTOV Num Radial = 3000Whisky PerturbTOV::TOV Npoints Theta = 50

9 Whisky PerturbTOV::TOV Npoints Phi = 50

11 #−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−# PERTURBATION TYPE ELL=2

13 #−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−

15 Whisky PerturbTOV::l mode = 2Whisky PerturbTOV::constraints = ”yes”

17 Whisky PerturbTOV::pert type = ”matter”Whisky PerturbTOV::matter pert type = ”enthalpy”

19

#−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−21 # PERTURBATION AMPLITUDE

# LAMBDA = [0.001, 0.01, 0.05, 0.1]23 # −> 0 1 2 3

#−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−25

Whisky PerturbTOV::pert amplitude = 0.0127

#−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−29 # PERTURBATION NODES

#−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−31

Whisky PerturbTOV::Nnodes = 0.5 # 0 nodes33 Whisky PerturbTOV::aleph = 0.0

Listing 5.3: File dei parametri, PerturbTOV#=========================================================

2 ## IO 2D, Psi4 scomposto

4 ##=========================================================

6

# for waves via Psikadelia8

IOASCII::out2D dir = ”data/2d”10 IOASCII::out2D criterion = divisor

IOASCII::out2D every = 412 IOASCII::out2D vars = ”SphericalHarmonics::decomposed vars”

63

5.9. Simulazioni

Nel seguito riporto i confronti tra le densita ottenute con i vari modelli.

Figura 5.2: Variazione di ρmax nel tempo per la tabella (5.3).

Figura 5.3: Variazione di ρmax nel tempo per la tabella (5.4).

Ovviamente l’ampiezza delle oscillazioni della densita massima au-menta all’aumentare dell’ampiezza della perturbazione. Si evidenzia dasubito pero un problema, come si puo osservare la configurazione iniziale

64

5.9. Simulazioni

Figura 5.4: Profilo densita iniziale e finale per il Test 3. In verde la configurazione iniziale.

Figura 5.5: Confronto densita iniziali per i test della tabella (5.3).

65

5.9. Simulazioni

presenta un “dente” che con il proseguire dell’evoluzione viene limato.Questo sta a significare che la perturbazione iniziale non e stata implemen-tata correttamente e quindi non si va piu a perturbare il sistema con unautofunzione di modo l = 2, m = 0, ma con qualche cosa di lievementedifferente. Cio si e ripercosso, come mostrero nel seguito, sull’estrazionedel segnale delle onde cercato.

Le figure 5.6 visualizzano il comportamento dei constraints. In partico-lare si osserva una dissipazione del constraint del momento lungo l’asse x euna crescita di quello hamiltoniano, tanto piu rapida quanto piu piccola e laperturbazione utilizzata. Le simulazioni svolte con il raffinamento peggioremostrano una maggiore instabilita, la crescita del constraint hamiltonianoavviene piu precocemente..

(a) Verifica constraints hamiltoniani per latabella (5.3).

(b) Verifica constraints hamiltoniani per latabella (5.4).

(c) Verifica momento constraints per latabella (5.3), lungo x.

(d) Verifica momento constraints per latabella (5.4), lungo x.

Figura 5.6: Controllo dei constraints del momento e hamiltoniani per le tabelle(5.3) e (5.4) .

Un altro elemento di controllo e visualizzare p/ρ2, che permette dicapire se l’EOS utilizzata viene mantenuta per tutta la simulazione. Effetti-vamente cio che il grafico 5.7 mette in evidenza e che nei punti dentro lastella K = 100, ovvero p segue l’equazione di stato politropica con Γ = 2

e K = 100. Le creste che si osservano sono situate al bordo della stella,nell’atmosfera che si viene a creare.

Visionati i grafici di controllo possiamo passare all’analisi dell’estrazio-ne delle onde vera e propria.

66

5.9. Simulazioni

Figura 5.7: Controllo dell’ EOS per l’ultimo step di dell’evoluzione.

I grafici 5.8, 5.9, 5.10, 5.11, mettono in evidenza le funzioni di Zerillitrovate partendo dai dati di WaveExtract.

−100 0 100 200 300 400 500 600 700 800 900

−0.15

−0.1

−0.05

0

0.05

0.1

u=t−robs

*

Ψ(e)

20

A= 0.1; radius4 (115)

A= 0.05; radius4 (115)

A= 0.01; radius4 (115)

A= 0.001; radius4 (115)

Figura 5.8: Confronto Ψe ottenute dalla tabella (5.3) con il metodo di Abrahams e Price.

Sono state messe a confronto le Ψe per le diverse ampiezze della pertur-bazione a raggio fissato e per diversi raggi ad ampiezza fissata. I confronti araggi differenti sono stati plottati sfruttando la variabile u = t− r∗ tempo ri-tardato, dove r∗ e la coordinata della tartaruga di Regge-Wheeler. Tutto ciocon i due differenti raffinamenti. Si evidenzia, per le perturbazioni maggiori

67

5.9. Simulazioni

−200 −100 0 100 200 300 400 500 600−0.2

−0.15

−0.1

−0.05

0

0.05

0.1

0.15

u=t−robs

*

Ψ(e)

20

A= 0.1; radius4 (115)

A= 0.05; radius4 (115)

A= 0.01; radius4 (115)

A= 0.001; radius4 (115)

Figura 5.9: Confronto Ψe ottenute dalla tabella (5.4) con il metodo di Abrahams e Price.

−200 0 200 400 600 800 1000 1200−0.2

−0.15

−0.1

−0.05

0

0.05

0.1

0.15

u=t−robs

*

Ψ(e)

20

A= 0.1; radius1 (30)

A= 0.1; radius2 (60)

A= 0.1; radius3 (90)

A= 0.1; radius4 (115)

Figura 5.10: Confronti Ψe ottenute per la tabella (5.3) al variare del raggio di estrazione.

68

5.9. Simulazioni

−200 −100 0 100 200 300 400 500 600 700−0.2

−0.15

−0.1

−0.05

0

0.05

0.1

0.15

u=t−robs

*

Ψ(e)

20

A= 0.1; radius1 (30)

A= 0.1; radius2 (60)

A= 0.1; radius3 (90)

A= 0.1; radius4 (115)

Figura 5.11: Confronti Ψe ottenute per la tabella (5.4) al variare del raggio di estrazione.

uno sfasamento a circa meta della simulazione. Tale shift nelle frequenze estato causato da una variazione della massa della stella, a sua volta causatodal fatto che abbiamo cercato di far evolvere il sistema con una determinataTOV purtroppo resa inefficace, avendo la cattiva perturbazione inizialemodificato le condizioni da noi poste.

Le funzioni di Zerilli presentano un andamento non fisico nella parteiniziale dell’evoluzione. Questo lo si evince dal fatto che l’ampiezza diquesto junk aumenta all’aumentare del raggio di estrazione come mostrala figura 5.12, invece di decrescere e raggiungere un valore costante.

L’immagine 5.13 confronta le H := h+− ih× ottenute dalla Ψ4. In realtaqui non e ancora stato tolto il fondo con andamento quadratico discusso neiparagrafi precedenti. Tramite un fit delle funzioni sopra con un polinomiodi grado due e stato possibile eliminare il “fondo” dovuto all’integrazionedoppia (vedere precedenti paragrafi), ottenendo la funzione H effettiva-mente cercata. E’ stato quindi possibile confrontare il segnale ricercato conil metodo della Ψ4 e con il metodo di Abrahams e Price, usando per questiultimi i dati procurati con WaveExtract. Le figure 5.14, 5.15, 5.16, 5.17, met-tono in risalto tale comparazione. Appare subito che il segnale dell’ondaottenuto con il metodo degli scalari di Weyl risulta essere maggiormen-te pulito, rispetto a quello ottenuto dalle funzioni di Zerilli-Moncrief. Ilrumore sovrapposto al segnale di queste ultime e dovuto anch’esso dallepessime condizioni iniziali. Tuttavia come mostra il grafico ottenuto conla perturbazione di ampiezza minore, anche il segnale ottenuto con la Ψ4

evidenzia una modulazione a frequenza inferiore. Si evidenzia anche unjunk iniziale che si ottiene solo con il metodo di Abrahams e Price. Tale ef-

69

5.9. Simulazioni

20 40 60 80 100 120 140 160 180 200

−0.1

−0.05

0

0.05

0.1

0.15

t

Ψ(e)

20

A= 0.1; radius1 (30)

A= 0.1; radius2 (60)

A= 0.1; radius3 (90)

A= 0.1; radius4 (115)

Figura 5.12: Andamento lineare dell’ampiezza del junk al variare del raggio diestrazione.

0 200 400 600 800 1000 1200−16

−14

−12

−10

−8

−6

−4

−2

0

2x 10

−3

t

Ψ42

,0

Ψ4 test

A = 0.1

A = 0.05

A = 0.01

A = 0.001

Figura 5.13: Controllo delle H := h+ − ih× ottenute con il metodo degli scalaridi Weyl.

70

5.9. Simulazioni

Figura 5.14: Ampiezza perturbazione λ = 0.1.


71

5.9. Simulazioni



72

5.9. Simulazioni

fetto diminuisce al diminuire dell’ampiezza della perturbazione, potrebbequindi anch’esso essere considerato una manifestazione dell’inaccuratezzadel set-up iniziale. La differenza delle ampiezze ottenute con i due diffe-renti metodi, per quanto sia difficile valutare le funzioni di Zerilli ricche dirumore, si aggira intorno al 10%.

Infine vengono mostrate le trasformate di fourier per i segnali sopra, inmodo da rendere piu marcate e pulite, le differenze riscontrate.


La frequenza fondamentale per il metodo di estrazione degli scalari diWeyl risulta essere sempre maggiore, circa 9−10% rispetto a quella ottenutacon lo sviluppo in multipoli. La tabella 5.5 riporta i valori delle frequenzefondamentali per i diversi modelli di perturbazione ed estrazione. Sonoriportate esclusivamente quelle ottenute dalla simulazione con il migliorraffinamento.

Tabella 5.5: Frequenze dominanti delle H per i vari modelli perturbati.

A WaveExtract [Hz] Ψ4 [Hz]

0.1 1624 16260.05 1574 15740.01 1576 15800.001 1579 1594

73

5.9. Simulazioni



74

5.9. Simulazioni


75

CAPITOLO 6

CONCLUSIONI

Forse il nostro universo si trova dentro al dente diqualche gigante.

ANTON CECHOV

Sono stati comparati due metodi di estrazione del segnale di onde gra-vitazionali ad oggi popolari nelle simulazioni numeriche di relativita gene-rale: (i) il metodo di Abrahams-Price basato sulla teoria delle perturbazionigauge-invarianti di Regge-Wheeler-Zerilli-Moncrief dello spazio-tempodi Schwarzschild; (ii) metodo di estrazione basato sugli scalari di Weyl,nota la funzione Ψ4. E’ stata utilizzata come sorgente delle onde una stelladi neutroni oscillante non rotante con equazione di stato politropica. E’stata aggiunta una perturbazione alla pressione iniziale. Le simulazionisono state eseguite sfruttando il codice CACTUS-CARPET-WHISKY, chepermette la risoluzione delle equazioni non lineari di Einstein in relativitagenerale in tre dimensioni. Il formalismo adottato e quello detto “3+1” chepermette di riscrivere le equazioni di Einstein nel modo chiamato BSSNOK.

La perturbazione della pressione doveva essere un autofunzione ap-prossimata della stella. La nostra attenzione e caduta solamente sul modol = 2 e m = 0.

L’accuratezza della frequenza fondamentale per i vari modelli di per-turbazione e raffinamento e stata valutata prendendo a confronto i valoririportati in [1]. I valori ottenuti sono in perfetto accordo con quelli riportatinell’articolo citato. Tuttavia si sono riscontrati problemi dovuti all’inaccu-ratezza della perturbazione iniziale che ha generato molto rumore sullefunzioni di Zerilli oltre ad averne addirittura variato la frequenza (neimodelli con perturbazione piu grande λ = 0.1 e λ = 0.05), facendo variarela massa della stella da 1.44 a 1.52.

La funzioni di Zerilli presentano un andamento non fisico nella parteiniziale dell’evoluzione. Questo lo si evince dal fatto che l’ampiezza di

77

questo junk aumenta all’aumentare del raggio di estrazione come mostrala figura 5.12, invece di decrescere e raggiungere un valore costante.

Il segnale estratto tramite la Ψ4 sembra piu pulito, non presenta il junkiniziale, ma quando l’ampiezza della perturbazione e piccola si evidenziauna modulazione che ne distorce l’andamento.

78

APPENDICE A

SPIN-WEIGHTED SPHERICAL HARMONICS

I fisici vengono addestrati a studiare un problema primadi arrivare a una decisione. Gli avvocati, i pubblicitari ealtri sono addestrati a fare esattamente il contrario:cercare dati per confermare una determinazione che egia stata assunta.

ROBERT CREASE

Quando lavoriamo con funzioni non scalari definite sulla sfera intro-duciamo le cosi chiamate spin-weighted spherical harmonics, una generaliz-zazione delle ordinarie armoniche sferiche. Le spin-weighted sphericalharmonics furono introdotte per la prima volta da Newman e Penrose [26]per lo studio di radiazione gravitazionale, ma possono essere usate ancheper lo studio delle equazioni di Maxwell, per l’equazione di Dirac e ingenerale per equazioni con campi di spin arbitrario.

Consideriamo una funzione complessa f sulla sfera che potrebbe cor-rispondere ad una qualche combinazione di componenti di un oggettotensoriale (o spinoriale) nella base ortonormale (r, θ, ϕ) diremo che f haspin-weight s se, sotto rotazione della base angolare (eθ, eϕ) di un angolo Ψ,essa trasforma come f → e−isΨf .

Le spin-weighted spherical harmonics, denotate con sYl,m(θ, ϕ), formano

una base per lo spazio delle funzioni con definito spin-weight s. Possonoessere introdotte in vari modi, per esempio partendo dalla definizione deglioperatori

ðf := −sinsθ(∂θ +

i

sinθ∂ϕ

)(fsin−sθ)

= −(∂θ +

i

sinθ∂ϕ − s cotθ

)f,

(A.1)

ðf := −sinsθ(∂θ −

i

sinθ∂ϕ

)(fsin−sθ)

= −(∂θ −

i

sinθ∂ϕ + s cotθ

)f,

(A.2)

79

dove s e lo spin-weight di f . Le spin-weighted spherical harmonics allora sonodefinite per |m|≤ l e l ≥ |s| in termini delle armoniche sferiche standard

sY l,m :=[(l − s)!

(l + s)!

]1/2ðs(Y l,m), +l ≥ s ≥ 0, (A.3)

sY l,m := (−1)s[(l + s)!

(l − s)!

]1/2ð−s(Y l,m), −l ≤ s ≤ 0. (A.4)

In particolare si ha 0Yl,m = Y l,m. La definizione sopra implica che

ð(sYl,m) = +[(l − s)(l + s+ 1)]1/2 s+1Y

l,m, (A.5)

ð(sYl,m) = −[(l + s)(l − s+ 1)]1/2 s−1Y

l,m. (A.6)

Per questo motivo gli operatori ð e ð vengono detti operatori di abbassa-mento e innalzamento. Si trova anche che

ðð(sYl,m) = +[l(l + 1)− s(s+ 1)]1/2 sY

l,m, (A.7)

ðð(sYl,m) = +[l(l + 1)− s(s− 1)]1/2 sY

l,m, (A.8)

(A.9)

cosi che le sYl,m sono autofunzioni di ðð e ðð con gli autovalori di L2.

Dalle definizioni sopra e possibile mostrare che

sYl,m(θ, ϕ) = (−1)s+m −sY

l,−m(θ, ϕ). (A.10)

80

APPENDICE B

CODICE

Ci sono solamente 10 tipi di persone nel mondo: chicomprende il sistema binario e chi no.

ANONIMO

Nel seguito riporto la configurazione usata per compilare il codice perle nostre simulazioni, la nostra Thornlist e un esempio di file di parametri.

B.1 Configurazione

Configurazione usata per la compilazione dell’Einstein Toolkit (fileconfig):

# CONFIGURATION : perturbTOVPhtread# CONFIG-DATE : Mon Mar 21 10:40:22 2011 (GMT)# CONFIG-HOST : einstein.pr.infn.it# CONFIG-STATUS : 0# CONFIG-OPTIONS :BLAS=yesBLAS_DIR=/opt/External/gcc4.1.2/libBLAS_LIBS=blasCC=gccCXX=g++F77=gfortranF90=gfortranGSL=yesGSL_DIR=/opt/External/gcc4.1.2HDF5=yesHDF5_DIR=/opt/External/gcc4.1.2LAPACK=yesLAPACK_DIR=/opt/External/gcc4.1.2/libLAPACK_LIBS=lapackLIBSZ_DIR=/opt/External/gcc4.1.2/lib/LIBZ_DIR=/opt/External/gcc4.1.2/lib/MPI=OpenMPIOPENMP=noOPENMPI_DIR=/usr/lib64/openmpi/1.4-gcc/PTHREAD=yes

81

B.2. Thornlist

B.2 Thornlist

ThornList utilizzata per la compilazione (file ThornList):

CactusArchive/ADMCactusBase/BoundaryCactusBase/CartGrid3DCactusBase/CoordBaseCactusBase/FortranCactusBase/IOASCIICactusBase/IOBasicCactusBase/IOUtilCactusBase/InitBaseCactusBase/LocalInterpCactusBase/LocalReduceCactusBase/SymBaseCactusBase/Time

CactusConnect/HTTPDCactusConnect/HTTPDExtraCactusConnect/Socket

CactusElliptic/EllBaseCactusElliptic/EllSOR

CactusNumerical/MoLCactusNumerical/Cartoon2DCactusNumerical/NoiseCactusNumerical/NormsCactusNumerical/PeriodicCactusNumerical/ReflectionSymmetryCactusNumerical/RotatingSymmetry180CactusNumerical/RotatingSymmetry90CactusNumerical/SlabCactusNumerical/SlabTestCactusNumerical/SpaceMaskCactusNumerical/SphericalSurface

CactusPUGH/PUGHCactusPUGH/PUGHInterpCactusPUGH/PUGHReduceCactusPUGH/PUGHSlabCactusPUGHIO/IOHDF5UtilCactusPUGHIO/IOHDF5

CactusUtils/FormalineCactusUtils/NaNCheckerCactusUtils/NiceCactusUtils/TerminationTriggerCactusUtils/TimerReport

# The Einstein Toolkit!TARGET = $ARR!TYPE = svn!URL = https://svn.einsteintoolkit.org ...

/cactus/$1/$2/branches/$ET_RELEASE!CHECKOUT =EinsteinAnalysis/ADMAnalysisEinsteinAnalysis/ADMConstraintsEinsteinAnalysis/AHFinderEinsteinAnalysis/AHFinderDirectEinsteinAnalysis/CalcKEinsteinAnalysis/EHFinderEinsteinAnalysis/Extract

82

B.2. Thornlist

EinsteinAnalysis/Hydro_AnalysisEinsteinAnalysis/MultipoleEinsteinAnalysis/WeylScal4EinsteinBase/ADMBaseEinsteinBase/ADMCouplingEinsteinBase/ADMMacrosEinsteinBase/ConstantsEinsteinBase/CoordGaugeEinsteinBase/EOS_BaseEinsteinBase/EOSG_BaseEinsteinBase/HydroBaseEinsteinBase/StaticConformalEinsteinBase/TmunuBaseEinsteinEOS/EOSG_HybridEinsteinEOS/EOSG_IdealFluidEinsteinEOS/EOSG_PolytropeEinsteinEOS/EOS_HybridEinsteinEOS/EOS_IdealFluidEinsteinEOS/EOS_PolytropeEinsteinEvolve/LegoExcisionEinsteinEvolve/NewRadEinsteinInitialData/DistortedBHIVPEinsteinInitialData/ExactEinsteinInitialData/IDAnalyticBHEinsteinInitialData/IDAxiBrillBHEinsteinInitialData/IDAxiOddBrillBHEinsteinInitialData/IDBrillDataEinsteinInitialData/IDConstraintViolateEinsteinInitialData/IDFileADMEinsteinInitialData/IDLinearWavesEinsteinInitialData/Meudon_Bin_BHEinsteinInitialData/Meudon_Bin_NSEinsteinInitialData/Meudon_Mag_NSEinsteinInitialData/NoExcisionEinsteinInitialData/RotatingDBHIVPEinsteinInitialData/TwoPuncturesEinsteinUtils/SetMask_SphericalSurfaceEinsteinUtils/TGRtensor

# Additional Cactus thorns!TARGET = $ARR!TYPE = svn!AUTH_URL = https://svn.cactuscode.org ...

/projects/$1/$2/branches/$ET_RELEASE!URL = http://svn.cactuscode.org ...

/projects/$1/$2/branches/$ET_RELEASE!CHECKOUT =ExternalLibraries/BLASExternalLibraries/GSLExternalLibraries/HDF5ExternalLibraries/LAPACKExternalLibraries/LORENEExternalLibraries/zlib

TAT/TATelliptic

# Various thorns from the AEI!TARGET = $ARR!TYPE = svn!URL = http://svn.aei.mpg.de/numrel/ ...

$1/$2/branches/$ET_RELEASE!AUTH_URL = https://svn.aei.mpg.de/numrel/ ...

$1/$2/branches/$ET_RELEASE!CHECKOUT =

83

B.2. Thornlist

AEIThorns/AEILocalInterpAEIThorns/BSSN_MoLAEIThorns/ManualTerminationAEIDevelopment/PsiKadeliaAEIDevelopment/WaveExtract

# Various thorns from LSU!TARGET = $ARR!TYPE = svn!URL = https://svn.cct.lsu.edu/repos/ ...

numrel/$1/$2/branches/$ET_RELEASE!CHECKOUT =LSUThorns/QuasiLocalMeasuresLSUThorns/SummationByPartsLSUThorns/SphericalHarmonics

# From Kranc (required e.g. by McLachlan)!TARGET = $ARR!TYPE = git!URL = git://github.com/ianhinder/Kranc.git!AUTH_URL = [email protected]:ianhinder/Kranc.git!REPO_PATH= Auxiliary/Cactus!REPO_BRANCH = ET_2010_06!CHECKOUT =KrancNumericalTools/GenericFD

# McLachlan, the spacetime code!TARGET = $ARR!TYPE = git!URL = git://carpetcode.dyndns.org/McLachlan!AUTH_URL = [email protected]:McLachlan!REPO_PATH= $2!REPO_BRANCH = ET_2010_06!CHECKOUT = McLachlan/doc McLachlan/m McLachlan/parMcLachlan/ML_BSSNMcLachlan/ML_BSSN_HelperMcLachlan/ML_BSSN_O2McLachlan/ML_BSSN_O2_HelperMcLachlan/ML_BSSN_TestMcLachlan/ML_ADMConstraintsMcLachlan/ML_ADMQuantities

# Carpet, the AMR driver!TARGET = $ARR!TYPE = git!URL = git://carpetcode.dyndns.org/carpet!AUTH_URL = [email protected]:carpet!REPO_BRANCH = ET_2010_06# For experimental version, comment previous# three lines and uncomment next two# Note that you will need to get access# first (ask Erik Schnetter)#!TYPE = hg#!URL = ssh://carpetmercurial@carpetcode ...

.dyndns.org/carpet!CHECKOUT = Carpet/docCarpet/CarpetCarpet/CarpetEvolutionMaskCarpet/CarpetIOASCIICarpet/CarpetIOBasicCarpet/CarpetIOHDF5Carpet/CarpetIOScalar

84

B.3. Parametri iniziali

Carpet/CarpetInterpCarpet/CarpetInterp2Carpet/CarpetLibCarpet/CarpetMaskCarpet/CarpetReduceCarpet/CarpetRegridCarpet/CarpetRegrid2Carpet/CarpetSlabCarpet/CarpetTrackerCarpet/LoopControl

Whisky_Exp/WhiskyWhisky_Exp/Whisky_MagStarWhisky_Exp/Whisky_DissipationWhisky_Dev/Whisky_RNSIDWhisky_Dev/Whisky_TOVSolverCWhisky_Dev/Whisky_PerturbTOVWhisky_Dev/Whisky_ModePowerWhisky_Dev/Whisky_AnalysisWhisky_Dev/Whisky_HydroAnalysisWhisky_Exp/Whisky_Init_DataWhisky_Dev/Whisky_InitExcision

B.3 Parametri iniziali

Esempio di file di parametri utilizzato per una simulazione con pertur-bazione (file PerTOVA001L4d250.par):

# ’data/PerTOVA1L4d375.par’ automatically generated# by Cactus version 4.0.b17# Original parameter file was ’PerTOVA1L4d375.par’# Run time/date was 11:55:51+0200 Oct 25 2011 on# host ’grid-wn84.pr.infn.it’ with 8 processor(s)

ActiveThorns = "ADMBaseADMConstraintsADMCouplingADMMacrosAEILocalInterpBoundaryBSSN_MoLCarpetCarpetInterpCarpetIOASCIICarpetIOBasicCarpetIOHDF5CarpetIOScalarCarpetLibCarpetReduceCarpetRegrid2CarpetSlabCartGrid3DConstantsCoordBaseCoordGaugeEOS_BaseEOS_PolytropeFormalineFortranHydroBaseInitBaseIOUtil

85


LocalInterpLoopControlManualTerminationMoLNaNCheckerPsiKadeliaReflectionSymmetrySlabSpaceMaskSphericalHarmonicsSphericalSurfaceStaticConformalSummationByPartsSymBaseTimeTmunuBaseWaveExtractWhiskyWhisky_HydroAnalysisWhisky_Init_DataWhisky_PerturbTOVWhisky_TOVSolverC"

# Parameters of thorn ADMBase (implementing ADMBase)ADMBase::evolution_method = "adm_bssn"ADMBase::initial_data = "tov"ADMBase::initial_lapse = "tov"ADMBase::initial_shift = "tov"ADMBase::lapse_evolution_method = "1+log"ADMBase::metric_type = "physical"ADMBase::shift_evolution_method = "gamma0"

# Parameters of thorn ADMConstraints# (implementing admconstraints)ADMConstraints::constraints_persist = "yes"ADMConstraints::constraints_timelevels = 3

# Parameters of thorn ADMMacros# (implementing ADMMacros)ADMMacros::spatial_order = 4

# Parameters of thorn Boundary# (implementing boundary)Boundary::radpower = 3

# Parameters of thorn BSSN_MoL# (implementing adm_bssn)BSSN_MoL::advection = "upwind2"BSSN_MoL::AlphaDissip = 0.020000000000000000416BSSN_MoL::BetaDriver = 3.2000000000000001776BSSN_MoL::bound = "newrad"BSSN_MoL::harmonic_f = 2BSSN_MoL::LapsePsiPower = 2BSSN_MoL::lapsesource = "modified"BSSN_MoL::ShiftAlpPower = 2BSSN_MoL::ShiftGammaCoeff = 0.75BSSN_MoL::ShiftPsiPower = 4BSSN_MoL::stencil_size = 2BSSN_MoL::timelevels = 3

# Parameters of thorn Cactus# (implementing Cactus)Cactus::cctk_initial_time = 0Cactus::cctk_itlast = 26200

# Parameters of thorn Carpet# (implementing Driver)

86


Carpet::convergence_level = 0Carpet::domain_from_coordbase = "yes"Carpet::enable_all_storage = "no"Carpet::ghost_size = 3Carpet::init_each_timelevel = "no"Carpet::max_refinement_levels = 4Carpet::num_integrator_substeps = 4Carpet::prolongation_order_space = 3Carpet::prolongation_order_time = 2Carpet::regrid_during_initialisation = "no"Carpet::regrid_in_level_mode = "yes"Carpet::use_buffer_zones = "yes"

# Parameters of thorn CarpetIOASCII# (implementing IOASCII)CarpetIOASCII::out0D_dir = "data/0d"CarpetIOASCII::out0D_every = 16CarpetIOASCII::out0D_vars =

"whisky_hydroanalysis::hydro_total_rest_masswhisky_hydroanalysis::hydro_angular_momentumwhisky_hydroanalysis::hydro_center_masswhisky_hydroanalysis::hydro_quadrupolewhisky_HydroAnalysis::hydro_modes"

CarpetIOASCII::out1D_dir = "data/1d"CarpetIOASCII::out1D_every = 128CarpetIOASCII::out1D_vars = "admbase::metric

admbase::lapseadmbase::shifthydrobase::rhohydrobase::epshydrobase::presswhisky::denswhisky::tauwhisky::sconadmconstraints::momentumadmconstraints::hamPsiKadelia::WeylComponentsPsiKadelia::IJinvariants"

CarpetIOASCII::out2D_criterion = "divisor"CarpetIOASCII::out2D_dir = "data/2d"CarpetIOASCII::out2D_every = 16CarpetIOASCII::out2D_vars =

"SphericalHarmonics::decomposed_vars"

# Parameters of thorn CarpetIOBasic# (implementing IOBasic)CarpetIOBasic::outInfo_every = 1CarpetIOBasic::outInfo_vars =

" hydrobase::rho admbase::alp "

# Parameters of thorn CarpetIOHDF5#(implementing IOHDF5)CarpetIOHDF5::checkpoint = "yes"CarpetIOHDF5::out_dir = "data/HDF5"CarpetIOHDF5::out_every = 1024CarpetIOHDF5::out_vars =" hydrobase::rho compression_level=9 refinement_levels=3

hydrobase::presscompression_level=9 refinement_levels=3hydrobase::eps compression_level=9 refinement_levels=3hydrobase::vel compression_level=9 refinement_levels=3ADMBase::lapse compression_level=9 refinement_levels=3ADMBase::shift compression_level=9 refinement_levels=3ADMBase::curv compression_level=9 refinement_levels=3ADMBase::metriccompression_level=9 refinement_levels=3"

CarpetIOHDF5::use_reflevels_from_checkpoint = "yes"

# Parameters of thorn CarpetIOScalar

87


# (implementing IOScalar)CarpetIOScalar::outScalar_dir = "data/Scalar"CarpetIOScalar::outScalar_every = 8CarpetIOScalar::outScalar_reductions =

"minimum maximum norm1 norm2"CarpetIOScalar::outScalar_vars = "hydrobase::rho

hydrobase::presshydrobase::epsadmconstraints::hamadmconstraints::momentumadmbase::lapseadmbase::curvadmbase::metricadmbase::shift"

# Parameters of thorn CarpetRegrid2# (implementing CarpetRegrid2)CarpetRegrid2::num_centres = 1CarpetRegrid2::num_levels_1 = 4CarpetRegrid2::position_x_1 = 0CarpetRegrid2::radius_1[1] = 60CarpetRegrid2::radius_1[2] = 30CarpetRegrid2::radius_1[3] = 15

# Parameters of thorn CartGrid3D# (implementing grid)CartGrid3D::type = "coordbase"

# Parameters of thorn CoordBase# (implementing CoordBase)CoordBase::boundary_shiftout_z_lower = 1CoordBase::boundary_size_z_lower = 3CoordBase::domainsize = "minmax"CoordBase::dx = 2CoordBase::dy = 2CoordBase::dz = 2CoordBase::xmax = 120CoordBase::xmin = -120CoordBase::ymax = 120CoordBase::ymin = -120CoordBase::zmax = 120CoordBase::zmin = 0

# Parameters of thorn EOS_Polytrope# (implementing EOS_2d_Polytrope)EOS_Polytrope::eos_gamma = 2EOS_Polytrope::eos_k = 100

# Parameters of thorn HydroBase# (implementing HydroBase)HydroBase::evolution_method = "whisky"HydroBase::timelevels = 3

# Parameters of thorn InitBase# (implementing InitBase)InitBase::initial_data_setup_method =

"init_all_levels"

# Parameters of thorn IOUtil (implementing IO)IOUtil::checkpoint_dir = "CHECKPOINT"IOUtil::checkpoint_every = 1024IOUtil::checkpoint_ID = "no"IOUtil::checkpoint_keep = 2IOUtil::checkpoint_on_terminate = "yes"IOUtil::out_dir = "data"IOUtil::out_mode = "proc"IOUtil::out_single_precision = "yes"

88


IOUtil::out_unchunked = "no"IOUtil::parfile_write = "generate"IOUtil::recover = "autoprobe"IOUtil::recover_dir = "CHECKPOINT"

# Parameters of thorn ManualTermination# (implementing ManualTermination)ManualTermination::check_file_every = 16ManualTermination::max_walltime = 96ManualTermination::on_remaining_walltime = 120ManualTermination::termination_file = "terminate_job"ManualTermination::termination_from_file = "yes"

# Parameters of thorn MoL# (implementing MethodOfLines)MoL::initial_data_is_crap = "yes"MoL::MoL_Intermediate_Steps = 4MoL::MoL_Num_Scratch_Levels = 1MoL::ODE_Method = "rk4"

# Parameters of thorn NaNChecker# (implementing NaNChecker)NaNChecker::action_if_found = "terminate"NaNChecker::check_after = 6000NaNChecker::check_every = 128NaNChecker::check_vars = "admbase::gxx"

# Parameters of thorn PsiKadelia# (implementing GRscalars)PsiKadelia::psif_vec = "standard-radial"PsiKadelia::psikadelia_persists = "yes"PsiKadelia::ricci_prolongation_type = "none"PsiKadelia::weyl_timelevels = 3

# Parameters of thorn ReflectionSymmetry# (implementing ReflectionSymmetry)ReflectionSymmetry::avoid_origin_x = "no"ReflectionSymmetry::avoid_origin_y = "no"ReflectionSymmetry::avoid_origin_z = "no"ReflectionSymmetry::reflection_x = "no"ReflectionSymmetry::reflection_y = "no"ReflectionSymmetry::reflection_z = "yes"

# Parameters of thorn SpaceMask#(implementing SpaceMask)SpaceMask::use_mask = "yes"

# Parameters of thorn SphericalHarmonics# (implementing SphericalHarmonics)SphericalHarmonics::ex_radii[0] = 30SphericalHarmonics::ex_radii[1] = 35SphericalHarmonics::ex_radii[2] = 40SphericalHarmonics::ex_radii[3] = 45SphericalHarmonics::ex_radii[4] = 50SphericalHarmonics::ex_radii[5] = 55SphericalHarmonics::ex_radii[6] = 60SphericalHarmonics::ex_radii[7] = 65SphericalHarmonics::ex_radii[8] = 70SphericalHarmonics::ex_radii[9] = 75SphericalHarmonics::ex_radii[10] = 80SphericalHarmonics::ex_radii[11] = 85SphericalHarmonics::ex_radii[12] = 90SphericalHarmonics::ex_radii[13] = 95SphericalHarmonics::ex_radii[14] = 100SphericalHarmonics::ex_radii[15] = 105SphericalHarmonics::ex_radii[16] = 110SphericalHarmonics::ex_radii[17] = 115

89


SphericalHarmonics::grid_type = "cart3d"SphericalHarmonics::interp_integration_order = 4SphericalHarmonics::InterpPointsPhi = 156SphericalHarmonics::InterpPointsTheta = 49SphericalHarmonics::lmax = 6SphericalHarmonics::number_of_radii = 18SphericalHarmonics::number_of_vars = 2SphericalHarmonics::SH_spin_weight[0] = -2SphericalHarmonics::SH_spin_weight[1] = -2SphericalHarmonics::vars[0] = "PsiKadelia::psi4re"SphericalHarmonics::vars[1] = "PsiKadelia::psi4im"

# Parameters of thorn SphericalSurface# (implementing SphericalSurface)SphericalSurface::maxnphi = 40SphericalSurface::maxntheta = 21SphericalSurface::nghostsphi[0] = 2SphericalSurface::nghoststheta[0] = 2SphericalSurface::nphi[0] = 40SphericalSurface::nsurfaces = 6SphericalSurface::ntheta[0] = 13SphericalSurface::symmetric_z[0] = "yes"

# Parameters of thorn SummationByParts# (implementing SummationByParts)SummationByParts::order = 4

# Parameters of thorn Time# (implementing time)Time::dtfac = 0.25

# Parameters of thorn TmunuBase# (implementing TmunuBase)TmunuBase::prolongation_type = "none"TmunuBase::stress_energy_at_RHS = "yes"TmunuBase::stress_energy_storage = "yes"

# Parameters of thorn WaveExtract# (implementing WaveExtract)WaveExtract::detector_radius[0] = 30WaveExtract::detector_radius[1] = 35WaveExtract::detector_radius[2] = 40WaveExtract::detector_radius[3] = 45WaveExtract::detector_radius[4] = 50WaveExtract::detector_radius[5] = 55WaveExtract::detector_radius[6] = 60WaveExtract::detector_radius[7] = 65WaveExtract::detector_radius[8] = 70WaveExtract::detector_radius[9] = 75WaveExtract::detector_radius[10] = 80WaveExtract::detector_radius[11] = 85WaveExtract::detector_radius[12] = 90WaveExtract::detector_radius[13] = 95WaveExtract::detector_radius[14] = 100WaveExtract::detector_radius[15] = 105WaveExtract::detector_radius[16] = 110WaveExtract::detector_radius[17] = 115WaveExtract::interpolation_operator =

"Lagrange polynomial interpolation"WaveExtract::interpolation_order = 2WaveExtract::l_mode = 6WaveExtract::m_mode = 6WaveExtract::maximum_detector_number = 18WaveExtract::out_dir =

"data/WaveExtraction"WaveExtract::out_every = 16WaveExtract::rsch2_computation =

90


"average Schwarzschild metric"WaveExtract::switch_output_format = 21WaveExtract::verbose = 1

# Parameters of thorn Whisky# (implementing Whisky)Whisky::bound = "flat"Whisky::initial_atmosphere_factor =

0.10000000000000000555Whisky::ppm_detect = "yes"Whisky::ppm_omega2 = 10Whisky::recon_method = "ppm"Whisky::rho_rel_min =

9.9999999999999995475e-07Whisky::riemann_solver = "Marquina"Whisky::whisky_eos_table = "2D_Polytrope"Whisky::whisky_eos_type = "Polytype"Whisky::whisky_mhd_on = "no"Whisky::whisky_stencil = 3

# Parameters of thorn Whisky_HydroAnalysis# (implementing Whisky_HydroAnalysis)Whisky_HydroAnalysis::compute_hydro_angular_momentum = "yes"Whisky_HydroAnalysis::compute_hydro_center_of_mass = "yes"Whisky_HydroAnalysis::compute_hydro_modes = "no"Whisky_HydroAnalysis::compute_hydro_quadrupole = "yes"Whisky_HydroAnalysis::modes_radius_max = 15Whisky_HydroAnalysis::modes_radius_min

= 0.10000000000000000555Whisky_HydroAnalysis::use_com_as_origin = 0

# Parameters of thorn Whisky_PerturbTOV# (implementing Whisky_PerturbTOV)Whisky_PerturbTOV::aleph = 0Whisky_PerturbTOV::constraints = "yes"Whisky_PerturbTOV::l_mode = 2Whisky_PerturbTOV::matter_pert_type = "enthalpy"Whisky_PerturbTOV::Nnodes = 0.5Whisky_PerturbTOV::pert_amplitude =

0.0010000000000000000555Whisky_PerturbTOV::pert_type = "matter"Whisky_PerturbTOV::PertTOV_Num_Radial = 3000Whisky_PerturbTOV::TOV_Npoints_Phi = 50Whisky_PerturbTOV::TOV_Npoints_Theta = 50

# Parameters of thorn Whisky_TOVSolverC# (implementing WhiskyTOVSolver)Whisky_TOVSolverC::TOV_Combine_Method = "average"Whisky_TOVSolverC::TOV_dr[0] = 0.125Whisky_TOVSolverC::TOV_Enforce_Interpolation = "yes"Whisky_TOVSolverC::TOV_Gamma[0] = 2Whisky_TOVSolverC::TOV_K[0] = 100Whisky_TOVSolverC::TOV_Num_Radial = 3000Whisky_TOVSolverC::TOV_Num_TOVs = 1Whisky_TOVSolverC::TOV_Position_x[0] = 0Whisky_TOVSolverC::TOV_Position_y[0] = 0Whisky_TOVSolverC::TOV_Rho_Central[0] =

0.0012800000000000001047Whisky_TOVSolverC::TOV_Velocity_x[0] = 0Whisky_TOVSolverC::TOV_Velocity_y[0] = 0

91

BIBLIOGRAFIA

[1] Baiotti L., Bernuzzi S., Corvino G., De Pietri R., Nagar A., Gravitational-wave extraction from neutron-star oscillations: Comparing linear andnonlinear techniques Phys. Rev. D. 79 2009

[2] Miguel Alcubierre (2008), Introduction to 3+1 Numerical Relativity,Oxford, New York.

[3] Misner C., Thorne k.S. e Wheeler J.A. 1973, Gravitation W. H. Freemanand Company, San Francisco.

[4] Regge T. e Wheeler John T., Stability of a Schwarzschild Singularity.Phys. Rev., 108:1063-1069, 1957.

[5] Zerilli F.J., Effective potential for even-parity Regge-Wheeler gra-vitational perturbation equations. Phys. Rev. Lett., 24(13):737-738,1970.

[6] Zerilli F.J., Gravitational field of a particle falling in a Schwarzschildgeometry analyzed in tensor harmonics. Phys. Rev. D, 2:2141, 1970.

[7] Martel K. e Poisson E., Gravitational perturbations of the schwarzschildspacetime: A practical covariant and gauge-invariant formalism. Phys.Rev., D71:104003, 2005.

[8] Alessandro Nagar e Luciano Rezzolla, Gauge-invariant non-sphericalmetric perturbations of Schwarzschild black-hole spacetime. Class.Quantum Grav. 22 R167-R192, 2005

[9] Moncrief V., Gravitational perturbations of spherically symmetricsystems. I. the exterior problem. Annals of physics, 88:323-342, 1974

[10] Isaacson R., Gravitational radiation in the limit of high frequen-cy.II.nonlinear terms and the effective stress tensor. Phys. Rev.,166:1272-1280, 1968.

[11] Baumgarte, T.W. e Shapiro, S.L., On the numerical integration ofEinstein’s field equations. Phys. Rev. D, 59:024007, 1999

93

Bibliografia

[12] Shibata, M. e Nakamura, T., Evolution of three-dimensionalgravitational waves: Harmonic slicing case. Phys. Rev. D, 52:5428,1995.

[13] Gerlach, U. e Sengupta, U., Gauge-invariant perturbations on mostgeneral spherically symmetric space-times. Phys. rev. D., 19:2268-2272,1979.

[14] Arnowitt, R., Deser, S., e Misner, C. W. The dynamics of general relati-vity. Witten, L., editore, Gravitation: An introduction to current research,pag. 227-265. John Wiley, New York, 1962.

[15] York, J. W. Kinematics and dynamics of general relativity. In Smarr, L.L., editore, Sources of gravitational radiation, pagine 83-126. CambridgeUniversity Press, Cambridge, UK, 1979.

[16] Smarr, L. e York J. W. Kinematical conditions in the construction ofspacetime. Phys. Rev. D, 17(10):2529-2552, 15 Maggio 1978.

[17] Smarr, L. e York J. W. Radiation gauge in general relativity. Phys. Rev.D, 17:1945, 1978.

[18] Bona, C., Massoo, J., Seidel, E., Stela, J. First order hyperbolicformalism for numerical relativity. Phys. Rev. D, 56:3405-3415, 1995.

[19] Hartle, J. B. Gravity an introduction to Einstein’s General Relativity.

[20] Dirac, P. A. M. Fixation of coordinates in the hamiltonian theory ofgravitation. Phys. Rev., 114:924, 1959.

[21] Banyuls F., Font J. A., Ibanez J. M., Marti J., Miralles J. A. Astrophys. J.476, 221 (1997).

[22] Gunnarsen L., Shinkai H. A., Maeda I. K., Classical Quantum Gravity12, 133, 1995.

[23] Abrahams A. M., Price R. H., Phys. Rev. D 53, 1963.

[24] Stergioulas N., Apostolatos T. A., Font J. A., Mon. Not. R. Astron. Soc.352, 1089, 2004.

[25] Nagar A., Diaz G., in Proceedings of the 27th Spanish Relativity Meeting(ERE 2003): Gravitational Radiation, Alicante, Spain, 2003 (University ofAlicante, Alicante, Spain 2004).

[26] Newman E. T., Penrose R., Note on the Bondi-Metzner-Sachs group. J.Math. Phys., 7(5):863-870, May 1996

[27] CACTUS Code Homepage. http://www.cactuscode.org/.

[28] CARPET Code Homepage. http://www.carpetcode.org/.

[29] Einstein Toolkit Homepage. http://www.einsteintoolkit.org/.

[30] WHISKY Code Homepage. http://www.whiskycode.org/.

94

http://www.cactuscode.org/

http://www.carpetcode.org/

http://www.einsteintoolkit.org/

http://www.einsteintoolkit.org/

http://www.whiskycode.org/

Bibliografia

95

NIVERSITA DEGLI ` STUDI DI PARMA - infn.it delle componenti del tensore di curvatura di Weyl...

Documents

Transcript of NIVERSITA DEGLI ` STUDI DI PARMA - infn.it delle componenti del tensore di curvatura di Weyl...