Szacowanie masy w stacjonarnej hydrodynamice ogólnorelatywistycznej

Uniwersytet Jagielloński w Krakowie

Wydział Fizyki, Astronomii i Informatyki Stosowanej

Wojciech Kulczycki

Szacowanie masy w stacjonarnej hydrodynamice ogólnorelatywistycznej

Rozprawa doktorska

Praca wykonana pod kierunkiem prof. dra hab. Edwarda Malca

Kraków 2019

Wydział Fizyki, Astronomii i Informatyki Stosowanej Uniwersytet Jagielloński

Oświadczenie

Ja niżej podpisany Wojciech Kulczycki (nr indeksu: 1074691) doktorant Wydziału Fizyki, Astronomii i Informatyki Stosowanej Uniwersytetu Jagiellońskiego oświadczam, że przedłożona przeze mnie rozprawa doktorska pt. „Szacowanie masy w stacjonarnej hydrodynamice ogólnore- latywistycznej” jest oryginalna i przedstawia wyniki badań wykonanych przeze mnie osobiście, pod kierunkiem prof. dra hab. Edwarda Malca. Pracę napisałem samodzielnie.

Oświadczam, że moja rozprawa doktorska została opracowana zgodnie z Ustawą o prawie autorskim i prawach pokrewnych z dnia 4 lutego 1994 r. (Dziennik Ustaw 1994 nr 24 poz. 83 wraz z późniejszymi zmianami).

Jestem świadom, że niezgodność niniejszego oświadczenia z prawdą ujawniona w dowol- nym czasie, niezależnie od skutków prawnych wynikających z ww. ustawy, może spowodować unieważnienie stopnia nabytego na podstawie tej rozprawy.

Kraków, dnia 21 marca 2019 r. ...

podpis doktoranta

Pragnę serdecznie podziękować mojemu promotorowi i opiekunowi naukowemu, profesorowi Edwardowi Malcowi, za pomoc, wsparcie i poświęcony mi czas w trakcie pisania niniejszej rozprawy, jak również przez całe studia na Uniwersytecie Jagiellońskim.

Dziękuję również doktorowi habilitowanemu Patrykowi Machowi za nieocenione rady i walny wkład w stworzenie kodu numerycznego, na którym opierają się obliczenia przedstawione w niniejszej rozprawie.

Swoje podziękowania pragnę złożyć także na ręce wszystkich członków grupy badawczej zaj- mującej się badaniem rozważanych w niniejszej rozprawie układów, a której mam przyjemność być członkiem. Do zespołu tego należą: prof. dr hab. Edward Malec, dr hab. Janusz Karkowski, dr hab. Andrzej Odrzywołek, dr hab. Patryk Mach, dr Michał Piróg, a od niedawna również mgr Wojciech Dyba. Dziękuję także wszystkim pracownikom Instytutu Fizyki Uniwersytetu Jegiellońskiego, a w szczególności członkom Zakładu Teorii Względności i Astrofizyki.

Pragnę również podziękować moim nauczycielom ze szkoły średniej, a w szczególności ma- gister Bożenie Szalewicz oraz magistrowi Stanisławowi Lisowi, którym zawdzięczam wybór kie- runku studiów i przygotowanie do ich podjęcia.

Serdeczne podziękowania kieruję również ku swoim najbliższym — Rodzicom oraz Siostrze

— za ich niezłomne wsparcie i pomoc.

Spis treści

1. Wstęp . . . 9

Oznaczenia i konwencje . . . 11

2. Konstrukcja numeryczna badanych układów ogólnorelatywistycznych . . . 15

2.1. Opis układu . . . 15

2.2. Metryka i równania Einsteina . . . 15

2.3. Prawo rotacji . . . 18

2.3.1. Wyprowadzenie prawa rotacji dla testowego dysku pyłowego . . . 19

2.3.1.1. Sposób I . . . 19

2.3.1.2. Sposób II . . . 21

2.3.2. Prawo rotacji dla ważkiego, samograwitującego dysku . . . 22

2.4. Symetrie równań . . . 22

2.4.1. Symetria inwersji . . . 22

2.4.2. Symetria skalowania . . . 23

2.5. Opis metody numerycznej . . . 23

2.5.1. Siatka . . . 24

2.5.2. Warunki brzegowe . . . 24

2.5.3. Dyskretyzacja równań . . . 25

2.5.4. Dane początkowe . . . 25

2.5.5. Procedura numeryczna . . . 26

2.5.6. Skalowanie . . . 27

2.6. Wielkości globalne . . . 27

3. Szacowanie masy układów newtonowskich z dyskiem . . . 30

3.1. Rotacja keplerowska w osiowosymetrycznych układach newtonowskich . . . 30

3.2. Opis układu . . . 30

3.3. Nierówności ograniczające masę . . . 31

3.3.1. Motywacja . . . 31

3.3.2. Przypadki szczególne . . . 31

3.3.2.1. Przybliżenie dysku testowego . . . 31

3.3.2.2. Granica $₁ → 0 przy ustalonych wartościach $₂, ρmax6= 0 . . . 31

3.3.2.3. Nieskończony dysku ze skończoną masą . . . 32

3.3.2.4. Granica $₁ → $₂ przy ustalonych wartościach $₂, ρmax6= 0 . . . 32

3.3.3. Hipoteza uogólniająca . . . 33

3.4. Skalowanie . . . 33

3.5. Przegląd wyników numerycznych . . . 34

3.5.1. Nierówności ograniczające masę . . . 34

3.5.2. Graniczna wartość ˜ω0 . . . 35

3.5.3. Graniczna wartość ˜M_d . . . 35

3.6. Wnioski . . . 35

4. Szacowanie masy układów w czasoprzestrzeni Kerra z pyłowym dyskiem testowym 37 4.1. Informacje wstępne . . . 37

4.2. Nierówności ograniczające masę . . . 38

4.2.1. a = 0 . . . 38

4.2.2. a > 0 . . . 38

4.2.3. a < 0 . . . 40

4.3. Nowa nierówność . . . 41

4.4. Nierówność uproszczona . . . 41

4.5. Wnioski . . . 41

5. Ogólnorelatywistyczne układy samograwitujące — szacowanie parametru rotacji 42 5.1. Hipotezy . . . 42

5.2. Parametr ˜w oraz masa dysku ˜M_T dla ˜M_h/ ˜M_T 1 . . . 43

5.3. Wnioski . . . 47

6. Masa, kręt i prędkość kątowa w układach samograwitujących . . . 48

6.1. Hipotezy . . . 48

6.2. Badanie wartości wyrażenia Ωr_c^3/2/√ M . . . 49

6.3. Wnioski . . . 51

7. Dalsze możliwości badań . . . 54

7.1. Masa dysku a gęstość maksymalna . . . 54

7.1.1. Motywacja . . . 54

7.1.2. Związek gęstości maksymalnej z masą dysku . . . 54

7.1.3. Maksimum gęstości ρ_max . . . 55

7.1.4. Maksimum promienia obwodowego w dysku . . . 56

7.1.5. Nierówności (5.1) i (6.1) dla bardzo ciężkich dysków . . . 56

7.2. Dyski kontrrotujące z brzegiem wewnętrznym pod ISCO . . . 58

7.3. Problem stabilności . . . 60

7.4. Wnioski . . . 61

8. Wnioski i podsumowanie . . . 62

A. Wyprowadzenie równań Einsteina (2.15) . . . 64

A.1. Wybrane elementy formalizmu 3 + 1 . . . 64

A.2. Transformacja konforemna równań więzów . . . 66

A.3. Równania Einsteina w metryce (2.1) . . . 67

A.3.1. Więz pędowy . . . 68

A.3.2. Więz hamiltonowski . . . 69

