Blog Post 37, Illustris Simulation

The Illustris Simulation is a numerical cosmological simulation which includes the effects of dark energy, dark matter, baryonic matter, and the relevant underlying physics (radiation, fluid dynamics, gravity, etc).

To look at one aspect of large scale structure in this simulation, we will consider the halo mass function (HMF). The HMF is just a fancy way of counting the humber of halos within some mass bin. The halo is the extended dark matter structure which surrounds galaxies individually, and galaxy clusters.

Create a histogram of the data with bins that have a width of 0.5 in log(M). Are the low mass halos more numerous, or are the high mass halos?

We first used the Illustris Explorer tool to examine subhalos within a 400.0 kpc/h radius. We performed this search several times, with each search returning the position, halo mass, and star mass of up to 20 subhalos. We combined the output for several searches to get the following histogram

As you can see in the histogram, low mass halos are significantly more numerous than the high mass halos.

On average, what fraction of the halo mass is the stellar mass? 

On average, the stellar mass makes up around 90% of the halo mass.

Exploring Structure and Reionization in the Illustris Simulation

Next, we compared the structures between the gas density and dark matter density in the simulation. On the large scale (full box), the simulation looks like the following, with gas density on the right and dark matter density on the left:

On the small scale, the simulation looks like the following, again with dark matter density on the left, and gas density on the right:

Compare the Gas Density and the Dark Matter Density on both large (full box) and small scales (single cluster).  Describe how the Gas and Dark Matter are similar/different on each scale and speculate as to why.

As you can see in the photos above, while the gas density and dark matter density both follow similar patters, the dark matter density is much more structured, while the gas density is more diffuse, both on the large scale and on the small scale. This is particularly noticeable on the small scale, as the dark matter appears to be much more clustered than the gas density, which is more evenly spread out.


Which, the gas or the dark matter, is more confined to the filamentary structure and why?

The dark matter is more confined to the filamentary structure, as you can see on the large scale image of the dark matter density. This could possibly be because of the increased density of the dark matter; as dark matter is thought to constitute around 80% of the matter in our galaxy, this greater mass results in a larger gravitational attraction between the dark matter particles than the gas particles.

In most of medium to large galaxies, is the gas densest towards the nucleus or the disk?

In most the medium to large-scale galaxies, the gas appears to be denser towards the nucleus of the galaxy, as shown in the following snapshot:

Are the most massive galaxies in the field or in clusters?

The most massive galaxies are in the cluster. As this is where the majority of gas is concentrated, it makes sense that the galaxies that form within the clusters are larger on average than those in the field.

Next, we watched the following video, which shows the evolution of the dark matter and gas within our universe in the time since the big bang:

http://www.illustris-project.org/movies/illustris_movie_cube_sub_frame.mp4

When do the first stars form in this simulation?

The first stars appear around 0.6 billion years after the Big Bang, when there is a redshift of about 8.

For what approximate redshift range is the rate at which stars are forming fastest?

The stars appear to be forming at the fastest rate when the redshift is between 1.0 and 2.0.

When, in years and redshift, does the gas temperature brighten (go from blue to having green?) This is the beginning of the "Epoch of Reionization" or the end of the "Dark Ages."

We begin to see the first ionized gases - or when the gas temperature first turns green - about 1.0 billion years after the Big Bang, with a redshift of 5.80.

In this simulation, when structures are forming, are smaller structures combining to form larger ones, or are large objects breaking up? Why do you think this is?

In the process of gas formation, we see structures form as large objects are breaking up. This makes sense, as structures gain more matter, they are likely to collapse under their own gravity in a supernovae.

Why do you think structures form along filaments?

As we saw above, dark matter forms along the filamentary structures. As dark matter makes up the majority of matter in the universe, gravitational attraction from dark matter results in the formation of structures along these filaments as well.

Blog Post 36, Worksheet 12.1, Problem 1 and 2(d)

1. Linear perturbation theory. In the early universe, the matter/radiation distribution of the universe is very homogeneous and isotropic. At any given time, let us denote the average density of the universe as \(\bar{\rho}(t)\). Nonetheless, there are some tiny fluctuations and and not everywhere exactly the same. So let us define the density at comoving position r and time t as \(\rho (x, t)\) and the relative density contrast as \[\delta (r, t) \equiv \frac{\rho (r, t) - \bar{\rho} (t)}{\bar{\rho}(t)}.\] In this exercise we focus on the linear theory, namely, the density contrast in the problem remains small enough so we need only consider terms linear in \(\delta\). We assume that cold dark matter, which behaves like dust (that is, it is pressureless) dominates the content of the universe at the early epoch. The absence of pressure simplifies the fluid dynamics equations used to characterize the problem. 

(a) In the linear theory, it turns out that the fluid equations simplify such that the density contrast \(\delta\) satisfies the following second-order differential equation \[\frac{d^2\delta}{dt^2}+\frac{2\dot{a}}{a}\frac{d\delta}{dt} = 4πG\bar{\rho}\delta.\] where a(t) is the scale factor of the universe. Notice that remarkably in the linear theory this equation does not contain spatial derivatives. Show that this means that the spatial shape of the density fluctuations is frozen is comoving coordinates, only their amplitude changes. Namely this means that we can factorize \[\delta(x, t) = D(t)\bar{\delta}(x),\] where \(\bar{\delta}(x)\) is arbitrary and independent of time, and D(t) is a function of time and valid for all x. D(t) is not arbitrary and must satisfy a differential equation. Derive this differential equation. 

We are given the equation in the form \[\frac{d^2\delta}{dt^2}+\frac{2\dot{a}}{a}\frac{d\delta}{dt} = 4πG\bar{\rho}\delta.\] We want to show that the equation for \(\delta\) can be written in the form \(\delta(x, t) = D(t)\bar{\delta}(x),\). To show this, we can use this equation as a substitution for \(\delta\) in above differential equation: \[\ddot{D}(t)\bar{\delta}(x) + \frac{2\dot{a}}{a}(\dot{D}(t)(\bar{\delta}(x)) = 4πG\bar{\rho}D(t)\bar{\delta}(x)\] This simplifies to: \[\ddot{D}(t) + \frac{2\dot{a}}{a}\dot{D}(t) = 4πG\bar{\rho}D(t),\] which is in the same form as our first equation.

(b) Now let us consider a matter dominated flat universe, so that \(\bar{\rho} = a^{-3}\rho_{c,0}\), where \(\rho_{c,0}\) is the critical density today, \(2H_0^2/8πG\) as in Worksheet 11.1 (aside: such a universe sometimes is called the Einstein-de Sitter model). Recall that the behaviour of the scale factor of this universe can be written \(a(t) = (3H_0t/2)^{2/3}\), which you learned in previous worksheets, and solve the differential equation for D(t). Hint: you can use the ansatz \(D(t) \propto t^q\) and plug it into the equation that you derived above; and you will end up with a quadratic equation for q. There are two solutions for q, and the general solution for D is a linear combination of two components: One gives you a growing function in t, denoting it is \(D_+(t)\); another decreasing function in t, denoting it as \(D\_(t)\). 

We are given that \(a(t) = (3H_0t/2)^{2/3}\), and thus we know that \(\dot{a}(t) = H_0 (\frac{3H_0t}{2})^{-1/3}\). We combine these to get the quantity: \[\frac{\dot{a}}{a} = \frac{H_0 (\frac{3H_0t}{2})^{-1/3}}{(3H_0t/2)^{2/3}}\] \[\frac{\dot{a}}{a} = \frac{2}{3t}\] We are also given that \(\bar{\rho} = a^{-3}\rho_{c,0}\). We can substitute in our expression for a above to get a new expression for \(\bar{\rho}\): \[\bar{\rho} = [3H_0t/2)^{2/3}]^{-3}(2H_0^2/8πG) = \frac{1}{6πGt^2}\] Now we can plug these two expressions (for \(\bar{\rho}\) and for \(\frac{\dot{a}}{a}\)) into our differential equation from part (a): \[\ddot{D}(t) + 2(\frac{2}{3t})\dot{D}(t) = \frac{4πGD(t)}{6πGt^2}\] \[\ddot{D}(t) + \frac{4}{3t}\dot{D}(t) = \frac{2}{3t}D(t)\] We can now use the substitution suggested in the hint above, \(D(t) \propto  t^q\) to make a substitution for \(D(t), \dot{D}(t)\), and \(\ddot{D}(t)\): \[q(q=1)t^{q-2} + \frac{4}{3t}qt^{q-1} - \frac{2}{33t^2}t^q = 0\] This simplifies to: \[\frac{1}{3}t^{q-2}[3q62 + q - 2] = 0\] We can use the quadratic equation to solve for our two possible answers for q, \(q = \frac{2}{3}\) and \(q = -1\). We plug these bag into the ansatz given in the question to get: \[D_+(t) \propto t^{2/3}\] \[D\_(t)\propto \frac{1}{t}\] (c) Explain why the \(D_+\) component is generically the dominant one in structure formation, and show that in the Einstein-de Sitter model, \(D_+(t) \propto a(t).\)

