The impact of radiation on the dynamics and appearance of the supermassive black hole in the Galactic center

(1)

University of Amsterdam

MSc Physics and Astronomy

Track: Astronomy & Astrophysics

Master Thesis

The impact of radiation on the dynamics and

appearance of the supermassive black hole in

the Galactic center

“Have you seen Interstellar?”

by

David Jacob van Eijnatten

10645462

July 29, 2018

60 ECTS

2017 - 2018

Supervisor:

prof. dr. S.B. Markoff

Daily supervisor:

Koushik Chatterjee, MSc

Examiners:

prof. dr. S.B. Markoff

prof. dr. M.B.M. van der Klis

(2)

Scientific abstract

Background: Sagittarius A* (Sgr A*) is the supermassive black hole of 4.6×106Min the

center of the Milky Way. It is accreting at a very low rate of 10−9− 10−7_M

/yr. Flares

are observed in the near-infrared and X-ray frequencies. The origin of the flares is the accretion flow near the black hole, but its physical mechanism is unknown. Recently, the Event Horizon Telescope imaged the flow around the black hole with a resolution that is on the order of the event horizon of the black hole. This work focuses on discovering the origin of these flares and the appearance of Sagittarius A* on these event horizon scales through radiative cooling.

Method: We use the general relativistic magnetohydrodynamic simulation code HAMRR in 2D combined with semi-analytical cooling functions for optically thin flows. The electrons are assumed thermal and the temperature is some fraction < 1 of the proton temperature. We combine this with the ray-tracing code BHOSS to generate images and spectra. The gas starts out in a disk with a radius of 50 RG (RG = GM_c2 ) around the black hole and is

allowed to be accreted for a minimum of one day.

Results: We find that even an accretion disk around a black hole with an accretion rate as low as Sgr A* is dynamically heavily affected by radiative cooling. The temperature can drop an order of magnitude compared to a non-cooled disk. Its spectral appearance is also affected. The synchrotron self-absorption frequency shifts, the flux decreases and a large portion of the high energy tail disappears from the synchrotron spectrum when cooling is incorporated. We find a best-fit accretion rate of ˙M = 1 × 10−8M/yr at an

inclination of 85 deg and a temperature ratio of 0.3. Unfortunately we were not able to statistically investigate flares as the MRI is not resolved anymore after a few days of simulation time.

Conclusion: At all accretion rates the dynamics of an accretion disk are affected by the loss of energy through radiation. When analyzing the results from the Event Horizon Telescope it is important to account for radiative losses. The accretion rate we find is higher than was found previously with other studies, although we are consistent when we lower our resolution. Flaring studies using GRMHD simulations are difficult without high-resolution 3D simulation with strongly magnetized disk where the MRI functions for tens of days.

(3)

Popular abstract

Black holes are some of the most fascinating objects in the cosmos. They can grow to masses millions or billions times heavier than the sun, they can impede the growth of galaxies and sling out matter to distances trillion times their diameter. Our Milky Way also contains one of these supermassive black holes, called Sagittarius A*. Its mass is about five million times the mass of the sun and we can observe it sucking up (“accreting”) small amounts of gas that is swirling around in its neighbourhood. We do not fully understand what we see, however: almost every day we see “flares”, burst of radiation that last a short time. To investigate this and to observe the black hole “shadow”, the shape the black hole imprints on the surrounding gas, astronomers are taking a picture of Sagittarius A* with a resolution that is about the diameter of the black hole with the Event Horizon Telescope. To understand what we will see in this picture and to explain the flaring behaviour of Sagittarius A* we use computer simulations to replicate the interaction between the black hole and its surrounding gas. This gas behaves like a fluid containing a magnetic field, so the simulated physics is the laws of fluid dynamics and general relativity. To this, I added optically thin radiation. Radiation allows energy to escape in the form of light. We find that this energy loss has a significant impact on how the gas behaves and on what the picture from the Event Horizon Telescope will look. We also found that that although our simulations were useful to determine the appearance of the accretion process, it was not possible to simulate this process long enough to gather enough statistics to say anything about the flares.

(4)

Samenvatting

Zwarte gaten behoren tot de meest fascinerende objecten in het heelal. Ze kunnen miljoe-nen tot miljarden keer zo zwaar zijn als de zon, ze bepalen de formatie en groei van sterrenstelsels en ze kunnen enorme hoeveelheden materiaal over lengten afschieten die biljarden mall zo lang zijn als hun diameter. Ons melkwegstelsel heeft ook een superzwaar zwart gat, genaamd Sagittarius A*. Het weegt 4,6 miljoen keer zo veel als de zon en het is geobserveerd dat het kleine hoeveelheden gas opzuigt (“accreteert”). Hoewel we dit weten, begrijpen we niet alles wat we zien: elke dag zien we flikkeringen van korte duur in het geobserveerde licht. Om deze flikkeringen te onderzoeken en om de “schaduw” van het zwarte gat te zien, nemen astronomen een foto van het zwarte gat met een resolutie hoger dan de diameter van het zwarte gat met de Event Horizon Telescope. De “schaduw” is de vertekening van het gas door het zwarte gat. Om deze foto en de flikkeringen te begrijpen draaien wij computer simulaties van de de relevante natuurkunde rond het zwarte gat. Het gas gedraagt zich als een gemagnetiseerde vloeistof en aldus simuleren wij de vloeistofwetten en zwaartekracht. We hebben diverse stralingsprocessen toegevoegd, waardoor energie kan ontsnappen aan het systeem in de vorm van licht. We hebben ondervonden dat dit energieverlies een significant effect heeft op het gedrag van het gas en op het uiterlijk van accretieproces dat wij op aarde kunnen observeren. Hoewel onze simulaties nuttig waren om dit uiterlijk te bepalen, was het lastig de origine van de flikkeringen te bepalen omdat onze simulaties het verzamelen van voldoende statistieken niet toelieten.

(5)

List of Figures

1.1 The M-σ relation. Black hole mass is correlated with the bulge velocity dispersion in galaxies through a powerlaw. This implies that black hole growth and galaxy formation are related. Taken from Gillessen et al. (2009a). 7

1.2 Image of galaxy cluster MS 0735.6+7421. In orange, optical data from the Hubble Space Telescope is shown, indicating the location of the stars. In purple, radio data from the Very Large Array is shown, indicating the location of the jet. In blue, X-ray data from the Chandra X-ray observatory is shown. The cavities in the in the X-ray emitting gas coincide with the jet, implying that the jet “blew bubbles” in the surrounding medium. This evacuation of gas could halt galaxy growth and is called AGN feedback. Credit: Chandra archive. . . 8

1.3 Combined radio image from the Very Large Array and Green Bank Tele-scope. In the center of the image the Galactic center, Sagittarius A (Sgr A) is visible. Some supernova remnants (SNRs) and non-thermal radio filaments (NRFs) are also visible. Figure taken from Yusef-Zadeh et al. (2004). . . 9

1.4 The orbit of S2 around Sgr A* as monitored by the Very Large Telescope (VLT) And Keck Observatory. S2 is on a very elliptical orbit (e = 0.88) with a pericentre of 120 AU. The best-fit mass of the central body is 3.7 ± 0.15 × 106_M

Credit: The MPE-IR Galactic Center Team. . . 10

1.5 Ray-traced images of a simulation of Sgr A*. Left: 1.3 mm image at an angle of 60 deg accreting at 4.2 × 10−8M/yr. Right: scatter-broadened

version of the left image. The color scale is normalized intensity. In both images the black hole shadow is easily visible as the dark area surrounded by radiation. Credit: Mo´scibrodzka et al. (2014) . . . 11

1.6 The three canonical accretion disk models. On the left the scale height H/R and optical depth τ are given. In the middle a schematic drawing is given of the system, where the arrows refer to an outflow. On the right the luminosity is given as a fraction of the Eddington luminosity, the theoretical luminosity at which no more matter can be accreted. Credit: Alexander Tchekhovskoy. . . 12