A.3.3. Równania ewolucji . . . 69

A.3.4. Przekształcenia układ równań . . . 72

A.3.5. Aplikacja formalizmu „puncture” . . . 73

B. Separacja zmiennych w nierównościach z parametrem rotacji . . . 74

Bibliografia . . . 76

1. Wstęp

Niech dany będzie związany grawitacyjnie, stacjonarny układ złożony z czarnej dziury i gazowego dysku, dla którego obserwator może zmierzyć rozmiary liniowe i, poprzez pomiar przesunięcia widma promieniowania ku czerwieni/fioletowi, wyznaczyć krzywą rotacji, tj. pręd- kość kątową płynu. Przy badaniu takich obiektów pojawia się naturalne pytanie, ile wynosi masa układu i/lub jego składowych. Niniejsza rozprawa poświęcona jest próbie odpowiedzi na to pytanie.

Przechodząc do szczegółów, celem niniejszej rozprawy jest analiza ogólnorelatywistycznych układów zbudowanych z czarnej dziury oraz ważkiego, samograwitującego, toroidalnego dysku rotującego keplerowsko pod kątem szacowania masy w takim systemie. Omówione zostaną także przypadki pyłowego dysku testowego w czasoprzestrzeni Kerra oraz układów newtonowskich.

Te drugie były przedmiotem prac [1, 2], a w niniejszej rozprawie wyniki wcześniej uzyskane zostaną uzupełnione o nowe elementy.

Układy złożone z czarnej dziury i rotującego dysku akrecyjnego są powszechnie obserwo- wane w aktywnych jądrach galaktyk. W wielu przypadkach z danych obserwacyjnych udało się wyznaczyć w nich krzywą rotacji materii. Jedną z najlepiej poświadczonych jest rotacja keplerowska [3], czego przykładem jest galaktyka NGC 4258 [4, 5]. Masę takiego układu, bądź też centralnej czarnej dziury można wyznaczyć z danych obserwacyjnych tylko w bardzo szcze- gólnych przypadkach, gdy obserwowane są ciała próbne na orbitach wokół niego. Przykładami takich obiektów są Sag A* [6], czy wspomniana już galaktyka NGC 4258 [5]. Motywacją do podjęcia zaprezentowanych w niniejszej rozprawie badań była próba odpowiedzi na pytanie, czy na podstawie wielkości obserwowanych związanych z dyskiem akrecyjnym, takich jak pręd- kość kątowa materii oraz odległość od masy centralnej można oszacować masę układu lub też czarnej dziury.

W niniejszej rozprawie badane są stacjonarne konfiguracje złożone z czarnej dziury oraz toroidalnego, samograwitującego, osiowosymetrycznego dysku akrecyjnego opisanego politropo- wym równaniem stanu i rotującym według ogólnorelatywistycznego uogólnienia keplerowskiego prawa rotacji. Rozprawa składa się z sześciu rozdziałów.

W rozdziale 2 omówiona została konstrukcja numeryczna badanych ogólnorelatywistycznych układów. Przedstawione zostały w niej równania Einsteina, jak również wyprowadzenia prawa rotacji. Zaprezentowano też użytą procedurę numeryczną oraz zdefiniowano potrzebne wielkości.

Rozdział 3 poświęcony jest szacowaniu masy w układach newtonowskich, złożonych z punk- towej masy centralnej oraz politropowego, samograwitującego, osiowosymetrycznego dysku ro- tującego keplerowsko. Bazuje on na wynikach przedstawionych w pracach [1, 2], uzupełniając je o nowe elementy związane z problemem szacowania masy.

W rozdziale 4 analizowane są pyłowe dyski testowe w czasoprzestrzeni Kerra. Zaprezento- wano w nim analityczny dowód nierówności ograniczających masę układu oraz masę nieredu- kowalną czarnej dziury, który w krótszej formie został opublikowany w pracy [30].

Rozdziały 5 i 6 dotyczą problemu szacowania masy w przypadku układów ogólnorelatywi- stycznych z ważkim, samograwitującym dyskiem akrecyjnym. Przedstawione są także wyniki numeryczne wspierające postawione hipotezy. Po części wyniki te były opublikowane w pracy [30].

W rozdziale 7 przedstawiono wyniki, które nie były jeszcze publikowane. Nakreśla się także dalsze perspektywy badań.

W zamieszczonym na końcu dodatku A zaprezentowano wyprowadzenie równań Einsteina opisujących układ ogólnorelatywistyczny będący przedmiotem niniejszej rozprawy. Z kolei w dodatku B uzasadniono możliwość separacji zmiennych w dwóch przedstawionych w rozdziale 5 nierównościach.

Analiza układów newtonowskich wskazuje, że masę centralną oraz masę dysku można ogra- niczyć od góry i od dołu parametrem pojawiającym się w prawie rotacji, maksymalną wartością gęstości barionowej materii w dysku oraz rozmiarem liniowym układu. W układach ogólnore- latywistycznych uzyskuje się analogiczne ograniczenia, wspierane przez wyniki numeryczne.

Analiza pyłowych dysków testowych wykazała istnienie dodatkowych ograniczeń masy układu, udowodnionych analitycznie. Wyniki numeryczne wskazują na prawdziwość uogólnień tych nie- równości na przypadek ważkich dysków przy dwóch dodatkowych założeniach.

Literatura na temat płynów rotujących wokół masy centralnej jest obszerna. Poniżej przed- stawione zostaną wybrane prace dotyczące tego zagadnienia. W przypadku układów newtonow- skich rozważano zarówno przypadki dysku testowego, np. w pozycji [4], jak i samograwitującego, m.in. artykułach [7, 8, 9]. Analiza układów złożonych z politropowego dysku rotującego keple- rowsko była przedmiotem prac [1, 2]. Zaprezentowano w nich także nierówności ograniczające masę centralną oraz masę dysku. Badane były również układy ogólnorelatywistyczne. W li- teraturze istnieje kilka podejść do tego zagadnienia. Na ogół rozważane są dyski testowe w ustalonej czasoprzestrzeni, np. [10] oraz przeglądy [11, 12, 13] z licznymi referencjami. Przypa- dek czasoprzestrzeni Kerra omawiany jest m.in. w pozycji [14]. Do pierwszych prac, w których w obliczeniach numerycznych uwzględniono samograwitację należą [15, 16, 17]. Odtworzono wówczas tzw. efekt wleczenia grawitacyjnego. O ile w przypadku teorii Newtona twierdzenie Poincarégo-Wavre’a pozwala zdefiniować szeroką klasę prostych praw rotacji, pośród których znajduje się rotacja keplerowska, to w ogólnej teorii względności poznano ich niewiele. Standar- dowo używano rotacji sztywnej [18, 19, 20], bądź też ze stałą gęstością krętu [21]. W tej ostatniej pozycji po raz pierwszy w konstrukcji numerycznej użytej do modelowania układów z ważkim, samograwitującym dyskiem zastosowano formalizm „puncture” [22]. Spośród alternatywnych sformułowań warto wspomnieć pionierskie podejście Nishidy i Eriguchiego [16] oraz tzw. me- todę pseudospektralną [20]. Analizę dysków keplerowskich w przybliżeniu postnewtonowskim zaprezentowano w pracach [23, 24]. Odkryto wówczas przeciwny do wleczenia grawitacyjnego efekt nazwany antywleczeniem grawitacyjnym. W 2012 r. odkryte zostało nieliniowe prawo rotacji [25], które w granicy newtonowskiej może w przybliżeniu odtwarzać klasę krzywych rotacji z prędkością kątową proporcjonalną do potęgi promienia cylindrycznego. Należy do niej także rotacja keplerowska. Prawo rotacji, które w granicy newtonowskiej w sposób dokładny odtwarza tę klasę potęgowych krzywych rotacji odkryto w 2015 r. [26]. Stosuje się ono w szcze- gólności do ogólnorelatywistycznej rotacji keplerowskiej wokół bezspinowych czarnych dziur.