As t increases over time, \(D\_(t)\) will get closer and closer to 0, as \(D\_(t)\propto \frac{1}{t}\), and thus the \(D_+\) component is generically the dominant one in structure formation.

In the Einstein-de Sitter model, we know that \(a(t) = (3H_0t/2)^{2/3}\), and thus \(a(t) \propto t^{2/3}\). Thus we have: \[D_+(t) \propto t^{2/3}\] \[a(t) \propto t^{2/3}\] \[D_+(t) \propto a(t)\]
2. Spherical collapse. Gravitational instability makes initial small density contrasts grow in time. When the density perturbation grows large enough, the linear theory, such as the one presented in the above exercise, breaks down. Generically speaking, non-linear and non-perturbative evolution of the density contrast have to be dealt with in numerical calculations.

(d) Plot r as a function of t for all three cases (i.e. use y-axis for r and x-axis for t), and show that in the closed case, the particle turns around and collapses; in the open case, the particle keeps expanding with some asymptotically positive velocity; and in the flat case, the particle reaches an infinite radius but with a velocity that approaches zero. 

Closed case

\[r = A(1 - \text{cos }\eta)\] \[t = B(\eta - \text{sin }\eta)\]

In this plot, t is on the x-axis, and r(t) is on the y-axis. (Scaled to not include the multiplication factors of A and B). As the time increases, the particle at first moves outward some distance r, and then turns around and returns to a distance of r = 0. 

Open case

