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₁

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

dysku jako funkcja promienia ri tak, aby dla każdego ri kąt θsec(ri) znajdował się w pierwszym węźle, w którym h = 1, licząc od równika. Dopasowania tego dokonywano pod koniec każdej iteracji.

Kolejnym krokiem było wyznaczenie entalpii właściwej h, korzystając z równania (2.20).

Następnie wyliczano stałą proporcjonalności K w równaniu stanu, korzystając z parametru ρ_max oraz maksymalnej wartości entalpii właściwej h, która wyszukiwana była w całej obję-tości dysku. Jest to ulepszenie sposobu Hachisu [36], który przeszukiwał jedynie płaszczyznę równikową. Dr Michał Piróg zauważył, że niniejsza modyfikacja znacznie poprawia zbieżność w metodzie punktu stałego. Wynika to z faktu, że w trakcie iteracji maksimum gęstości prze-suwa się tymczasowo poza płaszczyznę równikową, po czym stopniowo na nią wraca. Następnie wyliczano gęstość barionową materii w dysku ρ, korzystając ze wzoru:

ρ = γ − 1

Kγ (h − 1)

_γ−1¹ .

Kolejnym rachunkiem było wyznaczenie źródeł S_ϕ, S_B, S_βT oraz S_q, korzystając z równań (2.16), a także współczynników M₁, B₁, J₁ oraz q₁ w warunkach brzegowych (2.44) dla r = r_ex. Następnie rozwiązywano układ równań (2.15) metodą punktu stałego, uzyskując nowe war-tości funkcji ϕ, B, β_T, q, wykorzystywane w kolejnym kroku iteracyjnym. Wyliczano także funkcję β_K, całkując równanie (2.14).

Na koniec każdej iteracji dokonywano wspomnianego już dopasowania obszaru wyliczania prędkości kątowej Ω do kształtu dysku, a w części obliczeń dokonywano także modyfikacji parametru ρ_max celem uzyskania dysku o odpowiedniej masie lub kręcie, a także wartości r₁ i r₂ w celu otrzymania dysku o zadanych wartościach promienia obwodowego r_c na brzegach, odpowiednio r_c(r₁) oraz r_c(r₂).

2.5.6. Skalowanie

W przedstawionych w rozdziałach 5, 6 i 7 obliczeniach ogólnorelatywistycznych posłużono się zmiennymi bezwymiarowymi uzyskanymi po dokonaniu skalowania opisanego w rozdziale 2.4.2. W przeciwieństwie do omówionego w rozdziale 3 przypadku newtonowskiego z punktową masą centralną, w obecności czarnej dziury niemożliwe jest dokonanie przeskalowania, które uniezależniłoby rozwiązania zarówno od promienia r2, jak i gęstości maksymalnej ρmax. W niniejszej rozprawie we wszystkich obliczeniach przyjęto λ = r₂, wskutek czego ˜r₂ = 1.

Celem zagwarantowania obecności punktu ˜r = 1 w siatce numerycznej przyjęto następującą wartość stałej f we wzorze (2.39):

f ≡ ˜r⁻

1 i2−1

in ,

skąd ˜r₂ ≡ ˜r (i₂) = 1.Wartość i₂ dobrano tak, aby ˜r_ex ≈ 20.

Przeprowadzone w ten sposób skalowanie ułatwia porównanie wyników przedstawionych w rozdziale 5 z odpowiadającymi im rezultatami pochodzącymi z prac [1, 2], a zaprezentowanymi w rozdziale 3.