Jego numeryczną implementację w przybliżeniu postnewtonowskim przeprowadzono w pozycji [27]. W opisie w pełni ogólnorelatywistycznym dokonano tego w pracach z moim udziałem [28, 29]. Pozycje te zawierają również dwa nowe elementy. Najistotniejszym jest uzyskane w nich nowe prawo rotacji — ogólnorelatywistyczna wersja rotacji keplerowskiej wokół wirującej czarnej dziury. Ponadto, dokonano również jego numerycznej implementacji. To nowe prawo rotacji, w uzyskaniu którego brałem aktywny udział, i metoda numeryczna [28, 29] stały się podstawą do badań zrelacjonowanych w napisanej z moim udziałem pracy [30]. Przedstawiono w niej nierówności wiążące prędkość kątową, masę, kręt i rozmiar układu złożonego z wirującej czarnej dziury i dysku.

Oznaczenia i konwencje

Układy współrzędnych

— Ogólnorelatywistyczne:

◦ Współrzędne quasi-izotropowe (metryka g_µν): (t, r, θ, φ) (str. 15)

◦ Współrzędne Boyera-Lindquista (metryka g_µν^BL): (t, rBL, θ, φ) (tBL= t, θBL= θ, φBL = φ) (str. 17)

— Newtonowskie:

◦ Współrzędne cylindryczne: ($, φ, z)

◦ Współrzędne sferyczne: (r, θ, φ) Oznaczenia użyte w pracy

— Stałe fizyczne:

◦ c — prędkość światła w próżni

◦ G — stała grawitacji

— Wielkości wspólne dla układów newtonowskich i ogólnorelatywistycznych:

◦ ρ — gęstość barionowa materii dysku (str. 15)

◦ p — ciśnienie gazu w dysku (str. 15)

◦ ρ_max — maksymalna wartość gęstości barionowej ρ w dysku (str. 26)

◦ p_max — maksymalna wartość ciśnienia p w dysku (str. 55)

◦ K — stała proporcjonalności w politropowym równaniu stanu (str. 15)

◦ γ — wykładnik adiabatyczny w politropowym równaniu stanu (str. 15)

◦ h — entalpia właściwa (str. 15)

◦ Ω — prędkość kątowa (str. 16)

— Wielkości opisujące układ ogólnorelatywistyczny:

◦ m — w czasoprzestrzeni Kerra masa czarnej dziury (str. 17); w obecności samograwitu- jącego dysku parametr masy (str. 17)

◦ a — w czasoprzestrzeni Kerra spin czarnej dziury (str. 17); w obecności samograwitują- cego dysku parametr spinu (str. 17)

◦ ˆa — rzeczywisty spin czarnej dziury; w czasoprzestrzeni Kerra a = ˆa (str. 29)

◦ a_rot — parametr spinu w prawie rotacji; wprowadzony na potrzeby obliczeń numerycz- nych (str. 26)

◦ α, β, ψ, q — funkcje metryczne w metryce g_µν (str. 15)

◦ α_K, β_K, ψ_K, q_K — funkcje metryczne dla czasoprzestrzeni Kerra w metryce g_µν (str. 17);

β_K to również wkład do funkcji β pochodzący od czarnej dziury (str. 15)

◦ β_T — wkład do funkcji β pochodzący od dysku (str. 15)

◦ K_ij — krzywizna zewnętrzna, zwana też II formą fundamentalną (str. 16)

◦ Kˆ_ij — krzywizna zewnętrzna po transformacji konforemnej metryki (str. 16)

◦ H_E, H_F — funkcje proporcjonalne odpowiednio do składowych ˆK_rφoraz ˆK_θφtensora ˆK_ij (str. 16)

◦ Σ_BL, ∆_BL — funkcje pomocnicze w metryce g_µν^BL (str. 17)

◦ ϕ, B — przedefiniowane funkcje metryczne ψ, α (str. 18)

◦ r_c — (łac.: circumferentialis) promień obwodowy (str. 15)

◦ S_q, S_ϕ, S_B, S_βT — źródła w równaniach Einsteina (str. 18)

◦ A, b — wielkości pomocnicze występujące w źródłach S_q, S_ϕ, S_B, S_βT (str. 18)

◦ T_µν — tensor energii-pędu (str. 15)

◦ Ω₁ — wartość prędkości kątowej Ω na równiku na wewnętrznym brzegu dysku (str. 26)

◦ Ω₂ — wartość prędkości kątowej Ω na równiku na zewnętrznym brzegu dysku (str. 26)

◦ C₁, C₂ — stałe całkowania w dwóch wariantach równania Bernoulliego (str. 19)

◦ C — stała całkowania po przeformułowaniu równania Bernoulliego (str. 19)

◦ rh — promień horyzontu czarnej dziury (str. 17)

◦ r_BL,h — promień horyzontu czarnej dziury we współrzędnych Boyera-Lindquista (str.

37)

◦ rISCO — promień najbardziej wewnętrznej stabilnej orbity kołowej (str. 37)

◦ r_BL,ISCO — promień najbardziej wewnętrznej stabilnej orbity kołowej we współrzędnych Boyera-Lindquista (str. 38)

◦ rc,ISCO— promień najbardziej wewnętrznej stabilnej orbity kołowej wyrażony w promie- niu obwodowym (str. 38)

◦ r_c,ISCO,K — promień najbardziej wewnętrznej stabilnej orbity kołowej w czasoprzestrzeni Kerra wyrażony w promieniu obwodowym (str. 42)

◦ Z₁, Z₂ — funkcje pomocnicze przy definiowaniu r_ISCOw czasoprzestrzeni Kerra (str. 38)

◦ r, ¯¯ ρ, . . . — zmienne po transformacji inwersji (str. 22)

◦ r, ˜˜ ρ, . . . — zmienne po transformacji skalowania (str. 23)

◦ λ — parametr skalujący (str. 23)

◦ r_in — (łac.: internus) promień wewnętrznego brzegu siatki numerycznej (str. 24)

◦ rex — (łac.: externus) promień zewnętrznego brzegu siatki numerycznej (str. 24)

◦ N_r, N_θ — liczba punktów siatki numerycznej odpowiednio w kierunku radialnym i kąto- wym (str. 24)

◦ ri (i = 1, . . . , Nr) — współrzędna radialna węzłów siatki numerycznej (str. 24)

◦ θ_j (j = 1, . . . , N_θ) — współrzędna kątowa węzłów siatki numerycznej (str. 24)

◦ f, ∆µ — stałe użyte przy konstrukcji siatki numerycznej (str. 24)

◦ M1, B1, J1, q1 — funkcje pojawiające się w asymptotyce funkcji metrycznych φ, B, βT, q (str. 25)

◦ r₁ — promień współrzędnościowy wewnętrznego brzegu dysku (str. 25)

◦ r2 — promień współrzędnościowy zewnętrznego brzegu dysku (str. 25)

◦ r_c,1 — promień obwodowy geometrycznego wewnętrznego brzegu dysku (str. 29)

◦ r_c,2 — promień obwodowy geometrycznego zewnętrznego brzegu dysku (str. 29)

◦ i2 — numer węzła siatki zewnętrznego brzegu dysku (str. 27)

◦ θ_sec — (łac. sectio) kąt ucięcia: wartość graniczna kąta θ przy wyliczaniu prędkości kątowej Ω; wprowadzony na potrzeby obliczeń numerycznych (str. 26)

◦ u^µ — 4-prędkość (str. 16)