\[r = A(\text{cosh }\eta - 1 \] \[t = B(\text{sinh }\eta - \eta)\]

The following two plots show the open case with with relatively large and small limits for \(\eta\), respectively, with t on the x-axis, and r(t) on the y-axis:

As time increases, the distance r from the center of the sphere continues to increase. The velocity of the particle is always positive, and thus the particle keeps expanding with some asymptotically positive velocity. 

Flat case

\[r = A\eta^2/2\] \[t = B\eta^3/6\]

In this case, as time increases, while the distance of the particle continues to increase, the velocity of the particle decreases. Thus, the particle reaches an infinite radius but with a velocity that approaches zero.

Blog Post 35, Worksheet 11.1, Problem 3

3. Baryon-to-photon ratio of our universe. 

(a) Despite the fact that the CMB has a very low temperature (that you have calculated above), the number of photons is enormous. Let us estimate what that number is. Each photon has energy \(h\nu\). From equation (1), figure out the number density, \(n_\nu\), of the photon per frequency interval \(d\nu\). Integrate over \(d\nu\) to get an expression for total number density of photon given temperature T. Now you need to keep all factors, and use the fact that \[\int_0^\infty \frac{x^2}{e^x-1}dx \approx 2.4\]

In problem 1, we are given the following equation for the energy density per frequency interval \(d\nu\) of the black body radiation: \[u_{\nu}d\nu = \frac{8πh_p\nu^3}{c^3}\frac{1}{e^{\frac{h_p\nu}{e_BT}}-1}d\nu\] To get the number density, \(n_\nu\), from the above equation, we need to divide by the energy of a photon, \(h\nu\):
\[n_\nu = \frac{u_\nu d\nu}{h_p \nu} = \frac{1}{h_p \nu} \frac{8πh_p\nu^3}{c^3}\frac{1}{e^{\frac{h_p\nu}{e_BT}}-1}d\nu\] This gives us the following formula for \(n_\nu\): \[n_\nu = \frac{8π}{c^3}\frac{\nu^2}{e^{\frac{h_p\nu}{e_BT}}-1} d\nu\] Now we can use u substitution to simplify the integral, with \(u = \frac{h\nu}{kT}\) and \(du = \frac{h\:d\nu}{kT}\). We substitute for \(\nu^3\) and \(d\nu\) to get: \[n_\nu = \frac{8πk^3T^3}{c^3h^3} \frac{u^2}{e^u-1}du\] Now we integrate over du: \[n_\nu = \frac{8πk^3T^3}{c^3h^3} \int_0^\infty \frac{u^2}{e^u-1}du\] Using the given value for the integral, we get that \[n_\nu = \frac{19.2 πk^3T^3}{c^3h^3}\]

(b) Use the following values for the constants: \(k_B = 1.38 • 10^{-16} \text{erg K}^{-1}\), \(c = 3.00 • 10^{10} \text{ cm s}^{-1}\), \(h = 6.62 • 10^{27} \text{ erg s}\), and use the temperature of CMB today that you have computed from 2(d) to calculate the number density of photon today in our universe today (i.e. how many photons per cubic centimeter). 

This is just a matter of plugging in the given constants and our value for the temperature, 2.7 K: \[n_\nu = \frac{19.2π(1.38 • 10^{-16})^3(2.7)^3}{(3 • 10^{10})^3(6.62 • 10^{-27})^3} \approx 400 \text{ photons cm}^{-3}\]

(c) Let us calculate the average baryon number density today. In general, baryons refer to protons or neutrons. The present-day density (matter + radiation + dark energy) of our Universe is \(9.2 • 10^{-30} \text{g cm}^{-3}\). The baryon density is about 4% of it. The masses of proton and neutron are very similar \((\approx 1.7 • 10^{-25} \text{ g}\)). 
What is the number density of baryons?

The baryon density is about 4% of the present-day density, or\[ 0.4(9.2 • 10^{-30} \text{g cm}^{-3}) = 3.68 • 10^{-31} \text{g cm}^{-3}\] To find the number density of baryons, we divide this baryon density by the masses of each particle (proton or neutron): \[\frac{3.68 • 10^{-31} \text{g cm}^{-3}}{1.7 • 10^{-25} \text{ g}} = 2.16 • 10^{-7} \frac{\text{baryons}}{\text{cm}^3}\]

(d) Divide the above two numbers, you get the baryon-to-photon ratio. As you can see, our universe contains much more photons than baryons (proton and neutron):

Our ratio is obtained by dividing our answer from part (b) by our answer from part (c), which gives a ratio of approximately \(2.9 • 10^9 \frac{\text{photons}}{\text{baryon}}\).

Blog Post 34, Worksheet 11.1, Problem 1

1. Temperature of the Universe. Remember that, although the universe today is dominated by dark energy and matter (including ordinary matter and dark matter), much earlier on it was dominated by radiation. In this exercise we study the temperature evolution of a radiation dominated universe. 

When the electromagnetic wave is in equilibrium with the environment, its spectrum is uniquely determined by the temperature of the equilibrium. This state is called the blackbody radiation. The spectrum is called the Planck spectrum, named after the physicist who discovered it. The energy density per frequency interval dv of the black body radiation is given by \[u_{\nu}d\nu = \frac{8πh_p\nu^3}{c^3}\frac{1}{e^{\frac{h_p\nu}{e_BT}}-1}d\nu\] where \(h_p\) is the Planck constant, \(k_B\) is the Boltzmann constant, \(\nu\) is the frequency, and T is the temperature. 

(a) How is the equation for \(u_{\nu}d\nu\) difference from the equation for flux given in our previous worksheets?

From previous worksheets, we know that the equation for flux as a function of frequency is as below: \[F_{\nu} = \frac{2πh\nu^3}{c^2}\frac{1}{e^{\frac{h\nu}{kT}}-1}\] Therefore, \(u_{\nu}d\nu = \frac{4}{c}F_{\nu}d\nu\).

(b) Integrate the Planck spectrum over the frequency and figure out how the energy density u of the black body radiation depends on temperature T. Namely, figure out the power n in \(u \propto T^n\). (Since only the functional form of T is important here, in this exercise you do not have to figure out the exact value of the T-independent coefficient a). 

In Worksheet 2.1, we integrated \(F_{\nu}(T)\) over the frequency to determine that \(F(T) = \sigma T^4\), or \(F_\nu \propto T^4\). As we have just shown that \(u_\nu \propto F_\nu\), we can therefore conclude that \(u \propto T^4\).

(c) Remind yourself how the energy density of the radiation dominated universe depends on the scale factor a. 

In previous worksheets, we determined that the energy density of the radiation dominated universe changes with the scale factor a by a factor of \(a^{-4}\).

(d) Combine the two results to see how the temperature T of the universe depends on the scale factor a. Explain why this result implies that the early universe is very hot. 

From parts a, b, and c above, we know both of the following: \[u \propto a^{-4}\] \[u \propto T^4\] We can combine these two proportionalities to see that: \[a^{-4} \propto T^4\] \[T \propto \frac{1}{a}\] As the scale factor increases through time with the expansion of the universe, the temperature thus decreases. Therefore, the temperature of the universe increases as you go back in time, and the early universe was very hot.

Blog Post 32, Worksheet 10.1, Problem 3

3. Observable Universe. It is important to realize that in our Big Bang universe, at any given time, the size of the observable universe is finite. The limit of this observable part of the universe is called the horizon (or particle horizon to be precise, since there are other definitions of horizon.) In this problem, we'll compute the horizon size in a matter dominated universe in co-moving coordinates. 

To compute the size of the horizon, let us compute how far the light can travel since the Big Bang. 

(a) First of all - Why do we use the light to figure out the horizon size?

Light travels at the fastest (known) possible speed in the universe. Thus, the distance that light has traveled determines limit of the observable part of the universe. 

(b) Light satisfies the statement that \(ds^2 = 0\). Using the FRW metric, write down the differential equation that describes the path light takes. We call this path the geodesic for a photon. Choosing convenient coordinates in which the light travels in the radial direction so that we can set \(d\theta = d\phi = 0\), find the differential equation in terms of the coordinates t and r only. 

Setting ds = 0, we know that \[0 = -c^2dt^2 + a^2(t)[\frac{dr^2}{1-kr^2} + r^2(d\theta^2 + \text{sin}^2\theta d\phi^2)]\] After we make the given assumptions, we are left with: \[c^2dt^2 = a^2(t)\frac{dr^2}{1-kr^2}\] Now we rearrange the equation to yield our differential equation: \[\frac{c^2dt^2}{a^2(t)} = \frac{dr^2}{1 - kr^2}\]

(c) Suppose we consider a flat universe. Let's consider a matter dominated universe so that a(t) as a function of time is known (in the last worksheet). Find the radius of the horizon today (at \(t = t_0\). 

(Hint: move all terms with variable r to the LHS and t to the RHS. Integrate both sides, namely r from 0 to \(r_{\text{horizon}}\) and t from 0 to \(t_0\)).

In the case of a flat universe, k = 0, so the RHS of our above differential equation becomes just \(dr^2\). From our previous worksheet, we know that: \[a(t) = a_0\big(\frac{t}{t_0}\big)^{2/3}\] After we take the square root of both sides of our differential equations, we can sub in the above for a(t) to get: \[dr = \frac{cdt}{a(t)} = \frac{cdt}{a_0\big(\frac{t}{t_0}\big)^{2/3}}\] Now we integrate the LHS from r to \(r_{\text{horizon}}\) and the RHS from 0 to \(t_0\): \[\int_0^{r_{\text{horizon}}}dr =  \int_0^{t_0} \frac{cdt}{a_0\big(\frac{t}{t_0}\big)^{2/3}}\] This gives us: \[r_{\text{horizon}} = \frac{3c}{a_0}t_0\]

Blog Post 31, Worksheet 10.1, Problem 2

2. Ratio of circumference to radius. Let's continue to study the difference between closed, flat, and open geometries by computing the ratio between the circumference and radius of a circle. 

a. To compute the radius and circumference of a circle, we look at the spatial part of the metric and concentrate on the two-dimensional part by setting \(d\phi = 0\) because a circle encloses a two-dimensional surface. For the flat case, this part is just \[ ds_{2d}^2 = dr^2 + r^2 d\theta^2.\] The circumference is found by fixing the radial coordinate (r = R and dr = 0) and both sides of the equation (note that \theta is integrated from 0 to 2π). 

The radius is found by fixing the angular coordinate (\(\theta, d\theta = 0\)) and integrating both sides (note that dr is integrated from 0 to R). 

Compute the circumference and radius to reproduce the famous Euclidean ration 2π. 

As stated in the question, we are looking only at the spatial part of the metric: \[ ds_{2d}^2 = dr^2 + r^2 d\theta^2.\] To find the circumference of the circle, we need to set the radius equal to a constant, r = R, with dr = 0. Making these substitutions, we are left with: \[ds_{2d}^2 = R^2 d\theta^2\] After we take the square root of both sides, we have: \[ds_{2d} = R d\theta\] Now, we integrate both sides. We integrate \(\theta\) from 0 to 2π, and we integrate s from 0 to some arbitrary constant C: \[\int_0^C ds_{2d} = \int_0^{2π} Rd\theta\] This gives us our circumference, \[C = 2πR\]

To find the radius, we instead set \(\theta\) and \(d\theta\) equal to zero, as we are integrate outward from the center of the circle: \[ds_{2d}^2 = dr^2\] \[ds_{2d} = dr\] Now, we integrate r from 0 to R, and we integrate s from 0 to some arbitrary constant r: \[\int_0^r ds = int_0^R dr\] \[r = \text{radius} = R\] Now that we have the circumference and radius, we can calculate our circumference to radius ratio: \[\frac{\text{circumference}}{\text{radius}} = \frac{2πR}{R} = 2π\]

(b) For a closed geometry, we calculated the analogous two-dimensional part of the metric in Problem (1). This can be written as: \[ds_{2d}^2 = d\xi^2 + \text{sin}\: \xi^2 d\theta^2.\] Repeat the same calculation above and derive the ratio for the closed geometry. 

Compare your results to the flat (Euclidean) case; which ratio is larger? 

To find the circumference for a closed geometry, we now set \(d\xi = 0\). This is analogous to setting \(dr = 0\) in the above problem, but now we are working in the hyperspherical coordinate system. \[ds_{2d}^2 = d\xi^2 + \text{sin}\: \xi^2 d\theta^2\] \[ds_{2d}^2 = \text{sin}\: \xi^2 d\theta^2\] \[ds_{2d} = \text{sin}\: \xi d\theta\] Now we integrate \(\theta\) from 0 to 2π, and we integrate s from 0 to some arbitrary constant C: \[\int_0^C ds_{2d} = \int_0^{2π} \text{sin}\: \xi d\theta\] \[C = \text{circumference} = 2π\text{sin}\: (\xi)\] To calculate the radius, we set \(d\theta = 0\), leaving us with: \[ds_{2d}^2 = d\xi^2\] \[ds = d\xi\] Now we integrate s from 0 to some radius r, and \(\xi\) from 0 to \(\xi\): \[\int_0^r ds_{2d} = \int_0^{\xi} d\xi\] \[r = \xi\] circumference to radius ratio is thus: \[\frac{\text{circumference}}{\text{radius}} = \frac{2π \text{sin }(\xi)}{\xi}\] The value of \(\text{sin }(\xi)\), will always be between -1 and 1, and thus this ratio will always be less than the ratio for the flat geometry, 2π. 

(c) Repeat the same analyses for the open geometry, and comparing to the flat case. 

For the open geometry (k  = -1), the relevant equation becomes \[ds_{2d}^2 = d\xi^2 + \text{sinh }^2 d\theta^2\] Using the same logic as above, we derive the circumference by setting \(d\xi = 0\) and integrating \(\theta\) from 0 to 2π. To derive the radius, we set \(d\theta = 0\) and integrate \(\xi\) from 0 to \(\xi\). Our ratio thus becomes: \[\frac{\text{circumference}}{\text{radius}} = \frac{2π \text{sinh }(\xi)}{\xi}\] The value of \(\frac{\text{sinh }\xi }{\xi}\) will always be greater than one, so the ratio for the open geometry will always be greater than that of the flat case. 

(d) You may have noticed that, except for the flat case, this ratio is not a constant value. However, in both the open and closed case, there is a limit where the ratio approaches the flat case. Which limit is that?

Both the open and closed case approach the flat case ratio, 2π, in the limit as \(\xi\) goes to 0. 

Blog Post 29, Worksheet 9.1, Problem 2

2. GR modification to Newtonian Friedmann Equation: 

In Question 1, you have derived the Friedmann Equation in a matter-only universe in the Newtonian approach. That is, you now have an equation that describes the rate of change of the size of the universe, should the universe be made of matter (this includes stars, gas, and dark matter) and nothing else. Of course, the universe is not quite so simple. In this question we'll introduce the full Friedmann equation which describes a universe that contains matter, radiation and/or dark energy. We will also see some correction terms to the Newtonian derivation. 

(a) The first Friedmann equation: \[(\frac{\dot{a}}{a})^2 = \frac{8π}{3}G\rho - \frac{kc^2}{a^2} + \frac{\Lambda}{3}\] The second Friedmann equation: \[\frac{\ddot{a}}{a} = -\frac{4πG}{3c^2}(\rho c^2 + 3P) + \frac{\Lambda}{3}\]
In these equations, \(\rho\) and P are the density and pressure of the content, respectively. k is the curvature parameter; k = -1, 0, 1 for open, flat, and closed universe, respectively. \(\Lambda\) is the cosmological constant. Note that in GR, not only density but also pressure are the sources of energy. 

Starting from these two equations, derive the third Friedmann equation, which governs the way average density in the universe changes with time: \[\dot{\rho}c^2 = -3\frac{\dot{a}}{a}(\rho c^2 + P)\] First, we multiply both sides of the first Friedmann equation by \(a^2\): \[\dot{a}^2 = \frac{8πG\rho}{3}a^2-kc^2+\frac{\Lambda}{3}a^2\] Next we take the time derivative of each side. This gives us: \[2\dot{a} \ddot{a} = \frac{8πG}{3}a^2 \dot{\rho} + \frac{16}{3}πG\rho a \dot{a} + \frac{2\Lambda}{3}a\dot{a}\] \[\ddot{a} = \frac{4πGa^2\dot{\rho}}{3\dot{a}} + \frac{8πG\rho a}{3} + \frac{\Lambda}{3}a\] Now we can plug this expression for \(\ddot{a}\) into the second Friedmann equation: \[\frac{4πGa\dot{\rho}}{3\dot{a}} + \frac{8πG\rho}{3} + \frac{\Lambda}{3} = -\frac{4πG}{3c^2}(\rho c^2 + 3P) + \frac{\Lambda}{3}\] Several terms cancel, and we are left with: \[\frac{\dot{\rho}a}{\dot{a}} + 2\rho = -\rho - \frac{3P}{c^2}\] We can now rearrange the equation to give us the third Friedmann equation: \[\dot{\rho}c^2 = -3\frac{\dot{a}}{a}(\rho c^2 + P)\]
(b) If the matter is cold, its pressure P = 0, and the cosmological constant \(\Lambda = 0\). Use the third Friedmann equation to derive the evolution of the density of the matter \(\rho\) as a function of the scale factor of the universe a. You can leave this equation int terms of \(\rho, \rho_0, a, \text{ and } a_0\), where \(\rho_0\) and \(a_0\) are current values of the mass density and scale factor. 

In the third Friedmann equation, P becomes zero, leaving us with: \[\dot{\rho}c^2 = -3\frac{\dot{a}}{a}(\rho c^2)\] A \(c^2\) cancels on both sides, giving us \[\frac{\dot{\rho}}{\rho} = -3\frac{\dot{a}}{a}\] Now we can integrate: \[\int_{\rho_0}^{\rho} \frac{\dot{\rho}}{\rho} = -3\int_{a_0}^{a} \frac{\dot{a}}{a}\] We know that \(\dot{\rho} = \frac{d\rho}{dt}\), and \(\dot{a} = \frac{da}{dt}\). If we substitute both of these in, the dts cancel, and we have: \[\int_{\rho_0}^{\rho} \frac{1}{\rho}d\rho = -3\int_{a_0}^{a} \frac{1}{a}da\] \[\text{ln}(\frac{\rho}{\rho_0}) = -3\text{ln}(\frac{a}{a_0})\] If we raise e to the power of both sides, we get: \[\frac{\rho}{\rho_0} = e^{-3\text{ln}(\frac{a}{a_0})} = (e^{\text{ln}(\frac{a}{a_0})})^{-3} = (\frac{a}{a_0})^{-3}\] This gives us: \[\frac{\rho}{\rho_0} = (\frac{a}{a_0})^{-3} = (\frac{a_0}{a})^3\] \[\rho = \frac{\rho_0 a_0^3}{a^3}\]
Using the relation between \(\rho\) and a that you just derived and the first Friedmann equation, derive the differential equation for the scale factor a for the matter dominated universe. Solve this differential equation to show that \(a(t) \propto t^{2/3}\). This is the characteristic expansion history of the universe if it is dominated by matter. (Hint: When solving this final differential equation, recall that at time t = 0, a = 0 (the Big Bang). At time \(t = t_0, a = a_0 = 1\) (present day).) 

From above, we know that \(\rho \propto a^{-3}\). The first Friedmann equation tells us that \[\dot{a}^2 = \frac{8π}{3}G\rho a^2 - kc^2 + \frac{\Lambda a^2}{3}\] We set k = 0 and \(\Lambda\) = 0, as per parts a and b, giving us: \[\dot{a}^2 = \frac{8π}{3}G\rho a^2\] We can plug in \(\rho = \frac{\rho_0 a_0^3}{a^3}\) to give us: \[\dot{a}^2 = \frac{8πG\rho_0 a_0^3}{3a}\] Take the square root of both sides to get: \[\dot{a} = (\frac{8πG\rho_0 a_0^3}{3})^{1/2}a^{-1/2}\] \[\dot{a} \propto a^{-1/2}\] We know that \(\dot{a} = \frac{da}{dt}\), so we can separate the variables to get: \[\frac{da}{dt} \propto a^{-1/2}\] \[a^{1/2}da \propto dt\] Now we integrate: \[\int a^{1/2}da \propto \int dt\] This gives us \(a^{3/2} \propto t\), or \(a \propto t^{2/3}\).

(c) Radiation dominated universe. Let us repeat the above exercise for a universe filled with radiation only. For radiation, \(P = \frac{1}{3}\rho c^2\) and \(\Lambda = 0\). Again, use the third Friedmann equation to see how the density of the radiation changes as a function of scale factor. 

We repeat the process above, except we replace P = 0 with \(P = \frac{1}{3}\rho c^2\). Friedmann's third equation becomes: \[\dot{\rho}c^2 = -3(\frac{\dot{a}}{a})(\frac{4}{3}\rho c^2)\] This simplifies to: \[\frac{\dot{\rho}}{\rho} = -4 \frac{\dot{a}}{a}\] Once again we substitute for \(\dot{\rho}\) and \(\dot{a}\) and integrate to get \[\rho = \frac{a_0 \rho_0}{a^4}\] \[\rho \propto a^{-4}\]
Again using the relation between \(\rho\) and a and the first Friedmann equation to show that \(a(t) \propto t^{1/2}\).

Plugging this relation into the first Friedmann equation, we get that \[\dot{a}^2 \propto \rho a^2\] \[\dot{a}^2 \propto a^{-2}\] \[\dot{a} \propto a^{-1}\] Substitute \(\dot{a} = \frac{da}{da}\): \[\frac{da}{dt} \propto a^{-1}\] \[ \int a da \propto \int dt\] \[a^2 \propto t\] \[a \propto t^{1/2}\]

(d) Cosmological constant/dark energy dominated universe. Imagine a universe dominated by the cosmological-constant-like term. Namely in the Friedmann equation, we can set \(\rho = 0\) and P = 0 and only keep \(\Lambda\) nonzero. 
As a digression, notice that we said "cosmological-constant-like" term. This is because the effect of the cosmological constant may be mimicked by a special content of the universe which has a negative pressure \(P = \rho c^2\). Check that the effect of this content on the right-hand-side of the third Friedmann equation is exactly like that of the cosmological constant. To be general we call this content the dark energy. How does the energy density of the dark energy change in time? 

Show that the scale factor of the cosmological-constant-dominated universe expands exponentially in time. What is the Hubble parameter of this universe?

This is the same process again as in (c) and (d), with P = 0, \(\rho\) = 0, \(\Lambda \neq 0\), which yields \[a \propto e^t\] The energy density of the dark energy decreases exponentially in time. The Hubble parameter is thus: \[H(t) = \frac{\dot{a}}{a} = \frac{e^t}{e^t} = 1\]
(e) Suppose the energy density of a universe at its very early time is dominated by half matter and half radiation. (This is in fact the case for our universe 13.7 billion years ago and only 60 thousand years after the Big Bang). As the universe keeps expanding, which content, radiation or matter, will become the dominant component? Why?

The scale factor for matter (\(a \propto t^{2/3}\)) is larger than the scale factor for radiation (\(a \propto t^{1/2}\)). As a result, as the universe keeps expanding, matter will become the dominant component.


(f) Suppose the energy density of a universe is dominated by similar amount of matter and dark energy. As the universe keeps expanding, which content, matter or the dark energy, will become the dominant component? Why? What is the fate of our universe?

Similar to above, the scale factor for dark energy (\(a \propto e^t\)) grows at a faster rate than the scale factor for matter (\(a \propto t^{2/3}\)), and thus as the universe keeps expanding, dark energy will become the dominant component.

Blog Post 28, Worksheet 9.1, Problem 1

1. A Matter-only Model of the Universe in Newtonian approach

In this exercise, we will derive the first and second Friedmann equations of a homogeneous, isotropic and matter-only universe. We use the Newtonian approach. 

Consider a universe filled with matter which has a mass density \(\rho (t)\). Note that as the universe expands or contracts, the density of the matter changes with time, which is why it is a function of time t.

Now consider a mass shell of radius R within this universe. The total mass of the matter enclosed by this shell is M. In the case we consider (homogeneous and isotropic universe), there is no shell crossing, so M is a constant. 

(a) What is the acceleration of this shell?

We know from Newton's second law that the net force acting on an object is equal to its mass times it acceleration: \[F = m\dot{v}\] We know that the force is coming from gravity, so we have: \[F = -G\frac{Mm}{R^2} = m\dot{v}\] The small ms (mass of a particle) cancel, and we are left with: \[\dot{v} = \frac{-GM}{R^2}\]
(b) To derive an energy equation, it is a common trick to multiply both sides of your acceleration equation by v. Turn your velocity, v, into dR/dt, cancel dt, and integrate both sides of your equation. This is an indefinite integral, so you will have a constant of integration; combine these integration constants and call their sum C. Convince yourself the equation you've written down has units of energy per unit mass. 

From above, we have: \[\dot{v} = \frac{-GM}{R^2}\] Multiplying both sides by v, we get: \[\dot{v}v = \frac{-GM}{R^2}v\] We know that \(\dot{v} = \frac{dv}{dt}\) and \(v = \frac{dR}{dt}\). Substituting these in, we get: \[\frac{dV}{dt}v = \frac{-GM}{R^2}\frac{dR}{dt}\] \(\frac{1}{dt}\) cancels on both sides, and we are left with: \[v dv = \frac{-GM}{R^2} dR\] Now we can integrate both sides: \[\int v dv = \int \frac{-GM}{R^2} dR\] \[\frac{v^2}{2} = \frac{GM}{R} + C\] We know that \(v = \dot{R}^2\), so we can rewrite our equation as: \[\frac{1}{2}\dot{R}^2 - \frac{GM}{R} = C\]
(c) Express the total mass M using the mass density, and plug it into the above equation. Rearrange your equation to give an expression for \((\frac{\dot{R}}{R})^2\), where \(\dot{R}\) is equal to \(\frac{dR}{dt}\).

We know that the total mass, M, is equal to \(\rho V\), where V is the total volume. As we are considering a sphere, \(V = \frac{4}{3}πR^3\). Plug this into our equation for C, and we get:\[ \frac{1}{2} \dot{R}^2 - \frac{G(\rho\frac{4}{3}πR^3)}{R} = C\] \[\frac{1}{2}\dot{R}^2- \frac{4G\rhoπR^2}{3} = C\] Rearranging the equation, we get: \[ (\frac{\dot{R}}{R})^2 = \frac{8G\rhoπ}{3} + \frac{2C}{R^2}\]
(d) R is the physical radius of the sphere. It is often convenient to express R as R = a(t)r, where r is the comoving radius of the sphere. The comoving coordinate for a fixed shell remains constant in time. The time dependence of R is captured by the scale factor a(t). The comoving radius equals to the physical radius at the epoch when a(t) = 1. Rewrite your equation in terms of the comoving radius, r, and the scale factor, a(t). 

Replacing R with a(t)r in our above equation, we get: \[ (\frac{\dot{R}}{a(t)r})^2 = \frac{8G\rhoπ}{3} + \frac{2C}{(a(t)r)^2}\] We know that \(\dot{R} = \frac{dR}{dt} = \frac{d(a(t)r)}{dt} = \dot{a}(t)r\), so we can substitute \(\dot{a}(t)r\) for \(\dot{R}\): \[ (\frac{\dot{a}(t)r}{a(t)r})^2 = \frac{8G\rhoπ}{3} + \frac{2C}{(a(t)r)^2}\]
(e) Rewrite the above expression so that \((\frac{\dot{a}}{a})^2) appears alone on the left side of the equation.

Both \(r^2\)s cancel on the left side of the equation, leaving us with: \[ (\frac{\dot{a}}{a})^2 = \frac{8G\rhoπ}{3} + \frac{2C}{(a(t)r)^2}\]

(f) Derive the first Friedmann Equation. From the previous worksheet, we know that \(H(t) = \frac{\dot{a}}{a}\). Plugging this relation into your above result and identifying the constant \(2C/r^2 = -kc^2\) where k is the "curvature" parameter, you will get your first Friedmann equation. The Friedmann equation tells us about how the shell expands or contracts; in other words, it tells us about the Hubble expansion (or contraction) rate of the universe. 
Plugging in  \(H(t) = \frac{\dot{a}}{a}\) into our above equation, we get: \[ (H(t))^2 = \frac{8G\rhoπ}{3} + \frac{2C}{a(t)^2r^2}\] The constant \(2C/r^2 = -kc^2\) can be substituted in the term on the far right, leaving us with: \[ (H(t))^2 = \frac{8G\rhoπ}{3} + \frac{kc^2}{a(t)^2}\]
(g) Derive the second Friedmann Equation: Now express the acceleration of the shell in terms of the density of the universe, and replace R with R = a(t)r. You should see that \(\frac{\ddot{a}}{a} = \frac{-4π}{3}G\rho\), which is known as the second Friedmann equation. 

From part a, we know that \[\dot{v} = \frac{-GM}{R^2}\] We also know that \[M = \rho V = \rho (\frac{4}{3}πR^3)\] Plugging this in for M, we get: \[\dot{v} = -\frac{4}{3}G\rho πR = -\frac{4}{3}G\rho π a(t)r\] We can replace \(\dot{v}\) with \(\frac{d^2R}{dt^2}\). However, we know that \(\frac{d^2R}{dt^2} = \frac{d^2(a(t)r)}{dt^2} = \ddot{a}(t)r\). So we're left with: \[\ddot{a}(t)r = -\frac{4}{3}G\rho π a(t)r\]An r cancels on both sides, and we divide both sides by a(t) to get: \[\frac{\ddot{a}}{a} = \frac{-4Gπ\rho}{3}\]

Blog Post 27, Worksheet 8.1, Problem 2

2. In 1929, astronomer Edwin Hubble discovered that almost all distant galaxies exhibit a positive redshift. Furthermore, it appeared that the farther the galaxy, the larger its redshift. Here we will rediscover Hubble's Law using modern spectroscopic data and supernovae Ia as a our distance indicator to these galaxies. The data we will use come from the Sloan Digital Sky Survey (SDSS), a project that aims to comprehensively map the universe. You can access the relevant data products for this exercise at http://goo.gl/fmIvqc

(a) Below is a list of supernovae observed between 2004 and 2007 and their positions in RA and Dec. You can find the images and spectra of their host galaxies by entering their coordinates in the respective fields. Explore the functions available, including magnifying the image, reading off the photometric measurements (magnitudes in wavebands u, g, r, i, z) of your selected object, and using the 'Explore' button to access more quantitative measurements for these objects. In particular, familiarize yourself with the 'interactive spectrum' feature. 

(b) One of the features for determining distances to Type Ia supernovae is its peak absolute magnitude. You explored the peak bolometric luminosities of SN Ia's in Worksheet 7.1. The peak V-band magnitude for SN Ia's is about -19.3. Use the apparent peak magnitude given in the table above to calculate the distance of these supernovae in unit of Mpc. 

The distance modulus gives the relationship between apparent magnitude (m), absolute magnitude (M), and distance (d): \[m = M = 5\text{log}(d) - 5\] \[d = 10^{\frac{m-M+5}{5}}\] Using an absolute magnitude of M = -19.3 and the apparent magnitudes given, we get the following distances:

Supernova	Right Ascension (RA)	Declination (dec)	extinction-corrected V-magnitude	Distance (Mpc)
SN 2004hu	18.76	0.263	17.28	207
SN 2005gi	13.97	0.505	17.38	216.7
SN 2006rz	56.528	0.39	16.29	131.2
SN 2007hx	31.614	-.899	18.32	334.1

(c) We can use the absorption or emission lines of the host galaxy to find their redshifts which, as you found in Question 1), roughly equals the recessional velocity as a fraction of the speed of light. To measure the redshift to each host galaxy, click on 'Explore' and then on the link 'Interactive Spectrum.' Uncheck the boxes Best Fit and Mark Emission Lines. Zoom in on the absorption line labeled \(H \alpha\), and move your mouse over to the center of the line to read the observed wavelength in Angstroms. The \(H\alpha\) has a rest (i.e. emitted) wavelength of 6563.0 Angstroms. Calculate the redshift, and then derive the radial velocity in kilometers per second, using the relation z = v/c. How close does your redshift measurements compare to the one SDSS reports in the table under the Interactive Spectrum link? Repeat for all galaxies. 