1.7 A typical multiwavelength spectra of Sgr A*. The dot-dashed line is the thermal synchrotron and SSC emission, the short-dashed lines in the non-thermal synchrotron emission, the dotted line is the total synchrotron and SSC emission, the long-dashed line is the bremsstrahlung emission and the solid line is the total quiescent emission. The upper bowtie is a strong X-ray flare. The data was not taken simultaneous but over a period of 5 years. Taken from Yuan et al. (2003). . . 14

(7)

1.8 Three hardness-intensity diagrams for different compact objects from Fender (2010). Hardness is defined as the ratio of fluxes in two X-ray bands, where high hardness (A in the diagram) indicates a nonthermal spectrum and low hardness (B in the diagram) indicates a thermal spectrum. The sources travel counterclockwise around the diagram and spend most of their time in the lower right corner, where they are faintest. For neutron stars and black holes hardness is used on the x-axis while for white dwarfs a definition analogous to the DFL is used: powerlaw luminosity over total luminosity. 15

1.9 The DFLD from K¨ording et al. (2006). This diagram is similar to the HID of stellar mass compact objects. Powerlaw luminosity over total luminosity is plotted versus total luminosity. R is the radio loudness, the ratio of radio over optical luminosity. The data is binned and only bins with more than 10 counts are shown. The gap between the LLAGN on the lower right and the quasars on the top is caused by various selection effects in the used survey. As you can see, LLAGN statistically show a similar pattern as XRBs even with the quiescent state in the lower right being the most heavily populated. . . 16

1.10 The fundamental plane of black hole accretion. There is a powerlaw rela-tion between radio and X-ray luminosity for all low-luminosity black hole sources when corrected for mass. The only exception is Sgr A*, even its most powerful flares shown here can not boost its X-ray emission to the to sufficient heights. The light green datapoints show other black hole sources. Taken from Markoff (2005). . . 17

2.1 A schematic of the MRI. Two particles in a differentially rotating flow connected by a spring exercise forces on each other. Angular momentum is transferred from the inner particle to the outer particle. The spring in this case is a magnetic field line. Credit: Nick Murphy. . . 24

4.1 The initial setup of the simulation. The torus is plotted in density, with the black hole in the middle. The black lines in the disk give the poloidal magnetic field loops. The gray lines give contours of the grid. Not all cells are plotted. . . 44

4.2 The initial temperature profile of the simulation. These are temperatures typical for an ADAF, going up to > 1011 K. The black dotted lines give H/R = 1. The inner regions of the ADAF also agree with this. . . 44

4.3 The accretion rate for three different accretion rates with Te/Ti=0.3 with

a resolution of 10242. Their average values are very close to the intended accretion rate. They also remain constant over the whole simulation domain. 45

4.4 Angle-averaged, radial profiles of density, magnetic field strength and tem-perature for three different accretion rates. The solid lines display the cooled simulations M9T3HR, M8T3HR and M7T3HR. The dashed lines display the non-cooled simulation HR scaled to the same density as the three cooled simulations. The density and magnetic field scale with the accretion rates as expected. The profiles are clearly affected by the cooling. We do not expect the temperature to explicitly scale with accretion rate, although the temperature is severely lowered by radiative cooling. . . 46

(8)

4.5 A plot of the angle-dependent density profile for M9T3HR, M8T3HR and M7T3HR. As the accretion rate rises, the flow starts to collapse and the equatorial density rises with two orders of magnitude. . . 47

4.6 A plot of the angle-dependent distribution of the positron to proton ratio z for M9T3HR, M8T3HR and M7T3HR. This ratio times two tells you how much leptons there are due to pair-production compared to the leptons from ionized gas. The collapsed high density, low temperature areas of the disk have a reduced value for z. . . 48

4.7 The angle-dependent distribution of the synchrotron self-absorption quency for M9T3HR, M8T3HR and M7T3HR. The self-absorption fre-quency rises with the accretion rate, making the optically thick inner flow appear increasingly larger. . . 48

4.8 Three different time series. On top, the total spatially integrated cooling rate in a 10 RGsphere around the black hole is given in arbitrary units. In

the middle, the accretion rate over the event horizon is given in arbitrary units. Below, the average ratio of the MRI wavelength over the cell size is given separately for each dimension. The horizontal line is drawn at a ratio of unity where the MRI wavelength drops below the resolution of the simulation. The φ-component of the MRI can never be resolved due to the axisymmetric nature of these simulations.. . . 50

A.1 The 230 GHz images and radio spectra at i = 5 deg, i = 45 deg and i = 85 deg for M9T3HR. Observations taken from Falcke et al. (1998); Sch¨odel et al. (2003); Bower et al. (2015). . . 64

A.4 The 230 GHz image at i = 45 deg and spectra of M9T3HR compared to HR for which we assumed ˙M = 1 × 10−9M/yr. Observations taken from

Falcke et al. (1998); Sch¨odel et al. (2003); Bower et al. (2015). . . 67

Falcke et al. (1998); Sch¨odel et al. (2003); Bower et al. (2015). . . 68

(9)

List of Tables

2.1 The radiative and pair-producing processes in a non-magnetized plasma. Taken from Svensson (1982a). . . 31

4.1 The simulations discussed in this work. We ran 2D simulations at three different accretion rates within the observational constraints of Sgr A*. We run each of these for two different temperature ratios. at a low and high resolution of 5122 and 10242, respectively. We ran one very high resolution purposes for a convergence test. For each resolution we also ran a noncooled test to benchmark the cooling effects against. The naming works as follows: M is followed by the negative logarithm of the accretion rate, T is followed by ten times Te/Ti. This is followed by the resolution.