◦ V — prędkość liniowa (str. 19)

◦ Ω_w — prędkość wleczenia (str. 37)

◦ A˜1, ˜A2 — wyrażenia pomocnicze w nierównościach ograniczających masę (str. 48)

◦ j — gęstość krętu (str. 18)

◦ ξ, X, Z, N , D, N₁, N₂, D₁, D₂, D₃ — wielkości pomocnicze przy wyprowadzeniu prawa rotacji dla dysków testowych w czasoprzestrzeni Kerra (str. 20)

◦ M — masa asymptotyczna układu (str. 28)

◦ M_h — masa czarnej dziury (str. 28)

◦ M_T — masa dysku (str. 29)

◦ J — kręt asymptotyczny układu (str. 28)

◦ J_h — kręt czarnej dziury (str. 28)

◦ J_T — kręt dysku (str. 28)

◦ A_h — powierzchnia horyzontu czarnej dziury (str. 28)

◦ M_irr — (łac.: irreducibilis) masa nieredukowalna czarnej dziury (str. 28)

◦ w — stała w ogólnorelatywistycznym keplerowskim prawie rotacji (str. 19)

◦ fˆ_h, ˆf_T — funkcje pomocnicze w nierównościach ograniczających odpowiednio M_h i M_T (str. 45)

◦ f_h, f_T — zależne jedynie od wartości γ części odpowiednio funkcji ˆf_h, ˆf_T (str. 45)

— Wielkości opisujące układ newtonowski:

◦ ω0 — stała proporcjonalności w newtonowskim keplerowskim prawie rotacji (str. 30)

◦ U — potencjał grawitacyjny (str. 30)

◦ U_d — potencjał grawitacyjny pochodzący od dysku (str. 30)

◦ Ud,min — minimalna wartość potencjału Ud w obrębie dysku (str. 32)

◦ U_c — potencjał odśrodkowy (str. 31)

◦ C — stała całkowania w równaniu Bernoulliego (str. 31)˜

◦ v — pole prędkości cząstek w dysku (str. 30)

◦ $₁, $₂— promień cylindryczny odpowiednio wewnętrznego i zewnętrznego brzegu dysku (str. 31)

◦ Mc — masa centralna (str. 30)

◦ M_d — masa dysku (str. 30)

◦ r, ˜˜ $, ˜z, ˜$₁, ˜$₂, ˜ρ, ˜h, ˜U , ˜U_d, ˜U_c, ˜M_c, ˜M_d, ˜ω₀ — przeskalowane wielkości r, $, z, $₁, $₂, ρ, h, U , Ud, Uc, Mc, Md, ω0 (str. 33)

◦ ∆ — laplasjan w przeskalowanych współrzędnych (str. 33)˜

◦ f_c, f_d — funkcje pomocnicze w nierównościach ograniczających odpowiednio M_c i M_d (str. 36)

Przyjęte konwencje:

— Indeksy greckie przyjmują wartości 0, 1, 2, 3.

— Indeksy łacińskie przebiegają wartości 1, 2, 3.

— Sygnatura tensora metrycznego g_µν: (−, +, +, +).

— Przy analizie układów ogólnorelatywistycznych: G = c = 1.

— Dla dysku kontrrotującego względem czarnej dziury: Ω > 0, a < 0.

2. Konstrukcja numeryczna badanych układów ogólnorelatywistycznych

2.1. Opis układu

Rozważany układ składa się z czarnej dziury oraz rotującego keplerowsko, toroidalnego dysku. Zakłada się także jego stacjonarność oraz osiową symetrię. W konstrukcji uwzględniona została grawitacja pochodząca od materii, dzięki czemu można badać również ważkie dyski.

Niniejszy rozdział bazuje na pracach [28, 29].

2.2. Metryka i równania Einsteina

W konstrukcji badanego układu użyto stacjonarnej quasi-izotropowej metryki danej wzorem:

g_µνdx^µdx^ν = −α²dt²+ ψ⁴e^2q dr²+ r²dθ²
+ r²sin²θ (dφ + βdt)² , (2.1) gdzie funkcje metryczne α, β, ψ, q w przypadku osiowej symetrii zależą jedynie od współrzęd- nych r i θ, α jest funkcją upływu czasu, tzw. „lapsem” (ang.: lapse), β składową φ wektora przesunięcia (ang.: shift vector), a ψ czynnikiem konforemnym.

Współrzędna radialna r nie jest wielkością geometryczną, w związku z czym nie jest ob- serwablą. Dlatego też z punktu widzenia obserwacji astronomicznych znaczenie ma promień obwodowy rc zdefiniowany następująco:

r_c≡ q

|gφφ| = ψ²r sin θ.

Funkcję β można przedstawić następująco:

β = β_K+ β_T, (2.2)

gdzie βK i βT są częściami pochodzącymi odpowiednio od czarnej dziury oraz dysku.

Przyjęto, że materię opisuje tensor energii pędu gazu doskonałego dany wzorem:

T_µν = ρhu_µu_ν + pg_µν, (2.3)

gdzie ρ, h, p są odpowiednio gęstością barionową, entalpią właściwą i ciśnieniem. W literatu- rze częściej spotyka się postać tensora T_µν, w której w miejscu iloczynu ρh pojawia się suma ciśnienia oraz gęstości energii. Obydwie te formy są równoważne w przypadku barotropowych równań stanu.

W badanych układach przyjęto politropowe równanie stanu o postaci:

p = Kρ^γ, (2.4)

gdzie K i γsą stałymi.

Wówczas entalpia właściwa h wyraża się następująco:

h = 1 + Kγ γ − 1ρ^γ−1,

przy czym przyjęto tu konwencję, zgodnie z którą w próżni (tj. dla ρ ≡ 0) h ≡ 1.

W rozważanym układzie 4-prędkość cząstek dysku dana jest wzorem:

u^µ= u^t, 0, 0, u^φ

(2.5) i jest ona znormalizowana:

u^µu_µ= −1. (2.6)

Prędkość kątowa Ω zdefiniowana jest jako:

Ω ≡ u^φ

u^t. (2.7)

Z równań (2.5), (2.6) i (2.7) uzyskuje się:

u^t = −g_tt− 2Ωg_tφ− Ω²g_φφ
^−1/2 , u^φ = Ω −g_tt− 2Ωg_tφ− Ω²g_φφ
−1/2

.

W celu wyprowadzenia równań Einsteina dla danych początkowych użyto formalizmu 3 + 1 ogólnej teorii względności, wykorzystując tzw. metodę „puncture”, sformułowaną pierwotnie w celu modelowania układów wielu czarnych dziur [22]. Równania te zostały wyprowadzone przez Masaru Shibatę w pracy [21]. Metoda ta gwarantuje odtworzenie czasoprzestrzeni Kerra w granicy zerowej masy dysku, a także umożliwia umieszczenie wewnętrznego brzegu siatki numerycznej na horyzoncie czarnej dziury. Szczegóły wyprowadzenia równań Einsteina przed- stawiono w dodatku A. Poniżej ograniczono się do ich przedstawienia.

Niech K_ij oznacza krzywiznę zewnętrzną na hiperpowierzchni stałego czasu (t = const).

Niech ˆK_ij będzie krzywizną zewnętrzną po transformacji konforemnej metryki:

Kˆ_ij ≡ ψ²K_ij.

W przypadku czasoprzestrzeni osiowosymetrycznej opisanej metryką (2.1) jedynymi niezero- wymi składowymi ˆK_ij są:

Kˆ_rφ = ψ⁶

2αr²sin²θ ∂_rβ, Kˆ_θφ = ψ⁶

2αr²sin²θ ∂_θβ.

Wówczas jedyne nietrywialne równanie więzu pędowego przyjmuje postać:

∂_r

r²Kˆ_rφ

r² +

∂_θ

sin θ ˆK_θφ

r²sin θ = 8παψ⁶e^2qT_φ^t. (2.8) Z uwagi na liniowość tego równania względem krzywizny ˆK_ij, można je rozseparować, ko- rzystając z równania (2.2), na części związane z β_K oraz β_T:

Kˆ_rφ = H_Esin²θ r² + ψ⁶

2αr²sin²θ ∂_rβ_T, Kˆ_θφ = H_Fsin θ

r + ψ⁶

2αr²sin²θ ∂_θβ_T,

gdzie funkcje H_E i H_F spełniają równanie uzyskane z więzu pędowego (2.8) [31, 32]:

r∂H_E

∂r sin³θ + ∂ H_Fsin²θ

∂θ = 0.

Funkcja βK zadana jest wówczas następującą relacją:

∂β_K

∂r = 2H_Eα

r⁴ψ⁶ . (2.9)

Przyjęto, że funkcje H_E i H_F wyrażają się wyrażeniami uzyskanymi dla czasoprzestrzeni Kerra w pracy [33]. W wyprowadzeniu tym skorzystano z metryki Boyera-Lindquista:

g_µν^BLdx^µdx^ν = −



1 − 2mr_BL Σ_BL



dt²+ Σ_BL

∆_BLdr²_BL+ Σ_BLdθ²+ + ∆_BLΣ_BL+ 2mr_BL(r_BL² + a²)

Σ_BL sin²θ dφ²− 4mr_BLa sin²θ

Σ_BL dtdφ, (2.10) gdzie: Σ_BL = r_BL² + a²cos²θ, ∆_BL = r²_BL− 2mr_BL+ a², m jest masą czarnej dziury, a jej spinem (|a| < m), zaś współrzędna radialna r_BL związana jest z promieniem w metryce (2.1) relacją:

r_BL= r + m + m²− a²

4r , (2.11)

przy czym pozostałe współrzędne są tożsame z tymi w metryce (2.1), tj.: t_BL = t, θ_BL = θ, φ_BL = φ.

Funkcje H_E i H_F wyrażają się wówczas następująco:

H_E = ma [(r_BL² − a²) Σ_BL+ 2r²_BL(r_BL² + a²)]

Σ²_BL ,

H_F = −2ma³r_BL∆^1/2_BL cos θ sin²θ

Σ²_BL .

Dzięki temu przy nieobecności dysku funkcja β_Kjest równa składowej φ wektora przesunięcia w czasoprzestrzeni Kerra, a równania Einsteina redukują się do rozwiązania Kerra we współ- rzędnych quasi-izotropowych. Funkcje metryczne α, β, ψ, q przybierają wówczas następującą postać:

α_K =

s Σ_BL∆_BL

(r²_BL+ a²) Σ_BL+ 2ma²r_BLsin²θ, (2.12a)

β_K = − 2mar_BL

(r²_BL+ a²)Σ_BL+ 2ma²r_BLsin²θ, (2.12b) ψ_K = ⁴

s

(r_BL² + a²)Σ_BL+ 2ma²r_BLsin²θ r²ΣBL

, (2.12c)

q_K = ln Σ_K

p(r²_BL+ a²)Σ_BL+ 2ma²r_BLsin²θ. (2.12d) Po wstawieniu tych funkcji do metryki (2.1) oraz transformacji r → r_BL, odtwarza się z niej metrykę Boyera-Lindquista (2.10).

Zakładając „puncture” w r = 0 oraz ustalając za [33] promień horyzontu czarnej dziury r_h w:

r_h ≡ 1 2

√

m² − a², (2.13)

(przy czym należy zaznaczyć, że w obecności ważkiego dysku wielkości m i a tracą swoje znaczenie fizyczne, a stają się jedynie parametrami związanymi odpowiednio z masą i spinem)

przedefiniowano funkcje metryczne ψ i α, wprowadzając pomocnicze funkcje ϕ i B, zdefiniowane następująco:

ψ ^def=  1 + r_h

r

 e^ϕ, α ^def= ψ⁻¹

1 − r_h r

 e^−ϕB.

Równanie (2.9) można wówczas zapisać następująco:

∂β_K

∂r = 2H_EBe^−8ϕ(r − r_h) r²

(r + r_h)⁷ . (2.14)

Z kolei równania Einsteina przyjmują wtedy postać czterech równań eliptycznych na funkcje q, ϕ, B oraz β_T:



∂_r²+ 2r

r² − r_h²∂_r+ 1

r²∂_θ²+ ctg θ r² ∂_θ



ϕ = S_ϕ, (2.15a)



∂_r²+ 3r²+ r²_h

r (r²− r_h²)∂_r+ 1

r²∂_θ²+2ctg θ r² ∂_θ



B = S_B, (2.15b)



∂_r²+ 4r²− 8rhr + 2r_h²

r (r²− r²_h) ∂_r+ 1

r²∂_θ²+3ctg θ r² ∂_θ



β_T = S_βT, (2.15c)



∂_r²+1

r∂_r+ 1 r²∂_θ²



q = S_q, (2.15d)

w których człony źródłowe S_q, S_ϕ, S_B, S_βT mają postać:

S_ϕ ≡ −2πe^2qψ⁴

 ρhu²_φ

ψ⁴r²sin²θ + ρh(αu^t)²− 2p



− A²

ψ⁸ − ∂_rϕ ∂_rb − (2.16a)

−1

r²∂_θϕ ∂_θb − 1 2

 r − r_h

r (r + rh)∂_rb + ctg θ r² ∂_θb

 ,

S_B ≡ 16πBe^2qψ⁴p, (2.16b)

S_βT ≡ 16πα²e^2qρhu^tu_φ

r²sin²θ − 8∂_rϕ ∂_rβ_T+ ∂_rb ∂_rβ_T− 8∂_θϕ ∂_θβ_T

r² +∂_θb ∂_θβ_T

r² , (2.16c) S_q ≡ −8πψ⁴e^2q

 ρhu²_φ

ψ⁴r²sin²θ − p



+ 3A² ψ⁸ + 2

 r − r_h

r (r + rh)∂_r+ctg θ r² ∂_θ



b + (2.16d) +

 8r_h

r²− r_h² + 4∂_r(b − ϕ)



∂_rϕ + 4

r²∂_θϕ ∂_θ(b − ϕ) , (2.16e) gdzie:

A² ≡

Kˆ_rφ² r²sin²θ +

Kˆ_θφ² r⁴sin²θ, b ≡ ln B.

2.3. Prawo rotacji

W celu dopełnienia układu równań (2.14)–(2.15) należy skorzystać z równań hydrodyna- micznych. W tym celu konieczne jest zadanie prawa rotacji j (Ω), gdzie j jest gęstością krętu.

Funkcja ta opisujące ruch cząstek w dysku, będąc składnikiem relatywistycznych równań Eu- lera ∇_µT^µν = 0. Wiadomo, że są one całkowalne pod warunkiem j = j (Ω) [18, 19]. Wewnątrz

dysku równania Eulera można scałkować na dwa sposoby, uzyskując dwie alternatywne definicje gęstości krętu j:

ˆ

dΩ u^tu_φ+ ln h

u^t = C₁, (2.17a)

ˆ

dΩ hu_φ+ h

u^t = C₂, (2.17b)

gdzie C₁, C₂ = const są stałymi całkowania równania Bernoulliego.

W niniejszej pracy posłużono się równaniem (2.17a), które implikuje następującą definicję gęstości krętu j:

j ≡ u^tu_φ. Można pokazać, że:

j = V²

(Ω + β) (1 − V²), (2.18)

gdzie V jest prędkością liniową cząstek dysku:

V ≡ ψ²

α (Ω + β) r sin θ.

Jedną z najlepiej poświadczonych obserwacyjnie krzywych rotacji jest rotacja keplerowska, jednak w układach złożonych z rotującego wokół czarnej dziury dysku jej ogólnorelatywistycz- nego uogólnienia dokonano dopiero w 2015 r. dla czarnej dziury bez spinu [26]. Z kolei w przypadku niezerowego spinu nastąpiło to w 2017 r., co przedstawiono w pracach [29, 28].

W pozycjach tych przedstawiono numeryczne modele dysków z powyższymi prawami rotacji.

Prawo obejmujące przypadki wirującej czarnej dziury ma następującą postać:

j(Ω) ≡ u^tu_φ= −1 2

d dΩlnn

1 −h

a²Ω²+ 3w^4/3Ω^2/3(1 − aΩ)^4/3io

, (2.19)

gdzie w jest parametrem związanym z masą układu, który dla dysku testowego spełnia równość w =√

m.

Przy założeniu prawa rotacji (2.19) równanie Bernoulliego (2.17a) przyjmuje postać:

hα√

1 − V²n 1 −h

a²Ω² + 3w^4/3Ω^2/3(1 − aΩ)^4/3io−1/2

= C, (2.20)

gdzie C ≡ e^C1.

Przy zadanym równaniu stanu (2.4) oraz prawie rotacji (2.19) równania (2.18) oraz (2.20) w połączeniu z równaniami Einsteina (2.14)–(2.15) pozwalają na pełny opis rozważanego w niniejszej rozprawie układu.

2.3.1. Wyprowadzenie prawa rotacji dla testowego dysku pyłowego

Postulowane prawo rotacji (2.19) bazuje na krzywej rotacji pyłowych dysków testowych leżących na płaszczyźnie równikowej w czasoprzestrzeni Kerra, której postać wyprowadzono analitycznie na dwa sposoby [28, 29].

2.3.1.1. Sposób I

W pierwszym wyprowadzeniu skorzystano z metryki (2.1). Przyjęto bez utraty ogólności m = 1. Prędkość kątowa dana jest wówczas wzorem:

Ω(r) ≡ u^φ

u^t = 8r^3/2

(2r + 1)²− a²^3/2

+ 8ar^3/2

. (2.21)

Gęstość krętu j w funkcji promienia r wyraża się zaś następująco:

j = −N

D, (2.22)

gdzie:

N = n

ar^9/2(2r + 1)²− a²^3/2

+ 8a²r⁶o

a⁴+ a² 8r²− 8r − 2
−

−16ar^3/2 q

(2r + 1)²− a²+ (2r + 1)⁴



·

·



a⁴− a² 8r²− 4r + 2
− 16ar^3/2q

(2r + 1)²− a²+ (2r + 1)² 4r²− 8r + 1

 , D = 2r⁵a²− (2r + 1)²²

a⁷− 3a⁵(1 − 2r)²+ a³8r² 6r²+ 20r + 3
− 24r + 3 −

−a (2r + 1)² 4r²− 8r + 1
2o .

Równania (2.21) i (2.22) zadają funkcję j(Ω) w postaci parametrycznej. W celu wyrugowania zmiennej r zauważono, że wzór (2.21) można zapisać następująco:

ξ² = (2r + 1)²− a²

4r , (2.23)

gdzie:

ξ ≡ Ω⁻¹− a
1/3

. (2.24)

Rozwiązanie równania (2.23) względem zmiennej r i wstawienie do wzoru (2.22) pozwala uzyskać relację j(Ω). Niestety, formuła ta zawiera zagnieżdżone pierwiastki i składa się z bar- dzo dużej liczby członów (rzędu 10³). Uproszczenie takiego wyrażenia jest problematyczne dla niewymiernych wartości a [34] i komputerowe systemy algebry nie są w stanie tego dokonać.

Z drugiej jednak strony, wyrażenie to ulega znacznemu uproszczeniu dla wymiernych wartości a, np. a = −2/3, a obecność wielu powtarzających się pierwiastków w liczniku i mianowniku wskazuje na istnienie znacznie prostszej formuły. Poniżej przedstawiono procedurę uproszczenia wyrażenia j(Ω). Niech:

X ≡ p

a²− 2ξ²+ ξ⁴,

Z ≡ p

1 − a²+ 2X(1 + X − ξ²).

Równanie (2.22) można zapisać jako:

j(ξ) = (a + ξ³) N₁N₂

ξ³D₁D₂D₃ , (2.25)

gdzie:

N₁ ≡ ξ ξ²− 3
Z²+ 2a 1 + X − ξ²
Z, (2.26)

N₂ ≡ ξ⁴+ a²
Z²+ 2aξ 1 + X − ξ²
Z, (2.27) D₁ ≡ ξ² ξ²− 3
2

− 4a², (2.28)

D₂ ≡ ξ²− 1 − X, (2.29)

D₃ ≡ 4ξ⁶− 1 − 3X + (9 + 8X) ξ²− 4 (3 + X) ξ⁴− a² 3 + X − 3ξ²
. (2.30) Można łatwo sprawdzić, że D₂D₃ = Z⁴, skąd:

j(ξ) = (a + ξ³) N₁N₂

ξ³Z⁴D₁ . (2.31)

Z kolei równanie (2.26) można zapisać następująco:

2a 1 + X − ξ²
Z = N₁− ξ ξ²− 3
Z²,

z którego poprzez obustronne podniesienie do kwadratu eliminuje się pierwiastek w Z. Uzyskaną równość można rozwiązać względem N1, uzyskując następujący fizyczny pierwiastek:

N₁ =ξ ξ²− 3
− 2a Z². (2.32)

Postępując analogicznie z równaniem (2.27), otrzymuje się fizyczny pierwiastek N₂:

N₂ = ξ⁴− 2aξ + a²
Z². (2.33) Z drugiej strony:

D1 =ξ ξ² − 3
− 2a ξ ξ²− 3
+ 2a . (2.34) Po wstawieniu równań (2.32), (2.33) i (2.34) do wzoru (2.31) uzyskuje się:

j(ξ) = (a + ξ³) (ξ⁴− 2aξ + a²)

ξ³[ξ (ξ²− 3) + 2a] , (2.35)

a po skorzystaniu z relacji (2.24) i ponownym wprowadzeniu parametru m można dojść do ostatecznego wyniku:

j(Ω) = −1 2

d dΩln

n 1 −

h

a²Ω²+ 3m^2/3Ω^2/3(1 − aΩ)^4/3 io

. (2.36)

2.3.1.2. Sposób II

Metoda ta, choć odkryta poźniej, jest znacznie prostsza. Korzysta się w niej z metryki we współrzędnych Boyera-Lindquista (2.10). Prędkość kątowa zadana jest wówczas relacją [35]:

Ω ≡ u^φ u^t =

√m r_BL^3/2+ a√

m, (2.37)

z której uzyskuje się:

r_BL = m^1/3(1 − aΩ)^2/3

Ω^2/3 . (2.38)

Gęstość krętu j w czasoprzestrzeni Kerra dana jest wzorem:

j ≡ u^tuφ= − g_tφ+ g_φφΩ g_tt+ 2g_tφΩ + g_φφΩ². Korzystając z metryki (2.10), uzyskuje się stąd dla θ = π/2:

j(Ω) = 2ma − r_BL(r²_BL+ a²) Ω − 2ma²Ω 2m (aΩ − 1)²+ r_BL[(r²_BL+ a²) Ω²− 1], skąd po skorzystaniu ze wzoru (2.38) otrzymuje się ostatecznie:

j(Ω) = −1 2

d dΩlnn

1 −h

a²Ω²+ 3m^2/3Ω^2/3(1 − aΩ)^4/3io .