We can find the redshift for each galaxy by using the redshift formula, \[z = \frac{\lambda_{\text{observed}}-\lambda_{\text{emitted}}}{\lambda_{\text{emitted}}}\] We know that \(\lambda_{\text{emitted}}\) for each supernova should be 6563.0 Angstroms, so to find the redshift, we plug in \(\lambda_{\text{emitted}}\) = 6563.0 Angstroms, and \(\lambda_{\text{observed}}\) = the observed wavelength in Angstroms as shown in the table.

To find the velocity, we multiply the redshift by the speed of light: \(v = cz\).

Supernova	Observed Wavelength	Redshift	Reported redshift	Velocity (km/s)
SN 2004hu	6878.599	0.04809	0.047806	14427
SN 2005gi	6897.632	0.05099	0.05078596994	15297
SN 2006rz	6767.06	0.03109	0.0309035	9327
SN 2007hx	7085.981	0.07967	0.0794432	23901

(d) Make a plot of your findings, with distance on the x-axis and velocity on the y-axis. Report the slope of your line in appropriate units. This is the Hubble Constant, \(H_0\). 

\[H_0 = 72.1 \frac{\text{km/s}}{\text{Mpc}}\]

(e) Write an equation for this line in the form of v = __ D, where v is an object's recessional velocity and D is the distance to that object. Express your Hubble Constant in terms of units km/s/Mpc. Congratulations, you have arrived at Hubble's Law!