LR refers to low resolution (5122_{), HR refers to higher resolution (1024}2₎

(10)

CHAPTER

1

Introduction

1.1 Supermassive black holes

When Albert Einstein presented the theory of general relativity in 1916 (Einstein,1916), one of the first solutions found to the Einstein equations was the vacuum solution near a static, neutral and spherical massive object, like a star or planet. However, the solution also seemed to indicate that any object that was sufficiently dense as to have a radius r < rS = GM_c2 , where rS is called the Schwartzschild radius, G is the gravitational

con-stant, M is the mass of the object and c is the speed of light, would collapse to a point or singularity and nothing inside rS could ever leave this sphere, not even light. What is

more, it is possible for some core collapse supernovae from massive stars to achieve this density. Black holes are now observed beyond a reasonable doubt across a huge range of masses: from a couple to ten times the mass of the sun (stellar-mass black holes) as detected in binary systems to millions or billions of solar masses (supermassive black holes or SMBHs) as detected in the cores of some galaxies. The formation paths of the stellar-mass black holes are clear: direct formation from supernovae, collapse of accret-ing neutron stars or binary neutron star mergers with a final core mass exceedaccret-ing the Tolman-Oppenheimer-Volkoff limit, the theoretical maximum mass a neutron star can have (Tolman, 1939; Oppenheimer & Volkoff, 1939). The formation paths of supermas-sive black holes are less well understood (Volonteri, 2010). Most galaxies are thought to to contain their very own SMBH and these galaxies also seem to be heavily influenced by their central black hole even though the black hole contains only a fraction of the total galaxy mass. This influence is most clearly expressed in the form of the M − σ relation, the strong correlation between black hole mass and stellar velocity dispersion in the bulge (Ferrarese & Merritt, 2000; Gebhardt et al., 2000), see Fig. 1.1. The most widely accepted explanation uses accretion powered outflows as the main influence on the neighbouring regions (Silk & Rees, 1998; King, 2003). More recently it was discovered that galaxy merger models could also reproduce the M − σ relation without any interac-tions of the black holes with the surrounding galaxies (Hirschmann et al., 2010).

Large scale cosmological simulations of galaxy formation tend to overpredict the amount of massive galaxies that form (Bower et al., 2006b). The only known mechanism that can quench galaxy growth on such scales are large scale jets produced by SMBHs.

(11)

Figure 1.1: The M-σ relation. Black hole mass is correlated with the bulge velocity dispersion in galaxies through a powerlaw. This implies that black hole growth and galaxy formation are related. Taken from Gillessen et al. (2009a).

Some supermassive black holes are “active”: they are accreting significant amounts of gas and occasionally driving an outflow. These sources are called “active galactic nuclei” (AGN) and produce some of the most impressive images in astronomy, see Fig. 1.2. Clearly these outflows must have a huge influence on their environment by heating and pushing surrounding gas. Other non-active or “quiescent” galaxies show a history of outflow activity (Lintott et al.,2009;Comerford et al.,2017). This suggests that SMBHs can switch on and off, a behaviour more clearly seen in X-ray binaries (XRBs), which undergo this switch in more human time scales of a few to tens of years (Falcke et al.,

2004; Remillard & McClintock, 2006). The influence of SMBH on their environment is usually called AGN feedback. This feedback can have a large impact on star formation rates in the galactic bulge, galaxy formation and galaxy evolution (Fabian,1994;Di Mat-teo et al.,2005;Croton et al.,2006;Fontanot et al.,2006;Hardcastle et al.,2007;Fabian,

2012). By pumping a huge amount of energy in their environment, AGN self-regulate their growth. Both low-density interstellar matter (ISM) is evacuated from the bulge and young stars commonly associated with winds do not form. This is observed as a proportionality between galaxy mass and SMBH mass (Magorrian et al.,1998).

(12)

Figure 1.2: Image of galaxy cluster MS 0735.6+7421. In orange, optical data from the Hubble Space Telescope is shown, indicating the location of the stars. In purple, radio data from the Very Large Array is shown, indicating the location of the jet. In blue, X-ray data from the Chandra X-ray observatory is shown. The cavities in the in the X-ray emitting gas coincide with the jet, implying that the jet “blew bubbles” in the surrounding medium. This evacuation of gas could halt galaxy growth and is called AGN feedback. Credit: Chandra archive.

1.2 Sagittarius A*

The Milky Way also has a supermassive black hole at its center. Because the Milky Way is the only galaxy that we can not look at from a distance, the story starts with the actual determination of its center.

Before the radio astronomy boom following the end of the second world war, the only way to determine the center of the Milky Way was through optical observations of stellar

(13)

motions (Goss & McGee,1996). The center had already appeared in one of Karl Jansky’s radio observations, but was not identified as the Galactic center (Jansky,1933). It could be inferred that our host Galaxy was a spiral galaxy with its center somewhere in the direction of Scorpius and Sagittarius. Because we can only look towards the center of the Milky Way through its disk, optical observations are almost impossible because of a 25 magnitude extinction (Genzel et al.,2010).

Figure 1.3: Combined radio image from the Very Large Array and Green Bank Telescope. In the center of the image the Galactic center, Sagittarius A (Sgr A) is visible. Some supernova remnants (SNRs) and non-thermal radio filaments (NRFs) are also visible. Figure taken from Yusef-Zadeh et al. (2004).

In the fifties radio emission from Sagittarius was quickly identified as the Galactic center and the source was dubbed Sagittarius A (Sgr A) (Piddington & Minnett,1951), see Fig.

1.3. A few years after the idea was proposed that a supermassive black could be present in the center of every galaxy by Lynden-Bell (1969) and Lynden-Bell & Rees (1971) a compact radio source was found inside Sgr A and was identified with our own SMBH. Both the radio source and the SMBH are referred to as Sagittarius A* or Sgr A* (Brown,

1982). In later studies, the size of Sgr A* was cut down to < 1 A.U. (Ghez et al., 2003;

Shen et al., 2005; Bower et al., 2006a; Doeleman et al., 2008). This indicates a density > 1021M/pc3. The only stable system on Galactic time scales with such a high density

(14)

Figure 1.4: The orbit of S2 around Sgr A* as monitored by the Very Large Telescope (VLT) And Keck Observatory. S2 is on a very elliptical orbit (e = 0.88) with a pericentre of 120 AU. The best-fit mass of the central body is 3.7 ± 0.15 × 106 _M

Credit: The

MPE-IR Galactic Center Team.

The most convincing evidence for the presence of a black hole at the center of our Galaxy comes from monitoring stellar orbits at the Galactic center. In Gillessen et al. (2009a,b) the orbits of 28 stars are presented over a period of 21 years. One of those stars, S2, has completed a full orbit since the monitoring began in 1995. From the orbit of S2 a mass for Sgr A* can be derived, 3.7 ± 0.15 × 106_M

, and a minimum size, 120 A.U., since S2 did

not collide with Sgr A*. Using the periastron approach of S2, the gravitational redshift of its light due to Sgr A* could actually be measured (GRAVITY Collaboration et al.,2018).

(15)

Figure 1.5: Ray-traced images of a simulation of Sgr A*. Left: 1.3 mm image at an angle of 60 deg accreting at 4.2 × 10−8 M/yr. Right: scatter-broadened version of the left

image. The color scale is normalized intensity. In both images the black hole shadow is easily visible as the dark area surrounded by radiation. Credit: Mo´scibrodzka et al.

(2014)

1.3 Flares, spectra and accretion flows

We know that Sgr A* has daily flaring activity in the near-infrared (NIR) and X-rays (Baganoff et al.,2001;Genzel et al.,2003;Ghez et al.,2004). Determining the origin and behaviour of these flares is imperative for understanding and interpreting the data that will be gathered by the Event Horizon Project, see Sec. 1.5. Many different explanations have been put forward for the cause of the increased emission. Due to observational constraints and the short timescales the possible flaring mechanisms are actually quite restricted. All possible physical mechanisms proposed below have in common that they occur in the inner tens of RG in either an inflow or outflow. Markoff et al. (2001)

sug-gested jet activity, Yusef-Zadeh et al.(2006a) suggested an expanding plasmon from the Van der Laan model and Zubovas et al. (2012) suggested tidal distruption of asteroids or exoplanets. Some explanations have to do with the accretion flow, from orbiting hot spots (Broderick & Loeb, 2005, 2006) and wave instabilities (Tagger & Melia, 2006) to nonthermal electrons produced by reconnection or turbulence (Yusef-Zadeh et al., 2004;

Dodds-Eden et al., 2009).

The accretion flow around Sgr A* is thought to be advection dominated (ADAF;Narayan et al., 1998; Yuan et al., 2003). This is a low accretion rate version of an accretion disk. Accretion disks are produced when gas surrounding a massive object becomes gravitation-ally bound while retaining some angular momentum. The gas starts to spiral around the massive objects and a disk-like object forms, usually just referred to as an accretion disk. Since the disk is not rigid, the gravitational field causes a differential rotation. This in

(16)

turn causes viscous dissipation of energy and gas slowly moves through the disk towards the central source. The importance of viscosity is usually measured by the Reynolds num-ber Re, the ratio of inertial forces to viscous forces. If Re >> 1, inertial forces dominate and the disk is stationary. If Re << 1, viscous forces dominate and the disk moves in-wards. From purely hydrodynamic considerations we obtain Re ∼ 1014_{, meaning the disk}

is completely stable (Frank et al.,1992). We do however see accreting accretion disks, so something else must cause viscous dissipation. This unknown effect was parametrized by the α-prescription by Shakura & Sunyaev (1973). This is simply the phenomenological formula for the kinmatic viscosity ν = αcsH. Later the magneto-rotational instability

(MRI) was discovered by Balbus & Hawley (1991), see Sect. 2.2.2. The appearance of accretion disks at different accretion rates is wildly different. Like stated before, one of the low accretion rate versions is referred to as an ADAF, a rather puffy disk with scale heights (H/R, where H is the height of the flow at radius R) of unity where the energy released by viscous dissipation can not easily be radiated away because of a decoupling between the hot protons and the radiating electrons and as such these flows tend to be extremely hot (up to 1010− 1011_{K) and underluminous (}_{Narayan & Yi}_,₁₉₉₄_{). The term}

advection dominated refers to the fact that most thermal energy is carried over the event horizon or away in an outflow and thereby “cools” the accretion flow instead of radia-tion. As the accretion rate increases, the higher density allows particles to radiate more efficiently and the disk collapses into a razor-thin disk or α-disk (Shakura & Sunyaev,

1973). These emit blackbody radiation and are therefore much brighter than ADAFs even when normalized with accretion rate. At higher accretion rates still, the escape time of a photon becomes longer than the accretion time scale of the parcel of gas that trapped it. This leads to a radiatively inefficient disk again often referred to as a slim disk (Begelman & Meier, 1982; Abramowicz et al., 1988). For a schematic visualization of these three accretion states, see Fig. 1.6.

Figure 1.6: The three canonical accretion disk models. On the left the scale height H/R and optical depth τ are given. In the middle a schematic drawing is given of the system, where the arrows refer to an outflow. On the right the luminosity is given as a fraction of the Eddington luminosity, the theoretical luminosity at which no more matter can be accreted. Credit: Alexander Tchekhovskoy.

Getting back to the flares, their appearance is a well-studied subject. X-ray flares appear as discrete “events”, with processes different to the steady state. The X-ray flares are strong with an amplification of a few to hundreds of times the quiescent flux and the

(17)

quiescent flux is very constant (Baganoff et al., 2003). The NIR flares are less separate from their background. A quasi-periodic oscillation or simply statistical outliers from Brownian noise could also be an explanation (Meyer et al., 2008; Do et al., 2009). This is contradictory to the fact that every X-ray flare seems associated with an NIR flare, al-though not every NIR flare can be associated with an X-ray flare (Hornstein et al.,2007). Due to the larger spatial resolution of X-ray instruments the thermal bremsstrahlung emitted around the Bondi radius could be obscuring the smaller X-ray flares (Witzel et al., 2012, 2018). NIR flares occur on average 4 times a day, while X-ray flares occur twice a day, with strong X-ray flares occurring daily (Baganoff et al.,2003; Eckart et al.,

2006; Porquet et al., 2008).

In 2012, a campaign was launched with the Chandra X-ray observatory to observe Sgr A* for 3 megaseconds or ∼ 35 days. This campaign was called the X-ray Visionary Project (XVP) (PIs: Baganoff, Markoff, Nowak). During the 3 megaseconds, 39 new flares were observed, effectively doubling the statistics from the decade prior. Because of the large number of flares, a statistical study could be done and many distribution of durations, fluences and luminosities were calculated, see Neilsen et al. (2013). It was also found that all of the variable emission is about 10 % of the total emission and was attributable to the inner ∼ 20 RG (Neilsen et al., 2015). In theory, one could try to replicate these

results using global simulation to find the origins of the flares.

A typical multiwavelength spectrum is given in Fig. 1.7. The spectrum has a strong turnover at 1012 _{Hz called the submillimeter peak. It can be explained by thermal}

syn-chrotron radiation with a self-absorpion frequency of 1012Hz (Yuan et al.,2003). Between 109 _{and 10}10_{Hz there is an excess attributable to synchrotron radiation from nonthermal}

electrons. The radio excess could also be explained by an outflow (Falcke & Biermann,

1996), and the variability profile of the spectrum is strong evidence for this (Bower et al.,

2015). The thermal synchrotron radiation is a necessary presence, because otherwise the NIR upper limits would be violated. The quiescent X-ray emission is produced at large scales out to the Bondi radius, while the flared X-rays have to come from close to the black hole (Wang et al.,2013; Neilsen et al., 2015).

(18)

8 10 12 14 16 18 20

log[ν(Hz)]

30 31 32 33 34 35 36 37

log[

ν

L

ν

(erg s

−1

)]

Sgr A

*

Figure 1.7: A typical multiwavelength spectra of Sgr A*. The dot-dashed line is the ther-mal synchrotron and SSC emission, the short-dashed lines in the nontherther-mal synchrotron emission, the dotted line is the total synchrotron and SSC emission, the long-dashed line is the bremsstrahlung emission and the solid line is the total quiescent emission. The upper bowtie is a strong X-ray flare. The data was not taken simultaneous but over a period of 5 years. Taken from Yuan et al.(2003).

1.4 The fundamental plane of black hole accretion and state transi-tions

Some stellar mass compact objects like X-ray binaries (XRBs) go through different accre-tion phases on time scales of months, years or decades. These phases are mostly clearly seen in the hardness-intensity diagram (HID), see Fig. 1.8. X-ray hardness is a measure of how non-thermal that X-ray spectrum is. High hardness means a spectrum that is dominated by a typical powerlaw from a corona or jet base (Markoff et al.,2005), usually accompanied by radio emission from a jet. Low hardness typically means a soft spectrum or a spectrum dominated by a thermal disk component close to the innermost stable cir-cular orbit (ISCO). These compact objects travel counterclockwise along the hysteresis curve of the HID, first increasing in flux while the hardness stays approximately constant, then travelling towards a soft spectrum where the disk dominates and the jet turns off. As the source gets fainter, it goes towards a quiescent hard spectrum. These sources have a very different formation path (supernovae versus mergers and accretion for XRBs) and accrete very differently (Roche lobe overflow in XRBs versus interstellar gas and tidal disruption events for SMBHs) than SMBHs. But there are some indications that their supermassive counterparts go through state transitions also. Time scales close to the

(19)

black hole should scale directly with the black hole mass (see Chap. 2) so that state transitions can never be observed directly since there is at least a factor of a few million between masses of Galactic and supermassive black holes. However, the HID of SMBHs could be inferred statistically by plotting the current sample. There is another problem: in standard disk theory, the peak temperature of the disc scales as follows: T ∝ ˙M−1/4 (Shakura & Sunyaev,1973). Where the disk shines in the X-rays for Galactic black holes, it shines in the optical or UV for SMBHs. In K¨ording et al. (2006) the HID is replaced by the disk fraction luminosity diagram or DFLD, see Fig. 1.9. Hardness is replaced by the powerlaw luminosity over the total luminosity and the intensity is replaced by the total luminosity. The used sample from K¨ording et al. (2006) is very similar to the conventional HIDs from XRBs. The same models that are used to fit hard-state XRBs can also be used to fit low-luminosity AGN (LLAGN), see for example Connors et al.

(2017).

Figure 1.8: Three hardness-intensity diagrams for different compact objects from Fender

(2010). Hardness is defined as the ratio of fluxes in two X-ray bands, where high hardness (A in the diagram) indicates a nonthermal spectrum and low hardness (B in the diagram) indicates a thermal spectrum. The sources travel counterclockwise around the diagram and spend most of their time in the lower right corner, where they are faintest. For neutron stars and black holes hardness is used on the x-axis while for white dwarfs a definition analogous to the DFL is used: powerlaw luminosity over total luminosity.

(20)

Figure 1.9: The DFLD from K¨ording et al. (2006). This diagram is similar to the HID of stellar mass compact objects. Powerlaw luminosity over total luminosity is plotted versus total luminosity. R is the radio loudness, the ratio of radio over optical luminosity. The data is binned and only bins with more than 10 counts are shown. The gap between the LLAGN on the lower right and the quasars on the top is caused by various selection effects in the used survey. As you can see, LLAGN statistically show a similar pattern as XRBs even with the quiescent state in the lower right being the most heavily populated. Another indication of a similarity between light and heavy black holes is the fundamental plane of black hole accretion (Merloni et al., 2003; Falcke et al., 2004), see Fig. 1.4. Galactic black holes in the low-hard state and their supermassive analogs show the same correlation between their radio and X-ray luminosity when corrected for mass. This is a very strong correlation, suggesting that the physics of sub-Eddington accretion flows and jet launching are similar across the mass scale. Once normalized for mass, the radiative processes at the inner accretion flow (producing X-rays) and the radiative processes in the core of the jet (producing the bulk of the radio emission) seem to be linked in the same way for black holes of all masses. Theoretically, this can be supported by the fact that jet launching physics might be scale-invariant, more specifically its equations can be scaled by RG = GM_c2 (Falcke & Biermann,1995;Heinz & Sunyaev,2003;Markoff et al.,2003), where

RG is the gravitational radius, see Sect. 3.5 for a further explanation of scale-invariance.

Although the jet launching plasma physics that dominates the dynamics might be scale-invariant, the associated radiative processes are not, hence the mass correction in the fundamental plane of black hole accretion. There is however one anomalous source, Sgr A*. Our very own SMBH lies below the correlation most of the time, while strong X-ray flares boost it within a few standard deviations of the fundamental plane of black hole accretion (Markoff, 2005). The statistics suggest that flares could be observed that put Sgr A* exactly on the fundamental plane of black hole accretion, at which point Sgr A* might start to follow it, for example by launching radiatively stronger jets. It is unlikely that Sgr A* is the one LLAGN that has a fundamentally different accretion/ejection mechanism. It might be that Sgr A* is in some intermediate, quiescent phase. There is some evidence to suggest this is that case. The Soltan relation (Soltan, 1982), a method

(21)

of estimating SMBH mass that works very well for its simplicity, assumes that all mass in SMBHs is accreted in gas. At the current accretion rate of Sgr A* it would take many Hubble times to get up to its measured mass. By looking at reflected hard X-rays from Sgr B2 (see Fig. 1.3) and other molecular clouds one can infer an X-ray luminosity 106

times brighter than a hundred years ago for Sgr A* (Revnivtsev et al., 2004;Ponti et al.,

2010; Clavel et al., 2013). Lastly, a few years ago the Fermi bubbles were discovered, two γ-ray emitting teardrop-shaped bubbles on either side of the galaxy, extending for 8-10 kpc, but originating from within a 100 pc zone near Sgr A* (Su et al.,2010). Many theories about the formation of these bubbles involve recent jet or general outflow activity of Sgr A* (Zubovas et al.,2011;Guo & Mathews, 2011;Cheng et al., 2011).

Figure 1.10: The fundamental plane of black hole accretion. There is a powerlaw relation between radio and X-ray luminosity for all low-luminosity black hole sources when cor-rected for mass. The only exception is Sgr A*, even its most powerful flares shown here can not boost its X-ray emission to the to sufficient heights. The light green datapoints show other black hole sources. Taken fromMarkoff (2005).

1.5 The Event Horizon Telescope

Out of all the discovered SMBHs, our own is by far the closest by a factor of 103 ( Fal-cke et al., 2000). Even when taking into account that SMBH masses can go above 1010 _{and a black hole’s radius scales linearly with its mass, Sgr A* also subtends the}

(22)

biggest angular size on the sky, with the nearby LLAGN M87 a close second. Its angu-lar size is also many times angu-larger than those of Galactic stelangu-lar-mass black holes, mak-ing Sgr A* the perfect candidate for black hole imagmak-ing. The angular size is given by ∠ = RH/D ≈ 0.2 nanoarcsec M_MBHkpc_D , giving microarcsecond scales for Sgr A* but

nanoarcsecond scales for all other Galactic black holes and most other extragalactic black holes (Goddi et al., 2017). The Event Horizon Telescope (EHT) (Doeleman et al., 2008,

2009; Fish et al., 2011) is going to actually resolve the black hole by looking at the accreting gas surrounding Sgr A*. From simulations one expects to see the black hole “shadow”, an emissionless patch on the event horizon (Falcke et al., 2000), see Fig. 1.5. The shape of this shadow can be used to test alternative theories of gravity in the strong limit. This “shadow” is thought to have an angular size of ∼ 50 microarcseconds. Most of the observations will be done at 230 GHz or 1.3 mm, where the angular resolution is ∼ 22 microacrseconds. At longer wavelengths, line-of-sight scattering would dominate any measurement and the gas becomes optically thick, potentially obscuring the shadow (Bower et al., 2004). Shorter wavelengths are technically very difficult, as the accuracy of the measurements needs to be a small fraction of the wavelength. The next phase of EHT observation can hopefully be done at 0.8 mm.

The results of these observations will need to be checked against theoretical predictions. The usual way this is done is by simulating the relevant fluid dynamics in general rela-tivity, then calculating the radiation properties this would produce. This profile is then fed through a scattering screen calculating to account for the ionized matter between us and Sgr A*. Finally, using VLBI simulations the actual image produced by EHT can be predicted. The effect of the radiation profile calculation and scattering screen can be seen in Fig. 1.5. Other things can also be calculated in a similar manner, like the spec-trum, lightcurve or polarization profile. The latter can be used to distinguish different magnetic field line configurations and as such provides another method to distinguish the jet from the disk. All of these different diagnostics can be fit to the eventual EHT data to determine which processes are exactly at play.

1.6 Introduction to simulations and research question

There are multiple ways to make predictions about these types of flows. One way is the analytical, the pen-and-paper approach, as in Shakura & Sunyaev (1973); Narayan & Yi (1994). In this method either explicit equations are derived or implicit equations are solved by hand in order to make predictions. Inevitably, one has to make many simplifica-tions to be able to get a solvable system. Another is semi-analytical, where the amount of simplifications are smaller because equations are solved by the computer. This may sound similar to a simulation, however, it is useful to distinguish the two. In semi-analytics you solve for a small number (∼ 10) of physically distinct parameters. In a simulation or numerical physics, you solve for any number of physically similar parameters, for exam-ple discretized in either space or time. By doing this, very little simplification of your equations is required. The results are however much harder to interpret because you now have millions to billions of values that need to be made sense of.

For some problems, analytic methods are adequate as can be seen by the huge amount of physics already done before the advent of the computer. For the specific type of problem we will be addressing in this thesis, the behaviour of gas in the proximity of a black hole, the equations are highly non-linear and are best described by simulations. The differ-ential equations are strictly local in space (only involves directly surrounding points) so

(23)

we can chop up the space around the block hole in cells and solve the equations between neighbouring cells. After the equations are solved, they are evolved in time and this process is repeated.

In the limit of infinite spatial and temporal resolution, this allows an accurate repro-duction of reality as long as the relevant physics is taken into account. However, this infinite resolution is never achieved due to limited machine precision and computational speed. In order to maintain stable solution and speed up the simulation other numerical methods are used, see Chap. 3. Including the relevant physics is not always straight-forward either. For the simulations in this thesis, the equations we will be solving are the ideal general relativistic magnetohydrodynamics (GRMHD) equations. These equa-tions describe how gas moves in a curved spacetime. The assumpequa-tions of these equaequa-tions are valid for astrophysical plasmas, see Chap. 2. In these equations a lot of physics is excluded, like resistivity, radiation or particle acceleration. Inclusion of these effects is obviously important for the observational signal but even when worrying only about the large-scale dynamics, these effects can be important, even for ADAFs where the very little radiation is produced (Dibi et al.,2012). The effects of radiation are also nonlinear in the simulations: a change in dynamics caused by the loss of internal energy through radiation can change the radiation profile.

Radiation is a manifestly non-local effect. Depending on the medium, photons can travel great lengths before interacting. In two limits, however, radiation can be treated as a local effect: for plasmas so tenuous that any produced radiation escapes the system im-mediately (optically thin) or plasmas so dense that all radiation interacts many times inside a single cell (optically thick), allowing again for the possibility of the local trans-port of a radiation field.

As stated before, the inclusion of radiation also allows us to compare simulation results with observations. In astronomy as a whole, the only way to “observe” in the most general sense is to detect either radiation, cosmic rays (highly energetic massive particles), neu-trinos and gravitational waves. Because of the large effective spatial resolution of particle experiements, it is difficult to locate the source of neutrinos, cosmic rays and gravitational waves, which leaves radiation as the main source of observational constraints for specific sources.

Ray tracing, a technique employed in the making of animated pictures by for example Pixar, allows for the generation of images of 3D objects. This can also be applied to GRMHD simulations (see Fig. 1.5) and can also be used to generate the temporal and spectral signature, see e.g. Dibi et al.(2012); Drappeau et al. (2013).

In this thesis, I willll be combining the GRMHD equations with an optically thin radia-tion field to determine the accreradia-tion dynamics and hopefully the origins of the flares of Sgr A*. This is done with the GRMHD code HAMR (Liska et al.,2018) a state-of-the-art, highly parallel, GPU-accelerated code that allows us to run bigger and longer simulations than ever before. The accretion flow around Sgr A* has been investigated before with GRMHD simulations. In most studies, radiation is not taking into account in the energy conservation equation (Mo´scibrodzka et al., 2009; Dexter et al., 2009, 2010). If we can observe these systems, some energy must be lost from these systems. This could be a negligible amount, but this might not always be the case as we will see in Chap. 4. I will try to show that subtracting the radiative energy from the total energy is important for the dynamics and observational appearance of Sgr A*. One study did self-consistently cool their simulations, but with a different code (COSMOS++ Anninos et al., 2005), with-out pair-production and at a much lower resolution, which determines the duration of

(24)

steady accretion and the amplification of the magnetic field by the MRI (Dibi et al.,2012;

Drappeau et al., 2013). This all makes a direct comparison to our results Our highest resolution run is a factor of 50 higher than theirs and we can run simulations for a much longer period. In (Drappeau et al., 2013) no images were generated, which we do. This all combined will allow us to study the dynamics, observational signature and flaring behaviour of Sgr A* better than ever before.

(25)

CHAPTER

2

Theory

2.1 General relativity and the Kerr metric

According to general relativity, the metric of a non-charged, spinning black hole is the Kerr metric (Kerr, 1963). There is observational evidence for a nonzero spin in both stellar-mass and supermassive black holes (Fabian et al., 1989;Zhang et al.,1997; Bren-neman & Reynolds,2006). Nonzero spin is achieved by black hole mergers and the accre-tion of gas with angular momentum of the same sign. The Kerr metric in Boyer-Lindquist coordinates is (Boyer & Lindquist, 1967):

ds2 = −1−2r Σ dt2+Σ ∆dr 2_+Σdθ2₊_r2_+a2₊2ra2 Σ sin θ

sin θ2dφ2−4ra sin θ

2

Σ dφdt (2.1) where (t, r, θ, φ) are the standard spherical coordinates. Σ = r2+a2cos θ2, ∆ = r2−2r+a2

and a is the dimensionless spin parameter, which ranges from 0 to 1. When considering accretion flows, it is common to let it range from -1 to 1 as to allow for prograde and retrograde disks or photon orbits, but for the metric itself is symmetric in this range. Spacetime itself actually rotates around a spinning black hole. Imagine a test particle at infinity, where gµν = ηµν, with no angular momentum with respect to the black hole. As

the particle approaches, spacetime is no longer flat and the test particle’s reference frame starts to rotate with Ω = −gφt

gφφ =

2rac

Σ(r2_+a2_)+2ra2_{sin θ}2. This effect is called frame dragging

and means that everything that approaches the black hole will experience a “force” along the rotation of the black hole (Misner et al.,1973).

A few surfaces can be distinguished from this metric. When ∆ → 0, grr → ∞. This is

also known as the ”event horizon”, located at rH = 1+

√

1 − a2 _{(we pick the higher-valued}

solution because this reduces to the Schwarzschild solution when a → 0). Nothing can leave the volume the horizon encompasses, thus providing a natural boundary condition for the radial coordinate. Another interesting surface exists outside the event horizon. When gtt switches sign, timelike paths can only stay timelike if they rotate in the same

direction as the black hole with a minimum velocity of Ω. Any retrograde flow is no longer possible. This surface is referred to as the “ergosphere”, located at rE = 1 +

√

1 − a2_{cos θ}

(we choose the highest-valued solution because the lower valued solution lies within the event horizon).

For the stability of the simulations and to make use of the aforementioned natural bound-ary condition, it would be useful to find a coordinate system in which the event horizon is no longer a coordinate singularity so we can integrate across it. As it turns out, this is possible with two coordinate transformations, following Meier (2012):

(26)

dt = dt0−4r

∆dr (2.2)

and

dφ = dφ0− 2a

∆dr (2.3)

These transformations give the line element ds2 = −(1 − 4r Σ2)dt 02 + (1 + 4r Σ2)dr 2 +Σ2dθ2+ (Σ2+ a2(1 + 2r Σ)) sin θ 2_dφ02 +4r Σ2drdt − 4ar sin θ2 Σ2 dφ 0 dt −2a(1 + 2r Σ2) sin θ 2_drdφ0 (2.4)

which only has the curvature singularity at Σ = 0, at the cost of 3 extra nonzero indepen-dent components. Notice that energy and angular momentum in the φ−direction are still conserved. This is the metric is called Kerr-Schild and used inside HAMR. For numerical reasons usually a simple (diagonal) transformation is made from (t0, r, θ, φ0) to (x0, x1,

x2, x3) where t = x0, r = ex1, θ = πx2+1₂(1 − h) sin 2πx2 and φ = x3. Collimation of the

grid toward the equatorial plan can be set by h = [0, 1]. 2.2 Magnetohydrodynamics

Plasmas appear throughout the astrophysics. From the low-density gas in the ISM or the solar wind to the dense plasmas of stellar interiors or the universe before recombination. Because plasmas are made up of charged particles, they can sustain magnetic fields and as such behave differently than ordinary gases. When the scales of neutral gases are much longer than typical microscopic scales of their constituent molecules or atoms, we can describe the gas by the macroscopic theory of fluid dynamics. So too, for plasmas, where the macroscopic theory is called magnetohydrodynamics or MHD. The equations of MHD are similar to those of fluid dynamics and because of this the theories share most assumptions.

2.2.1 Assumptions

The assumptions of hydrodynamics

There is only one basic assumption that hydrodynamics needs in its most general form, the rest of the assumptions are only needed to simplify the equations. This basic assumption is the continuum assumption: molecular properties are averaged on scales much smaller that characteristic length scales for the fluid. This is usually expressed as the Knudsen number Kn = _Lλ where λ is the mean free path of the molecules and L is some character-istic length scale. The same can be done for time scales.

Another assumption we will use is that of a perfect fluid, where we ignore all viscous, shear and conduction effects. We can do this if the particles can be approximated as point-like: the size of the particles is small compared to the mean separation (n−1/3 where n is the particle density).

(27)

The assumptions of MHD

Similar to the continuum assumption of fluid dynamics above, MHD has the quasi-neutrality assumption. Even though the plasma is made up of charge particles, each volume element is assumed neutral. This is because any net charge produces such a large electrostatic force that it quickly disappears. These net charges get removed on time scales of _ω1

p where ωp =

q

4πnee2

me is the plasma frequency, ne is the electron density, me is

the electron mass and e is the electron charge.

In most astrophysical scenarios, the plasma has a very high conductivity. MHD with infinite conductivity is called ideal MHD. The conductivity is σ = e2ne

meνc where νc is the

collisional frequency.

2.2.2 Magneto-rotational instability

An important consequence of ideal MHD is Alven’s theorem or the flux freezing theorem. This theorem says that magnetic field lines are “frozen” in the plasma: the plasma and magnetic flux move together. We start by deriving Ohm’s law for plasmas. A covariant version of a part of the derivation below is given in chapter 9.2 of Meier (2012). Ohm’s law is not clearly applicable to plasmas, especially without external field. To derive the generalized Ohm’s law, we start by writing down the momentum equation for a multi-component fluid. We do not derive this expression rigorously, but the usual shape of a Lagrangian derivative on the LHS and force terms on the RHS should instill some confidence.

msns

d~us

dt = ens( ~E + ~us× ~B) − ∇ ~P + ~Fc (2.5) where the subscript s denotes the particle species, ~P is the pressure, ~us is the velocity, ~E

is the electric field, ~B is the magnetic field and ~Fc denotes the force-vector holding the

collisional terms with the other species in the fluid. Summing the plugged-in equations for ions and electrons and noting that me<< mi and that ~ui= ~u +_mme~j

ien and ~ue = ~u −

~j en

for the center-of-mass velocity ~u = mi~ui+me~ue

mi+me gets us to me ne2 ∂~j ∂t = ~E + ~u × ~B − 1 ne(~j × ~B) − 1 ne∇Pe− 1 σ~j (2.6) where ~j is the current. The only density that appears is n since quasi-neutrality implies that ni = ne. The term on the left is known as the electron inertia for obvious reasons.

The first and second term on the right give the Lorentz force. The third term is known as the Hall effect. The fourth term is the electron pressure and the fifth term is a parametrization of the collisional force. As we will be performing global simulations, we will assume our fluid varies weakly on typical plasma time and distance scales. Using these assumptions, we can neglect the inertial, Hall and pressure term. Equation2.6now reduces to

~

E + ~u × ~B = 1

σ~j (2.7)

In the limit of perfect conductivity (σ → ∞) we get the ideal Ohm’s law ~

(28)

For completely ionizated plasma νc = 2.91 × 10−6neln ΛT−3/2 s−1 where ln Λ ≈ 10 is

the Coulomb algorithm. For typical values for our accretion flow this gives σ ≈ 1027_,

justifying our infinite conductivity assumption. Now we can start to derive the flux freezing theorem. The magnetic flux Ψ through a curve C spanned by a surface S is given by Ψ = Z S ~ B · dS (2.9)

We take the time derivative dΨ dt = Z S d ~B dt · dS + Z C ~ B · ~u × dl (2.10) The first term on the right we can rewrite using one of Maxwell’s law ∂ ~_∂tB = −∇ × ~E and the second term we can rewrite Stokes’ theorem, which gives

dΨ dt = − Z S ∇ × ~E · dS − Z S ∇ × (~u × ~B) · dS = − Z S ∇ × ( ~E + ~u × ~B) · dS = 0 (2.11) where we use the ideal Ohm’s law in the last step. So the flux through any curve moving with the plasma is conserved.

In a differentially rotating accretion flow, when two particles at the same radius, linked by a magnetic field line, get radially perturbed, the inner particle feels a decelerating force and the outer particle feels an accelerating force caused by magnetic stresses, see Fig.

2.2.2. This transport angular momentum to the outer particle and causes turbulence. This effect is called the magnetorotational instability and is thought to be the main driver of accretion.

Figure 2.1: A schematic of the MRI. Two particles in a differentially rotating flow con-nected by a spring exercise forces on each other. Angular momentum is transferred from the inner particle to the outer particle. The spring in this case is a magnetic field line. Credit: Nick Murphy.

(29)

2.3 General relativistic magnetohydrodynamics 2.3.1 Conserved fluid quantities

For a semi-rigorous derivation of the GRMHD equations as they are solved in GRMHD codes, we follow closely Gammie et al. (2003). Where warranted, we give justifications of assumptions and extra steps.

The conservation of particle number for a single-species fluid is given by

∇µJµ = 0 (2.12)

where Jµ _{= ρu}µ _{is the density current, where ρ is the rest-mass density and u}µ _{is the}

plasma velocity. After a choice of coordinates this becomes ∂ ∂xµJ µ_{+ Γ}λ λµJµ = ∂ ∂xµ( √ −gJµ_{) = 0} _(2.13)

where Γµ_νλ is the Christoffel symbol. This is true because Γλ_λµ = 1₂gδν∂µgδν = ∂µln g = 1

√ −g∂µ

√ −g.

The conservation of stress-energy (Bianchi identities) is given by

∇µTµν = 0 (2.14)

where ∇µ is the covariant derivative associated with the spacetime metric gµν and Tµν is

the stress-energy tensor. Mathematically this can be shown to be true since the covariant divergence of the Einstein tensor Gµν disappears and these two quantities can be related through the Einstein equation:

Gµν = 8πTµν (2.15) ∇µGµν = 0 (2.16)

Due to the symmetry of Tµν _{if we choose a coordinate system we can write equation}_2.14

as ∂ ∂xµT µν + Γµ_µλTνλ + Γν_µλTµλ = √1 −g ∂ ∂xµ( √ −gTµν) + Γν_µλTµλ = 0 (2.17) for the same reason given at equation 2.13.

2.3.2 The energy-momentum tensor

Any symmetric (0,2) tensor (in this case Tµν) can be decomposed into two scalar fields,

one vector field and one tensor field with respect to any time-like vector field. In this case that timelike vector field uµ. We obtain the decomposition by projecting Tµν along

uµ _{and the projection tensor h}

µν = uµuν + gµν (so that hµνuν = 0). Tµν = Tκλuκuλuµuν− 2hκµTκλuλuν+ 1 3h κλ_T κλhµν+ (hκµh λ νT(κλ)− 1 3h κλ_T κλhµν) (2.18)

We now have scalar fields Tκλuκuλ and 1₃hκλTκλ, vector field −2hκµTκλuλ and traceless

symmetric tensor field (hκ

µhλνT(κλ)− 1₃hκλTκλhµν).

The usual physical interpretation of Tµν _{for a fluid is the flux of the µth component of}

the energy-momentum vector through a hypersurface perpendicular to eν. This means

(30)

• For µ = ν = 0, we propagate energy density through time, so T00 _{= , where is}

the energy density.

• For µ = 0, ν = i, we propagate energy density through space, so T0i _{= T}i0 _{= q}i_,

where qi _{is the momentum density vector. Alternatively, for µ = i, ν = 0 we}

propagate momentum through time, leading to the same conclusion.

• For µ = ν = i, we propagate momentum along parallel to its direction, so Tii _{= p}i_,

where pi _{is the pressure density vector.}

• For µ = i, ν = j, i 6= j, we propagate momentum perpendicular to itself, so Ti6=j ₌

σij, where σij is the symmetric traceless shear stress density tensor.

Using these physical interpretations, we can now identify = Tκλuκuλ, qµ= −2hκµTκλuλ,

p = 1 3h

κλ_T

κλ and σµν = (hκµhλνT(κλ)− 1₃hκλTκλhµν), resulting in

Tµν = uµuν+ qµuν + phµν+ σµν (2.19)

Now we go to the “perfect fluid” approximation, where we ignore all viscous, shear and conduction effects such that qµ = 0 and σµν = 0 and truncate the moment equations at

second order, reducing our energy-momentum tensor to

Tµν = uµuν+ phµν = (ρ + u + p)uµuν + pgµν (2.20)

where ρ is the rest-mass energy and u is the internal energy density. The usual way equation2.20was constructed is by noting that for a fluid not to be affected by anisotropic (.eg. viscous) effects, its energy-momentum tensor needs to be diagonal where T11 = T22_{= T}33 _or

T = diag(, p, p, p) (2.21) If we find a tensorial equation that reduces to this expression in our rest-frame we are done. A relatively easy guess is

Tµν = ( + p)uµuν+ pηµν (2.22) where ηµν _{= −1, 1, 1, 1 is the Minkowski metric. Because u}µ _{= (−1, 0, 0, 0) in the}

rest-frame of the fluid, equation2.22 reduces to equation2.21. Generalizing the metric brings us back to equation 2.20.

Equations 2.17 and 2.13 give us 5 equalities for 6 unknowns (uµ_{, ρ, u and p). We can}

close this system of equations by introducing an equation of state of the form

p = p(ρ, u) (2.23) The most commonly used equation of state is that of an ideal gas

p = (γ − 1)u (2.24) where γ is the adiabatic index. For a non-relativistic gas, γ = 5/3 and for an ultra-relativistic gas, γ = 4/3.

(31)

2.3.3 Conserved electromagnetic quantities

Until now, we have ignored the “M”-part in “GRMHD”. In this Sect. we will explore the effects on coupling the equations we have already established with those of Maxwell. To remind everyone of their undergraduate physics, the homogeneous Maxwell equations are

∇ · ~B = 0 (2.25) ∇ × ~E −∂ ~B

∂t = 0 (2.26)

and the inhomogeneous Maxwell equations are

∇ · ~E = ρ (2.27) ∇ × ~B − ∂ ~E

∂t = ~J (2.28) These can be reduced to two equations

∂µFµν = Jµ (2.29)

and

∂[µFνλ] = 0 (2.30)

which can also be written as

∂µGµν = 0 (2.31)

where we introduced the Faraday tensor Fµν _{and its Hodge dual G}µν ₌ 1 2

µνκλ_F

κλ. The

contravariant representation of Fµν takes the form Fµν =

   

0 −Ex −Ey −Ez

Ex 0 −Bz By Ey Bz 0 −Bx Ez −By Bx 0     (2.32) Because of the homogeneous Maxwell equations 2.25 and 2.26 we can write

~ B = ∇ × ~A = 0 (2.33) ~ E = −∇φ −∂ ~A ∂t (2.34) ~

A and φ are not uniquely defined so we can choose our gauge for which we will choose the Lorentz gauge

∂φ

∂t + ∇ ~A = 0 (2.35) Grouping our new scalar and vector potential into a 4-vector we get

(32)

Aµ = (φ, ~A) (2.36) and it is straightforward to show that

Fµν = ∂µAν − ∂νAµ (2.37)

Analogous to the fluid energy-momentum tensor of equation 2.20 we can construct the electromagnetic energy-momentum tensor T_EMµν using less-than-rigorous arguments based on undergraduate physics.

• For µ = ν = 0, we propagate energy density through time, so T00 _{= , where is}

the energy density. From electromagnetism we know = E2_{+ B}2_.

• For µ = 0, ν = i, we propagate energy density through space, so T0i _{= T}i0 ₌

qi, where qi is the momentum density vector. Alternatively, for µ = i, ν = 0 we propagate momentum through time, leading to the same conclusion. From electromagnetism we know this as the Poynting flux, or qi _{= ~}_{E × ~}_B.

• For µ = i, ν = j, we need a tensor whose divergence gives the force vector and we can identify this with the Maxwell stress tensor σij =

δij

2 (E

2_{+ B}2_{) + E}

iEj+ BiBj.

In Minkowski space the arguments lead to the following tensor Tµν =     E2+ B2 EyBz− EzBy ExBz− EzBx ExBy− EyBx EyBz − EzBy −σxx −σxy σxz EzBx− ExBz −σyx −σyy σyz ExBy − EyBx −σzx −σzy σzz     (2.38) A more rigorous derivation can be done by employing the Lorentz 4-force density fµ. As this is the time-derivative of pµ, the spatial entries are given by the Lorentz force density

~

f = ρ ~E + ~J × ~B. The temporal entry is given by the time derivative of the energy or ~

f · ~u. This means we can write the 4-force density as

fµ= FµνJν (2.39)

Using equation 2.29 this becomes

fµ = Fµν∂λFλν = ∂ν(FνλFµλ) − Fνλ∂νFλµ (2.40)

The last term in equation2.40 can also be written as Fνλ∂νFλµ = 1 2F νλ_(∂ νFλµ+ ∂µFνλ) = − 1 2F νλ_∂ µFλν = 1 4∂µ(F νλ_F νλ) (2.41)

where in the first step we used the anti-symmetry of Fµν_{, in the second step we used}

equation2.30and in the third step we used the anti-symmetry again. Plugging this result into equation 2.40 gives

fµ= Fµν∂λFλν = ∂ν(FνλFµλ) −

1 4∂µ(F

νλ_F

νλ) = −∂νTµEMν (2.42)

where we recover our electromagnetic energy-momentum tensor again, this time formu-lated in a general spacetime as

(33)

T_EMµν = FµλF_λν −1 4g

µν

FλκFλκ (2.43)

With very laborious checking one can verify that this is equivalent to 2.38. You may notice that the covariant version of equation 2.42 is given by

fµ = ∇νT_EMµν (2.44)

but do not we also have equation 2.14? Peace is restored to the universe by realizing T_EMµν is only a part of the energy-momentum tensor, because the existence of ρ and ~J contained in fµ _{means we need matter, in our case given by T}µν

F . Equation 2.14 is thus

a kind of local version of Newton’s third law or conservation of momentum. This is the first time we see the deep relation between the fluid and the electromagnetic field, the relation that is so important in GRMHD.

We can rewrite equation 2.8 in the fluid frame to a form compatible with our other equations

uµFµν = 0 (2.45)

one of the most important results for ideal magnetohydrodynamics.

Because of the ideal MHD condition it is useful to define the magnetic 4-vector

bµ = µνκλuνFλκ (2.46)

which one can check reduces to b0 _{= 0, b}i _{= ~}_{B in the rest-frame of the fluid. Note that}

equation 2.45 implies that bµuµ = 0. We can use this to invert equation2.46, giving

Fµν = µνκλuκbλ (2.47)

and its dual

Gµν = bµuν − bνuµ (2.48) Using this last equation together with2.31 and the identity used in equation2.13 we get the following equations

∂i( √ −gBi) = 0 √ −g∂tBi+ ∂i[ √ −g(bj_ui_{− b}i_bj_{)] = 0} (2.49)

where Bi is implicitly defined as bt = Biuµgiµ, bi = Bi+b

t_ui

ut . The first equation is the

general relativistic equivalent of ∇ · ~B = 0 and the second equation is the induction equation. Using equation 2.47 we can reduce the electromagnetic energy-momentum tensor (equation 2.43) to

T_EMµν = b2uµuν +1 2b

2_gµν _{− b}µ_bν _(2.50)

Our total energy-momentum tensor, to conclude is

T_{M HD}µν = T_Fµν + T_EMµν = (ρ + u + p + b2)uµuν + (p +1 2b

2_)gµν _{− b}µ_bν _(2.51)

Now we have established our basic equations (equation2.13,2.17,2.49 and2.24) and our stress-energy tensor T_{M HD}µν = T_Fµν+ T_EMµν . How these equations are solved, we will discuss in chapter3.

(34)

2.4 Radiation

In order to cool the simulation, we need to subtract the radiated energy from the internal energy of the simulation. In exchange for the luxury of a low optical depth everywhere (see Sect. 2.4.1) we have to deal with non-blackbody radiation, see Tab. 2.1 for a summary of radiative and pair-producing processes in non-magnetized plasmas, to which we only need to add synchrotron radiation. Luckily, in Sgr A* we can ignore bremsstrahlung at radii < 103 _R

G (Narayan et al., 1998; Mo´scibrodzka et al., 2009), we do however take

bremsstrahlung into account for the cooling so our simulations can be used at some later point when larger simulation sizes might be considered. For the spectra and images we do not consider bremsstrahlung. The two radiative types we need to worry about for that are synchrotron radiation and Compton scattering. Compton scattering will be only used for the synchrotron radiation field and will be described by a correction factor η describing the energy amplification per synchrotron photon by scattering. The total cooling rate is given by

q_total− = ηq−_sync+ q_brems− (2.52) This is substracted from the source term belonging to the energy equation. Before we discuss the produced radiation, we first need to explore a key factor in our simulations: low optical depth.

2.4.1 The implications of low optical depth

In an ADAF, we can estimate the optical depth for Thomson scattering by τe = σTneH

where σT is the Thomson cross-section and H ∼ R is the scale-height of the flow (Esin

et al., 1996). H typically lies in the range 1010 _{to 10}13 _{for our grid. In our}

highest-accretion rate simulation, ne,max ≈ 1012, providing an upper limit of τe < 0.06 but

usually τe << 10−3. For other radiation types the optical depth is similarly low for

the frequencies that cool the most. Thus, we can assume that any radiation produced leaves the simulation without interacting with the fluid, meaning that radiation only removes energy. Low optical depth also means that the only pair-producing interaction that matters is ee ↔ eee+_e− _{while electron-positron pair annihilation still dominates the}

reverse reaction (Svensson, 1982a).

2.4.2 Approximations for Sagittarius A*

We assume the electrons get thermalized at some fraction of the proton temperature Te = fTTp (2.53)

fT = {0.1, 0.3, 1} (2.54)

where Te is the electron temperature, Tp is the proton temperature and fT is the ratio

as in Mo´scibrodzka et al. (2009); Esin et al. (1996); Dibi et al. (2012); Fragile & Meier

(2009). This is our main assumption, it allows for cheap analytic cooling functions that are imperative for GRMHD codes. It has been shown to agree with more sophisticated methods that track electron distributions (Dibi et al., 2012). This is an approximation of a weak coupling between protons and electrons. The protons heat up through viscous dissipation and in turn heat up the electrons, while the electrons cool through radiative processes and in turn cool the protons. There will always be some coupling between the

The impact of radiation on the dynamics and appearance of the supermassive black hole in the Galactic center

University of Amsterdam

MSc Physics and Astronomy

Track: Astronomy & Astrophysics

Master Thesis

The impact of radiation on the dynamics and

appearance of the supermassive black hole in

the Galactic center

“Have you seen Interstellar?”

by

David Jacob van Eijnatten

10645462

July 29, 2018

60 ECTS

2017 - 2018

Supervisor:

prof. dr. S.B. Markoff

Daily supervisor:

Koushik Chatterjee, MSc

Examiners:

prof. dr. S.B. Markoff

prof. dr. M.B.M. van der Klis

Contents

List of Figures

List of Tables

CHAPTER

1

Introduction

log[ν(Hz)]

log[

ν

L

(erg s

)]

Sgr A

CHAPTER

2

Theory