2.3.2. Prawo rotacji dla ważkiego, samograwitującego dysku

Prawo rotacji (2.19) zostało wyprowadzone dla testowego dysku pyłowego. Za pracami [28, 29] założono, że obowiązuje ono także dla ważkich torusów gazowych. W takim przy- padku parametr w przestaje oznaczać pierwiastek z masy czarnej dziury, a staje się nieza- leżnym parametrem wyliczanym dla każdego rozwiązania. Również wielkości m i a nie mogą być identyfikowane odpowiednio z masą i spinem czarnej dziury. Są one wówczas swobodnymi parametrami, przy czym |a| ≤ m, co wynika z konstrukcji numerycznej Shibaty. Należy przy tym zaznaczyć, że parametry w, m i a nie są wielkościami fizycznymi i nie są one wyznaczalne z obserwacji. Prawo (2.19) w granicy newtonowskiej odtwarza prawo Keplera (układy newtonow- skie omówiono w rozdziale 3), co, jeśli uwzględnić fakt, iż spin czarnej dziury nie daje wkładu do rzędów 0PN i 1PN rozwinięcia postnewtonowskiego, wykazano w pracy [26].

2.4. Symetrie równań

Poniżej przedstawione zostaną dwie symetrie charakteryzujące układ równań (2.14), (2.15), (2.18), (2.20). Jedna z nich — symetria skalowania — była wykorzystywana w niniejszej roz- prawie.

2.4.1. Symetria inwersji

Jeśli rozważyć transformację inwersji współrzędnej r względem powierzchni r = r_h:

¯ r ≡ r_h²

r ,

to rozważany układ równań pozostaje względem niej niezmienniczy [21]. W związku z tym formę przetransformowaną można uzyskać, przemianowując zmienną r na ¯r oraz korzystając z następujących podstawień:

ϕ (r) = ϕ (¯r) , B (r) = B (¯r) , β_K(r) = β_K(¯r) , β_T(r) = β_T(¯r) , q (r) = q (¯r) , ψ (r) = r_h

¯ r ψ (¯r) , r_BL(r) = r_BL(¯r) HE(r) = HE(¯r) , H_F(r) = −H_F(¯r) ,

A²(r) = rh

¯ r

12

A²(¯r) , α (r) = −α (¯r) ,

ρ (r) = ρ (¯r) , h (r) = h (¯r) , p (r) = p (¯r) , u_µ(r) = −u_µ(¯r) , u^µ(r) = −u^µ(¯r) , V (r) = −V (¯r) ,

przy czym znak w transformacji 4-prędkości oraz prędkości liniowej wynika z transformacji funkcji α.

Relacje z funkcjami φ, B, β_T i q Shibata wykorzystuje do sformułowania warunków brzego- wych na horyzoncie r = rhprzy rozwiązywaniu układu równań (2.15), a które zostaną przedsta- wione w rozdziale 2.5.2. Pozostałe zostały przeze mnie wyprowadzone z intencją przeprowadze- nia kompaktyfikacji przestrzeni, co mogłoby być użyteczne przy numerycznym rozwiązywaniu układu równań (2.14), (2.15), (2.18), (2.20).

2.4.2. Symetria skalowania

Równania opisujące badany układ są również niezmiennicze względem następującej trans- formacji skalowania:

r = λ˜r m = λ ˜m

a = λ˜a ϕ = ϕ,˜ B = B,˜ β_K = λ⁻¹β˜_K, β_T = λ⁻¹β˜_T,

q = q,˜ ψ = ψ,˜ A² = λ⁻²A²,

α = α,˜ ρ = λ⁻²ρ,˜ h = ˜h, p = λ⁻²p,˜ u^t = u˜^t, u_t = u˜_t, u^φ = λ⁻¹u˜^φ, u_φ = λ˜u_φ,

V = V ,˜

gdzie λ ma wymiar długości, zaś zmienne po przeskalowaniu są bezwymiarowe.

Symetria ta została wykorzystana w rozdziałach 5, 6 i 7, a szczegóły jej aplikacji przedsta- wiono w rozdziale 2.5.6.

2.5. Opis metody numerycznej

Wykorzystana w obliczeniach metoda numeryczna bazuje na pracy [21] z kilkoma modyfi- kacjami, z których najważniejsza wynika z zastosowania innego prawa rotacji, które wymaga rozwiązania równania (2.18) w celu wyznaczenia prędkości kątowej Ω. Kluczowy wkład w jej implementację wniósł dr hab. Patryk Mach. W metodzie tej kluczowe znaczenie ma rozwią- zanie Kerra, stanowiące dane początkowe dla funkcji metrycznych, a odzyskiwane w granicy nieważkiego dysku. Warto podkreślić, że użyte w niniejszej rozprawie prawo rotacji (2.19) po- zwala na uzyskanie rozwiązań z dowolnie lekkim dyskiem, co nie jest możliwe np. w przypadku rotacji sztywnej [20]. Rozwiązanie Kerra było również używane w celu testowania poprawności

działania kodu numerycznego. Z tego też powodu warunki brzegowe zostały dobrane tak, aby jak najlepiej odtworzyć to rozwiązanie przy nieobecności dysku.

2.5.1. Siatka

Użyta siatka numeryczna ma rozmiar N_r× N_θ, gdzie N_r i N_θ są liczbami węzłów względem zmiennych r i θ, oraz posiada następujące brzegi:

— rin ≡ rh — brzeg wewnętrzny znajduje się na horyzoncie czarnej dziury,

— r_ex — brzeg zewnętrzny siatki znajduje się bardzo daleko od zewnętrznego brzegu dysku, aczkolwiek pozostaje skończony,

— θ = 0 — oś symetrii układu,

— θ = π/2 — płaszczyzna równikowa, będąca z założenia płaszczyzną symetrii.

W kierunku radialnym punkty siatki zostały rozmieszczone geometrycznie:

r_i = r_infⁱ⁻¹, i = 1, . . . , N_r, (2.39) gdzie f jest stałą. Z kolei dla zmiennej kątowej posłużono się następującą definicją:

θ_j =







0, j = 1

arc cos1 + ³₂ − j
∆µ , j = 1, . . . , N_θ− 1

π

2, j = N_θ

,

gdzie ∆µ ≡ π/ [2 (N_θ− 1)].

O ile nie zaznaczono inaczej, przedstawione w niniejszej rozprawie wyniki uzyskano na siatce o rozmiarach: Nr = 802, N_θ = 102.

2.5.2. Warunki brzegowe

Z założenia symetrii względem płaszczyzny równikowej wynikają następujące warunki:

∂_θq|_θ=π/2= ∂_θϕ|_θ=π/2= ∂_θB|_θ=π/2= ∂_θβ_T|_θ=π/2= 0. (2.40) Na osi założono regularne warunki brzegowe:

q(θ = 0) = ∂_θϕ|_θ=0= ∂_θB|_θ=0= ∂_θβ_T|_θ=0= 0, (2.41) gdzie znikanie funkcji q na osi wynika z lokalnej płaskości metryki.

Z formalizmu „puncture” wynika istnienie symetrii inwersji na powierzchni r = r_h. Wynikają stąd następujące warunki brzegowe:

∂_rq|_r=r

h = ∂_rϕ|_r=r

h = ∂_rB|_r=r

h = ∂_rβ_T|_r=r

h = 0. (2.42)

W przypadku funkcji β_T konieczne jest zadanie dodatkowych warunków, za które przyjęto [21]:

βT(r = rh) = ∂_r²βT

r=rh = ∂_r³βT

r=rh = 0. (2.43)