This question just asks us to plug in the Hubble Constant we found above as the slope for this line: \[v = H_0 D = (72.1 \frac{\text{km/s}}{\text{Mpc}})D\]

Blog Post 26, Worksheet 8.1, Problem 1

1. Before we dive into the Hubble Flow, let's do a thought experiment. Pretend that there is an infinitely long series of balls sitting in a row. Imagine that during a time interval \(\Delta t\) the space between each ball increases by \(\Delta x\). 

(a) Look at the shaded ball, Ball C, in the figure above. Imagine that Ball C is sitting still (so we are in the reference frame of Ball C). What is the distance to Ball D after time \(\Delta t\)? What about Ball B?

After some \(\Delta x\), the distance between Ball D and Ball C will be \(\Delta x\), and the distance between Ball B and Ball C will also be \(\Delta x\), but in the opposite direction.

(b) What are the distances from Ball C to Ball A and Ball E?

As the distance between Ball C and Ball A was originally twice the distance from Ball C to Ball B, the distance between Ball C and Ball A is now \(2\Delta x\). The same is true for Ball C and Ball E.

(c) Write a general expression for the distance to a ball N balls away from Ball C after time \(\Delta t\). Interpret your finding.

After some time \(\Delta t\), a ball originally N balls away from Ball C will be \(N \Delta x\) away from Ball C. This means physically that the space between any two balls x and y a distance N apart will also be \(N \Delta x\) away from each other.

(d) Write the velocity of a ball N balls away from Ball C during \(\Delta t\). Interpret your finding. 

The velocity of a ball N balls away from Ball C during this time interval is simply given by \(\frac{\text{distance}}{\text{time}}\): \[v = \frac{N\Delta x}{\Delta t}\] As the distance between Ball C and some other ball increases, the velocity of that ball increases. This means that in a uniformly expanding universe, objects farther away from an object appear to be moving at a greater velocity.

Blog Post 25, Worksheet 7.2, Problem 2

2. There is actually a hard upper limit to the luminosity of this system - and to the luminosity of any accreting compact object. Consider that the photons being emitted in this scenario will interact with the surrounding material (which has yet to accrete onto the black hole). These photons will undergo Thomson scattering off of electrons in this material. In detail, the electric field of the incident lightwave (i.e., the photon) will accelerate an electron, causing it to then re-emit radiation. The photons are therefore able to transfer some of their momentum to the infalling gas. 

The energy flux of these photons at a distance r from the black hole is \[F = \frac{L}{4πr^2}.\] Then, recall that the momentum of a photon of energy E is simple p = E/c Therefore, the momentum flux at r from the black hole is \(\frac{L}{4πcr^2}\). 

Finally, the rate of momentum transfer to the surrounding electrons (or the force due to photons (\(f_{\text{rad}}\))) is modulated by the Thomson cross section, \(\sigma_t\) = 6.6524 * \(10^{-25} \text{ cm}^2\) (i.e., the effective area of an electron interacting with a photon): \[f_{\text{rad}} = \sigma_t \frac{L}{4πcr^2}\]