Wybór ten jest do pewnego stopnia arbitralny, a związany jest ze swobodą podziału funkcji β na część pochodzącą od czarnej dziury β_K oraz dysku β_T. Ma on wpływ na definicję oraz własności krętu czarnej dziury.

Warunki brzegowe dla r = rex wynikają z rozwinięcia multipolowego funkcji metrycznych przy r → ∞ oraz z warunku asymptotycznej płaskości czasoprzestrzeni:

ϕ ^r→∞−→ M₁

2r , (2.44a)

B ^r→∞−→ 1 − B₁

r², (2.44b)

β_T ^r→∞−→ −2J₁

r³ , (2.44c)

q ^r→∞−→ q₁sin²θ

r² , (2.44d)

gdzie:

M1 = −2 ˆ∞

rin

dr r²− r²_in
ˆπ/2

0

dθ sin θ Sϕ,

B1 = 2 π

ˆ∞

rin

dr(r²− r²_in)² r

ˆπ/2

0

dθ sin²θ SB,

J₁ = 4π ˆ∞

rin

dr r² ˆπ/2

0

dθ sin θ ραu^tψ⁶e^2qhu_φ,

q₁ = 2 π

ˆ∞

rin

dr r³ ˆπ/2

0

dθ cos(2θ) S_q− 4 πr_in²

ˆπ/2

0

dθ cos(2θ) q(r_in, θ) .

Ponieważ promień r_ex zewnętrznego brzegu siatki jest skończony, w kodzie numerycznym powyższe całki liczono w granicach od r_in do r_ex. Z tego również powodu założono przy cał- kowaniu równania (2.14), że β_K(r = r_ex) = 0, co jest przybliżeniem warunku β_K → 0 przy r → ∞.

2.5.3. Dyskretyzacja równań

W kodzie numerycznym równania Einsteina (2.15) zostały zdyskretyzowane z wykorzysta- niem wzorów trójpunktowych na pochodne. W przypadku warunków brzegowych sprawa jest o tyle delikatna, że użycie różnych formuł dyskretyzacji pochodnych może prowadzić do lepszego lub gorszego odtworzenia czasoprzestrzeni Kerra przy braku dysku. Zdecydowano się na wybór kombinacji, która prowadzi do najlepszej zgodności z rozwiązaniem Kerra. Użyto w niej wzorów trójpunktowych, z wyjątkiem warunków brzegowych na horyzoncie r = rh na funkcje ϕ oraz β_T:

— w równaniu ∂_rϕ|_r=r

h = 0 użyto wzoru czteropunktowego,

— w równaniach ∂rβ_T|_r=r

h = ∂_r²β_T|_r=r

h = ∂_r³β_T|_r=r

h = 0 użyto wzorów sześciopunktowych.

Wszystkie całki liczone były z użyciem reguły trapezu.

2.5.4. Dane początkowe

W omawianej procedurze zadawane były następujące dane początkowe:

— parametr spinu czarnej dziury a,

— parametr masy czarnej dziury m,

— wewnętrzny i zewnętrzny promień dysku: r₁ i r₂ (a dokładniej, zadaje się wartości wyzna- czające pierwsze następujące po nich węzły siatki),

— wykładnik adiabatyczny w równaniu stanu γ,

— maksymalną gęstość barionową materii w dysku ρ_max,

— zgadnięcie początkowe na wartości prędkości kątowej na wewnętrznym i zewnętrznym brzegu dysku na równiku: Ω1 i Ω2.

Jak już wspomniano, danymi początkowymi dla funkcji metrycznych było rozwiązanie Kerra.

Można było w tym celu również wykorzystać uzyskane już rozwiazanie, stosując procedurę tzw.

restartu, dzięki czemu można było uzyskiwać rozwiązania dla szerszej klasy wartości parame- trów wejściowych.

Przy rozpoczynaniu obliczeń od rozwiązania Kerra zgadnięcie początkowe parametrów Ω₁ i Ω2 należało dobrać tak, aby były one bliskie wartościom wynikającym z równania (2.21), choć, z uwagi na obraną metodę ich wyliczania, nieznacznie od nich mniejsze.

Z przyczyn technicznych konieczne było zadanie jeszcze dwóch parametrów, które zostaną teraz przedstawione.

Pierwszy z nich związany jest ściśle z prawem rotacji (2.19). Zaobserwowano, że w czasoprze- strzeni Kerra dysk z ustalonymi brzegami spełniający prawo rotacji (2.19) nie daje rozwiązania o rozsądnej grubości. Generowana jest wówczas pusta czasoprzestrzeń Kerra, bez możliwości uzyskania ważkich torusów. Z tego powodu w pierwszych iteracjach należy odstąpić od rotacji keplerowskiej. W tym celu w prawie rotacji (2.19) wprowadzono osobny parametr spinu a_rot. Jego wartość przy rozpoczynaniu obliczeń od czasoprzestrzeni Kerra w pierwszych iteracjach była różna od parametru a (w obliczeniach przyjęto początkową wartość a_rot = 0; jeżeli nato- miast docelowo a = 0, wówczas pierwsze iteracje liczono z a 6= 0), a następnie, po skorzystaniu z procedury restartu, zostawała z nim zrównana.

Drugi parametr związany jest z wyznaczeniem prędkości kątowej Ω. Wykorzystywano do tego równanie (2.18), które jednak w pobliżu osi θ = 0 ma rozbieżne rozwiązania. Z tego powodu wprowadzono tzw. kąt ucięcia θsec, który określa minimalną wartość kąta θ, przy któ- rym wyliczana jest prędkość kątowa Ω. Począwszy od drugiej iteracji, obszar ten może zostać dopasowany do kształtu dysku, aby zminimalizować ryzyko rozbieżności.

Należy także zaznaczyć, iż równania (2.18) nie można było rozwiązać dla Ω < 0, w związku z czym w przypadku dysków kontrrotujących względem czarnej dziury konieczne było przyjęcie konwencji, w której a < 0, a Ω > 0. Jednak, rozpoczynając obliczenia od czasoprzestrzeni Kerra, nie można było od razu użyć ujemnej wartości parametru a. Zatem w pierwszych ite- racjach przyjmowano a > 0, a następnie, po skorzystaniu z procedury restartu, dokonywano podstawienia a → −a.

2.5.5. Procedura numeryczna

Na początku obliczeń funkcje metryczne wyliczano według wzorów dla czasoprzestrzeni Kerra (2.12), bądź też wczytywano je z uzyskanego już rozwiązania, wykorzystując procedurę restartu.

Każdy krok iteracyjny rozpoczynał się od wyliczenia warunków brzegowych na równiku (2.40), osi (2.41) i horyzoncie (2.42)–(2.43) dla funkcji metrycznych ϕ, B, β_T, q.

Następnie z równań (2.18) oraz (2.20) w punktach (r₁, π/2) oraz (r₂, π/2), w których z definicji h = 1, wyznaczane były wartości Ω₁, Ω₂, w, C. Jest to układ czterech równań na cztery niewiadome, z których jedynie dwa równania są niezależne. Poprzez odpowiednie pod- stawienie uzyskano układ dwóch równań na Ω₁ i Ω₂, który można prosto rozwiązać metodą Newtona-Raphsona. Korzystając z uzyskanych rozwiązań, wyznaczano następnie stałe w oraz C.

Kolejnym rachunkiem było wyliczenie z równania (2.18) prędkości kątowej Ω. Obliczenia tego dokonywano w ograniczonym obszarze, zadanym przez następujące nierówności: r sin θ ≥ r₁, r ≤ r₂, θ_sec ≤ θ ≤ π/2, gdzie kąt ucięcia θ_sec w pierwszej iteracji był stałą, w kolejnych zaś mógł być (i w obliczeniach do niniejszej rozprawy zawsze był) na bieżąco dopasowywany do kształtu