(a) When this force from radiation pressure exceeds the force of gravity, accretion is halted and all the gas is blown away. For a black hole mass \(M_{\text{BH}}\), derive the maximum possible luminosity due to accretion. This is called the Eddington Luminosity. 

We know that the Eddington Luminosity will occur at the point where the force of gravity is equal to the force from radiation pressure, or when: \[\frac{GM_{\text{BH}}M_{\text{proton}}}{r^2} = \sigma_t \frac{L}{4πcr^2}\] If we solve for L, we get: \[L = L_{\text{Edd}} = \frac{4πcGM_{\text{BH}}M_{\text{proton}}}{\sigma_t}\]
(b) If the SMBH in Andromeda were accreting at 20% of its Eddington luminosity, how bright would it be? How does this value compare with the \(L_{\text{disk}}\) you calculated in 1(e)?

The Eddington luminosity for Andromeda's SMBH is given by \[L_{\text{Edd}} = \frac{4πcGM_{\text{BH}}M_{\text{proton}}}{\sigma_t}\] We know that \(M_\text{BH} = 2 * 10^8 M_\odot\), and \(M_{\text{proton}} = 1.7 * 10^{-24} g\). We can plug in these values to find 20% of the Eddington luminosity: \[ 0.2L_{\text{Edd}} = 0.2(\frac{4π(3 * 10^{10} \text{ cm/s})(G)(2 * 10^8 M_\odot)(2 * 10^{33} \text{g}/M_\odot)(1.7 * 10^{-24} g)}{6.6524 * 10^{-25} \text{ cm}^2})\] This gives us that the luminosity of the SMBH is \(5.2 * 10^{45}\) erg/s, which is almost exactly what we found \(L_\text{disk}\) to be in problem 1(e).

Blog Post 24, Worksheet 7.2, Problem 1

1. Recall that on Worksheet 6.1, you found that many galaxies host a central supermassive black hole and derived the mass of the super-massive black hole (SMBH) in Andromeda (M31). In galaxies that contain large amounts of gas rotating along with the stars in their disks, some of the gas will end up surrounding the black hole. As the galls falls inward, centrifugal forces will cause it to settle into a much smaller, denser disk which rotates around the central black hole, known as an "accretion disk."

The orbital period of the gas in the accretion disk changes with radius (just like the orbital periods of the planets). This means that adjacent packets of gas in the disk will rub against each other, generating frictional heat which gets released as radiation. The energy loss results in the material moving closer to the black hole, eventually falling onto it (hence, accretion).

Let's figure out how hot and bright this material ends up being around a black hole of mass \(M_{BH}\).

(a) First, imagine a gas packet of mass dM in the accretion disk. The gas packet falls from a radius r + dr to a radius r during the accretion process, losing potential energy. With the help of the Virial Theorem, write down the thermal energy the packet gains as a result (dE). Then, simply divide this dE by dt to express the luminosity of the packet, dL, in terms of of \(M_{BH}\), r, and the mass accretion rate dM/dt = \(\dot{M}\).

We know that potential energy of this system comes from gravitational potential energy. Thus, the initial potential energy between the mass dM and the total mass M is given by \[U_i = \frac{GM_{\text{BH}}dM}{r+dr}\] and the final potential energy is given by \[U_f = \frac{GM_{\text{BH}}dM}{r}\] The change in this potential energy is thus just given by: \[\frac{d}{dr}(U) = -\frac{GM_{\text{BH}}dM}{r^2}\] From the Virial Theorem, we know that \(K = -\frac{1}{2}U\), so: \[K = --\frac{1}{2}U = \frac{GM_{\text{BH}}dM}{2r^2} = dE\] Now, we divide each side by dt to get: \[\frac{dE}{dt} = dL = \frac{GM_{\text{BH}}dM}{2r^2dt} = \frac{\dot{M}GM_{\text{BH}}}{2r^2}\]
(b) Now let's assume the disk gas radiates like a blackbody at the same radius where the potential energy is released. Using the Stefan-Boltzmann law, write down an expression for the luminosity from a given annulus in the disk (between r and r + dr) in terms of the temperature in that annulus, T(r). 

The Stefan-Boltzmann law says that the luminosity of a blackbody is given by \(L = A\sigma T^4\), where A is the area of the object. Because we are looking at an annulus ring, we know that the area of the annulus is given by \(A = π((r+dr)^2-r^2)\). However, because we will be integrating over the annulus, each small ring of the annulus at some radius r will have an area of 2πr. Thus, the luminosity from a given annulus in the disk is given by: \[L = 2π\sigma T(r)^4 r\]
(c) Now, set dL from (a) equal to the expression you just derived in (b). Solve for T(r), and express in terms of the mass accretion rate, \(\dot{M}\). 
\[dL = 2π\sigma T(r)^4 r\] Now solve for T(r): \[T(r) = (\frac{dL}{2π\sigma r})^{1/4}\] From part a, we know that \[dL = \frac{\dot{M}GM_{\text{BH}}}{2r^2},\] which we can plug in to get: \[T(r) = (\frac{ \frac{\dot{M}GM_{\text{BH}}}{2r^2}}{2π\sigma r})^{1/4}\] \[T(r) = (\frac{GM_{\text{BH}}\dot{M}}{4π\sigma r^3})^{1/4}\]
(d) Finally, integrate your expression for dL over dr to find the total disk luminosity, \(L_{\text{disk}}\). Let's say here that the disk extends from inner radius \(r_{\text{in}}\) to outer radius \(r_{\text{out}}\). If \(r_{\text{out}}\) >> \(r_{\text{in}}\), what is \(L_{\text{disk}}\)?

From before, we know that \(dL = \frac{\dot{M}GM_{\text{BH}}}{2r^2}\). Now we integrate:  \[\int dL = \int_{r_{\text{out}}}^{r_{\text{in}}} \frac{\dot{M}GM_{\text{BH}}}{2r^2}dr\] \[L = \frac{GM_{\text{BH}}\dot{M}}{2} \int_{r_{\text{out}}}^{r_{\text{in}}} r^{-2}dr\] \[L = \frac{GM_{\text{BH}}\dot{M}}{2} [-r^{-1}]|_{r_{\text{out}}}^{r_{\text{in}}}\] \[L = \frac{GM_{\text{BH}}\dot{M}}{2r_{\text{in}}}\]
(e) Astronomers find it useful to express the \(L_{\text{disk}}\) you found in 1d) as a fraction, \(\eta\), of the total luminosity that would be released if the entire rest mass of the disk were converted to energy, \(\dot{M}c^2\): \[L_{\text{disk}} = \eta \dot{M} c^2\] Show that \(\eta\) is given by \[\eta = \frac{1}{2} \frac{GM_{\text{BH}}}{c^2r_{\text{in}}}\] We think of this as the radiative efficiency of the disk. And this efficiency apparently has a strong dependence on the inner disk radius (\(r_{\text{in}}\)). 

If \(L_{\text{disk}} = \eta\dot{M}c^2\), then \[\eta = \frac{L_{\text{disk}}}{\dot{M}c^2} = \frac{\frac{GM_{\text{BH}}\dot{M}}{2r_{\text{in}}}}{\dot{M}c^2} = \frac{1}{2}\frac{GM_{\text{BH}}}{c^2r_{\text{in}}}\]
(f) Accretion disks around massive black holes tend to have low efficiencies of \(\eta \approx 0.1\), and can accrete up to \(\dot{M} = 1 \:M_\odot \text{ yr}^{-1}\). What is the resulting disk luminosity? How hot would an accretion disk around the SMBH in Andromeda be?

We are given that \(\eta\)  = 0.1, and \(\dot{M} = 1 M_\odot \text{ yr}^{-1}\) . We plug this into the formula \(L_{\text{disk}} = \eta\dot{M}c^2\) to get:
\[L_{\text{disk}} = (0.1)(\frac{1 M_\odot}{1 yr})(\frac{1 yr}{3.2*10^7 \text{ s}})(\frac{2*10^{33} \text{ g}}{1 M_\odot})(3*10^10 \text{ cm/s})^2\] \[L_{\text{disk}} = 5.7 * 10^{45} \text{ erg/s}\]

From the Stefan-Boltzmann Law, we know that \(L = A \sigma T^4\), and therefore \(T = (\frac{L}{A\sigma})^{1/4}\). We know that the luminosity of the disk is given by \(5.7 * 10^{45} \text{ erg/s}\), and the radius of the black hole is given by \(R = \frac{2GM_{\text{BH}}}{c^2}\).  We therefore have \[T = (\frac{L_{\text{disk}}}{\frac{4}{3}π(\frac{2GM_{\text{BH}}}{c^2})^3\sigma})^{1/4}\] \[T = (\frac{5.7 * 10^45 \text{ erg/s}}{\frac{4}{3}π(\frac{2G(2 * 10^8 M_\odot)(2 * 10^{33} \frac{g}{M_\odot}}{(3*10^{10})^2})^3(5.7*10^{-5} \text{erg cm}^{-2}\text{ K}^{-4}\text{ s}^{-1})})^{1/4}\] \[T = 100 \:K\]

Blog Post 23, Astronomy and Religion

The intersection of astronomy and religion is one that goes back millennia – for as long as humans have been around, they have looked to the stars for the answers to our species’ most profound questions – Why are we here? How are we here? What are we looking at? This fall, I am taking a class called Humans in Space: Past, Present, and Future (History of Science 181), taught by Professor Matthew Hersch. In this course, we have covered many topics relating to the history of spaceflight and our views of the universe. A reoccurring theme that has appeared throughout this class, and through history, is the constant push-and-pull between astronomers and the religious views of the times. While this topic could spur an endless blog post covering many cultures as they have developed through time, I just want to focus on a few points in history that I found really interesting.

For many centuries, the universe was thought to be centered on the Earth, with all objects including the Sun, Moon, planets, and stars orbiting around the Earth. Known as the geocentric model or the Ptolemaic system, this model was accepted in many ancient civilizations, including the ancient Greeks. Prior to the Renaissance, the European Church had adopted this model, as developed by Aristotle and Ptolemy, as part of their own doctrine.

Source: https://upload.wikimedia.org/wikipedia/en/2/24/Ptolemaic-geocentric-model.jpg

Other theories were brewing, however. Although heliocentric theory had been first developed much earlier by Aristarchus (ca. 270 B.C.), Nicolaus Copernicus truly took the theory under his wing, publishing his book De revolutionibus orbium coelestium (On the Revolutions of the Celestial Spheres) shortly before his death in 1543, and forcing the world to question their views of the universe. Geocentric theory had agreed very nicely with Church doctrine, placing the Earth at the center of the universe as suggested in the Bible. Heliocentric theory called the Bible and those in power because of it into question – as you can imagine, the Church wasn’t all too happy with this new world view.

This unhappiness continued to worsen, and by the time of Giordano Bruno, they fully rejected the Copernican model. Bruno upset the Church even further by suggesting a plurality of worlds – not only a heliocentric universe, but a universe in which all of the stars we look up at are in fact their own universe. While this idea may seem harmless at first, Bruno framed it in a way that severely threatened the Church’s power – if the heavens to us on Earth are in fact individual universes, then the Earth must be the heavens to those other universes. And if the Earth is heaven to other universes, how can the Church have the power to determine who goes to heaven, if our world is already heaven to another world? Bruno went so far as to suggest that the Pope can’t in fact have any power, and was literally burned at the stake in 1600 for his ideas. Think about that – burned at the stake, for attempting to advance our scientific world view.

Eventually the Church did mold its accepted model closer to what we now know to be correct, but it took some steps to get there. Tycho Brahe (1546-1601) developed a model (the Tychonic system) in which the Sun continues to orbit the Earth, but all planets orbit around the Sun. This model agreed more closely with observation, yet also satisfied the Church's view that the Earth is at the center of the universe. For a while, the Church adopted the Tychonic system as their official accepted view of the universe. 

Source: https://upload.wikimedia.org/wikipedia/commons/a/a6/Tychonian_system.svg

Eventually the world came to accept heliocentrism as the correct model for our solar system, and to believe that there is no actual center of the universe, according to Einstein's principle of relativity. Yet looking at this history makes one really marvel at the path our world view took to develop - it makes me glad to be alive at a time and in a place where science is in a much more peaceful and separate relationship with religion, yet also raises a question: what do we believe now based on our observations, that is in actuality completely and utterly wrong? What will our next big shift in world view occur?