EB Lab Report

1. Introduction

This lab report details the data and findings from the observation of the members of the double-lined spectroscopic binary NSVS01031772. We refer to the paper NSVS01031772: A New 0.50+0.54 M⊙ Detached Eclipsing Binary by López-Morales et al., 2006.

Objective

The purpose of the lab is to:

• determine masses, radii, and separation between the members of the double-lined spectroscopic binary NSVS01031772 using both provided radial velocity (RV) data and the light curve data that we collect using the Clay Telescope
• familiarize ourselves with the Clay Telescopes, and telescopes in general, by using it to perform observations, learning how to use tools such as TCS, MaximDL, and The Sky.

Background

A binary star system is comprised of two stars orbiting a common point. In particular, a spectroscopic binary star system is one that is too far away for a telescope to resolve the two stars individually. Instead, we can measure the periodic change in the Doppler shifts of the spectral lines, as the two individual star revolve about their common center of mass. These Doppler shifts are then used to construct a radial velocity (RV) plot.

Furthermore, the changes in brightness over time can be used to construct a light curve. As one star crosses in front of the other, it covers some of the light from the other star, resulting in a change in brightness. This is illustrated in Figure 1:

Figure 1: An eclipsing binary star with light curve.
Image credits: https://www.physics.mun.ca/~jjerrett/binary/binary.html

Since the stars are not individually resolvable, we just observe a single point that changes in brightness over time. NSVS01031772 is a low-mass binary system, which are difficult to detect because of their low surface brightness. As a result, there is a dearth of observational data to inform the mass-radius relationship for stars less massive than 0.6 $M_{\odot}$. However, this information is important in constraining evolutionary models of stars, and also inform the masses and radii of planets that orbit these stars.

2. Methods

Radial Velocity Plot

A radial velocity (RV) is a plot of radial (line-of-sight) velocity versus time, constructed using the Doppler shifts of the spectral lines. López-Morales et al. included an RV plot in their paper, shown in Figure 2.

Figure 2: Radial velocity curves of NSVS0103.
Image credits: López-Morales et al., 2006

Positive radial velocity means that the star is going away from the observer, while negative radial velocity means that the star is coming towards the observer. The horizontal dashed line indicates the velocity of the center of mass of the system, which is about $19 \text{ km/s}$. This velocity arises from the relative motion between us and the center of mass of the binary star system.

Meanwhile, we can read the maximum radial velocities to determine the mass ratio of the system; the greater mass travels at lower velocity, while the lesser mass travels at a higher velocity. As shown in Figure 3, the movement of the stars can be inferred from the RV plot.

Figure 3: Relative motion of the stars

The solid line represents the radial velocity of Star 1, and the dotted line represents Star 2. From the plot, Star 1 has a maximum radial velocity of $-125 \pm 1 \text{ km/s}$, and Star 2 has a maximum radial velocity of $176 \pm 1 \text{ km/s}$. Then, we normalize these values by subtracting the velocity of the center of mass, 19 km/s:

$$v_1 = -125 \text{ km/s} - 19 \text{ km/s} = 144  \pm 1 \text{ km/s}$$

$$v_2 = 175 \text{ km/s} - 19 \text{ km/s} = 156 \pm 1 \text{ km/s}$$

From these data, we determine that Star 1 is more massive than Star 2. We can determine the mass ratio of the two stars using the center of mass equation:

$$M_1 a_1 = M_2 a_2$$

Where $M_1$, $M_2$ are the masses of Stars 1 and 2, respectively, and $a_1$, $a_2$ are the distances from the center of mass to Stars 1 and 2, respectively. We can relate distance $a$ to radial velocity $v$ using the equation for orbital period:

$$v = \frac{2\pi a}{P}$$

Now, we can rewrite the center of mass equation as follows:

$$M_1 v_1 = M_2 v_2$$

Rearrange and solve:

$$\frac{M_1}{M_2} = \frac{v_2}{v_1} = \frac{156 \text{ km/s}}{144 \text{ km/s}} \approx 1.08$$

Therefore, our mass ratio $\frac{M_1}{M_2}$ is about 1.08.

Light Curve

A light curve is generated by plotting the brightness versus time of an astronomical body. We can use them to determine the type of body that we are observing; for example, a binary star system has a periodic light curve, as shown in Figure 1.

We collected optical data from NSVS0103 using the Clay Telescope to generate our own light curves, plotting light intensity versus time. From Worksheet 5, Problem 2e, we found that luminosity $L(T) = 4\sigma\pi R^2 T^4$, which indicates that the brighter star is the one with a greater radius; conversely, the smaller star will be dimmer. When both stars are next to each other from our perspective, we get the light from both of them. However, when the smaller star (Star 2) transits in front of Star 1 (primary transit), it covers some of the light from Star 1, so we get reduced brightness, leading to a major dip. When Star 2 transits behind Star 1 (secondary transit), only the light from Star 1 is visible, so we again get reduced brightness, but only a minor dip.

These dips are periodic, meaning we can use them to determine the orbital period of the binary star system. In addition, we can use light curve to determine the depths of the dips, their durations, and the time of ingress and egress of the eclipsing star. We then use this information to calculate the masses and radii of the binary stars.

Mass

First, we use Kepler's third law to determine the masses of the stars:

$$P^2 = \frac{4\pi^2 a^3}{G(M_1 + M_2)}$$

We also know that the semimajor axis is the sum of $a_1$ and $a_2$:

$$a = a_1 + a_2$$

Furthermore, we can rewrite $M_2$ in terms of $M_1$ using what we know about their mass ratio:

$$M_2 = \frac{v_1}{v_2}M_1$$

Plug these in:

$$P^2 = \frac{4\pi^2 (a_1 + a_2)^3}{G(M_1 + \frac{v_1}{v_2}M_1)}$$

Next, we can use the equation for orbital period, assuming a circular orbit:

$$Pv_1 = 2\pi a_1$$

$$Pv_2 = 2\pi a_2$$

Then, we can express $a_1 + a_2$ as:

$$P(v_1 + v_2) = 2\pi (a_1 + a_2)$$

$$a_1 + a_2 = \frac{P}{2\pi}(v_1 + v_2)$$

Substitute this for $a$:

$$P^2 = \frac{4\pi^2 (\frac{P}{2\pi}(v_1 + v_2))^3}{G(M_1 + \frac{v_1}{v_2}M_1)}$$

Solve for $M_1$:

$$P^2 = \frac{4\pi^2 P^3(v_1 + v_2)^3}{8\pi^3 G(M_1 + \frac{v_1}{v_2}M_1)}$$

$$P^2 = \frac{P^3(v_1 + v_2)^3 v_2}{2\pi GM_1(v_2 + v_1)}$$

$$Pv_2 (v_1 + v_2)^2 = 2\pi G M_1$$

$$M_1 = \frac{P(v_1 + v_2)^2 v_2}{2\pi G}$$

Similarly, we can solve for $M_2$:

$$M_2 = \frac{P(v_1 + v_2)^2 v_1}{2\pi G}$$

Radius

To solve for the radii of the stars, we need the transit time and the transit depths.

The transit time is derived by taking a reference frame with respect to the larger star. This is illustrated in Figure 4.

Figure 4: Transit time of Star 2

If we hold Star 1, then Star 2 travels at a relative velocity of $v_1 + v_2$. A full transit takes a distance of $2R_1 + 2R_2$, which means that transit time can be expressed as:

$$t_{transit} = \frac{2R_1 + 2R_2}{v_1 + v_2}$$

According to the CfA's Transit Light Curve Tutorial, transit depth is given by:

$$\delta = (\frac{R_2}{R_2})^2$$

However, we need to account for the fact that the transiting Star 2 emits its own light, as the given equation assumes the transit of a planet with negligible luminosity. To compensate for this, we can subtract the luminosity of Star 2 from the light curve, which gives us a new equation:

$$(\frac{R_2}{R_1})^2 = \frac{\delta_1}{1- \delta_2}$$

We can rearrange this to isolate $R_2$ and express it in terms of $R_2$:

$$R_2 = R_1 \sqrt{\frac{\delta_1}{1-\delta_2}}$$

Then, we can rearrange our equation for transit time $t_transit$ and plug in $R_2$ to solve for $R_1$:

$$\frac{1}{2}t_{transit}(v_1 + v_2) = (R_1 + R_2)$$

$$\frac{1}{2}t_{transit}(v_1 + v_2) = (R_1 + R_1 \sqrt{\frac{\delta_1}{1-\delta_2}})$$

$$\frac{1}{2}t_{transit}(v_1 + v_2) = R_1(1 + \sqrt{\frac{\delta_1}{1-\delta_2}})$$

$$R_1= \frac{t_{transit}(v_1 + v_2)}{2(1 + \sqrt{\frac{\delta_1}{1-\delta_2}})}$$

Similarly, we can solve for $R_2$:

$$R_2= \frac{t_{transit}(v_1 + v_2)}{2(1 + \sqrt{\frac{1-\delta_2}{\delta_1}})}$$

Observations

We performed our observations using the Clay Telescope on top of the Science Center. We used TCS (Telescope Control System) to control the movements of the telescope, allowing us to slew to a specific RA and DEC. In addition, we also used The Sky program control the telescope, allowing us to locate an astronomical body from a large chart, and automatically slew the telescope to point at the desired object. The telescope automatically tracks an object, accounting for the Earth's rotation. Although the movement of the telescope is very fine, the dome moves in much coarser motions, and can get loud and frightening.

The observations were performed at night, split into different groups. Observations were performed on March 26, March 27, March 28, and April 5 of 2019. My group observed on March 26, around 12 AM to 2 AM. The sky was clear and the temperature was about 40 degrees Fahrenheit.

First we slewed the telescope to our field: RA, DEC = 13:45:35 +79:23:48. We used the R-band filter on the CCD. To make sure we were in the right place, we used a finder chart, shown in Figure 5.

Figure 5: Finder chart
Source: https://sites.fas.harvard.edu/~astrolab/finding_chart.pdf

Then, to make sure we could combine the data from various lab sections, we made sure that all images contained the same reference stars, shown in Figure 6. The system we are concerned with, NSVS0103, is labeled "Obj1".

Figure 6: Object field. Obj1 is the object we are observing, while the others are reference stars.
Source: https://sites.fas.harvard.edu/~astrolab/object_field.png

To determine the right exposure time, we use a bit of trial and error. We take a test exposure, and then view the result in MaximDL. We hover of the center of the object star and note the counts value. Since we want to keep our counts ideally 20,000-30,000 counts, we can increase the exposure time if our counts value was too low, or decrease it if it was too high. We used an exposure time of about 2 minutes.

After setting everything up, we leave the telescope to take a sequence of images, allowing us to go back into the lab and monitor the images as they came in.

3. Analysis

Note that we need to flat field our images, which removes noise from defects such as dust. This was done for us by the course staff. Flat frames were captured just after sunset, with the telescope skewed slightly between each frame, so that the frames can be averaged, making sure that consistent noise (like dust) would be removed. Figure 7 show an image taken from the telescope, pre-flat field. Figure 8 show the result after flat fielding.

Figure 8: Top: Before flat field. Bottom: After flat field. Much cleaner!
Source: Astro 16 Canvas

We used MaximDL to perform our photometry, which is the technique of measuring the intensity of electromagnetic radiation of an astronomical body. We use the reference stars to normalize our measurements by making sure they have the same values across our datasets. Each observation group will have different settings for image exposure, so the absolute measured intensity will vary. Instead, we detrend the measurements so that the measured intensity matches known values for the reference stars. Figure 9 shows the completed light curve from all 4 observations.

Figure 8: Final Light Curve

The original data was recorded with time expressed as hours after March 26 UT. To plot the light curve, we needed to "fold" the time based on the period of the binary star system, in order to get the curves to line up. Numerically, this is done by taking the unfolded time modulo period, which folds up the time to "fit" into period-wide chunks. We used a period of 8.83 hours, or 8 hours, 50 minutes, based on the results from the NSVS0103 paper, which reported an orbital period of 0.3681414 days.

The baseline value (flat part of the light curve) appears to be about $1.07 \pm 0.01$. The minimum of the primary transit (rightmost dip) is about 0.57, so the depth is the difference $1.07 - 0.57 = 0.50 \pm 0.02$. The minimum of the secondary transit (leftmost dip) is about 0.63, so the depth is the difference $\delta_2 = 1.07 - 0.63 = 0.44 \pm 0.02$. We need to normalize the transit depths with respect to the baseline, so we compute divide each transit depth by 1.07, giving us $\delta_1 = 0.47 \pm 0.02$ and $\delta_2 = 0.41 \pm 0.02$.

We measured four transit times, corresponding to each observation night: 1.464, 1.368, 1.416, and 1.392 hours. We take the average to get $t_{transit} = 1.41 \text{ hours}$ for an eclipse.

4. Results

Radius

We can solve for the radii $R_1$ and $R_2$, knowing the transit depths $\delta_1 = 0.47$, $\delta_2 = 0.41$. Transit time $t_{transit} = 1.41 \text{ hr}$.

$$R_1 = \frac{t_{transit}(v_1 + v_2)}{2(1 + \sqrt{\frac{\delta_1}{1-\delta_2}})}$$

$$R_1 = \frac{1.41 \text{ hr}\times \frac{3600 \text{ s}}{1 \text{ hr}}(144 \text{ km/s} + 156 \text{ km/s})(\frac{10^5 \text{ cm}}{\text{ km}})}{2(1 + \sqrt{\frac{0.47}{1-0.41}})}$$

$$R_1 = 4.02 \times 10^{10} \text{ cm} = 0.58 R_{\odot}$$

Similarly, we can solve for $R_2$:

$$R_2= \frac{t_{transit}(v_1 + v_2)}{2(1 + \sqrt{\frac{1-\delta_2}{\delta_1}})}$$

$$R_2 = \frac{1.41 \text{ hr}\times \frac{3600 \text{ s}}{1 \text{ hr}}(144 \text{ km/s} + 156 \text{ km/s})(\frac{10^5 \text{ cm}}{\text{ km}})}{2(1 + \sqrt{\frac{1-0.41}{0.47}})}$$

$$R_2 = 3.59 \times 10^{10} \text{ cm} = 0.52 R_{\odot}$$

Mass

Now we can go back and solve for the masses $M_1$ and $M_2$. $v_1 = 144  \pm 1 \text{ km/s}$ and $v_2 = 156 \pm 1 \text{ km/s}$. $P = 8.83 \text{ hours}$.

$$M_1 = \frac{P(v_1 + v_2)^2 v_2}{2\pi G}$$

$$M_1 = \frac{(8.83 \text{ hr} \times \frac{3600 \text{ s}}{1 \text{ hr}}) (144 \text{ km/s} + 156 \text{ km/s})^2 (156 \text{ km/s})(\frac{10^5 \text{ cm}}{\text{ km}})}{2\pi (6.67 \times 10^{-8} \text{ cm}^3 \text{ g}^{-1} \text{ s}^{-2})}$$

$$M_1 = 1.06 \times 10^{33} \text{ g} = 0.53 M_{\odot}$$

Similarly, we can solve for $M_2$:

$$M_2 = \frac{P(v_1 + v_2)^2 v_1}{2\pi G}$$

$$M_2 = \frac{(8.83 \text{ hr} \times \frac{3600 \text{ s}}{1 \text{ hr}}) (144 \text{ km/s} + 156 \text{ km/s})^2 (144 \text{ km/s})(\frac{10^5 \text{ cm}}{\text{ km}})}{2\pi (6.67 \times 10^{-8} \text{ cm}^3 \text{ g}^{-1} \text{ s}^{-2})}$$

$$M_2 = 9.83 \times 10^{32} \text{ g} = 0.49 M_{\odot}$$

Separation between the stars

Separation between the stars is the sum $a_1 + a_2$. We use the equation that relates $a$ to period:

$$a_1 + a_2 = \frac{P}{2\pi}(v_1 + v_2) = a$$

$$a = \frac{8.83 \text{ hr} \times \frac{3600 \text{ s}}{1 \text{ hr}}}{2\pi}(144 \text{ km/s} + 156 \text{ km/s})(\frac{10^5 \text{ cm}}{\text{ km}})$$

$$a = 1.51 \times 10^{11} \text{ cm} = 2.18 R_{\odot}$$

Error

Radius

The transit time values were taken from the spreadsheet, making them relatively accurate. $t_{transit} = 1.41 \pm 0.05 \text{ hours}$

The values of $v_1$ and $v_2$ were drawn from the RV plot, and lines on the computer were used to make sure they were accurate, resulting in $v_1 = 144  \pm 1 \text{ km/s}$ and $v_2 = 156 \pm 1 \text{ km/s}$.

$$\delta v = \sqrt{(1)^2 + (1)^2} = 1.414 \text{ km/s}$$

The transit depths are as follows:
$\delta_1 = 0.47 \pm 0.02$
$\delta_2 = 0.41 \pm 0.02$

To calculate the radius, the transit time is multiplied by the sum of the two velocities, divided by a square root of one ratio divided by another.

$$\delta R_1 = 4.02 \times 10^{10} \text{ cm}\sqrt{(\frac{0.05}{1.41})^2 + (\frac{1.414}{144+156})^2 + \frac{1}{2}(\frac{0.02}{0.47})^2 + \frac{1}{2}(\frac{0.02}{0.41})^2}$$
$$\delta R_1 = 2.23 \times 10^9 \text{ cm}$$

$$\delta R_2 = 3.59 \times 10^{10} \text{ cm}\sqrt{(\frac{0.05}{1.41})^2 + (\frac{1.414}{144+156})^2 + \frac{1}{2}(\frac{0.02}{0.47})^2 + \frac{1}{2}(\frac{0.02}{0.41})^2}$$
$$\delta R_2 = 2.00 \times 10^9 \text{ cm}$$

$$R_1 = 4.02 \times 10^{10} \pm 2.23 \times 10^9 \text{ cm}$$
$$R_2 = 3.59 \times 10^{10} \pm 2.00 \times 10^9 \text{ cm}$$

Mass

Mass is calculated from the product of the period, the square of the sum of the velocities, and a velocity. The error for the period is quite low, because it was confirmed with the result from the paper.

$P = 8.83 \text{ hours} \pm 0.005 \text{ hours}$
$v_1 = 144  \pm 1 \text{ km/s}$
$v_2 = 156 \pm 1 \text{ km/s}$

$$\delta v = \sqrt{(1)^2 + (1)^2} = 1.414 \text{ km/s}$$

$$\delta M_1 = 1.06 \times 10^{33} \text{ g} \times \sqrt{(\frac{0.005}{8.83})^2 + (\frac{1.414}{144+156})^2 + (\frac{1.414}{144+156})^2 + (\frac{1}{156})^2}$$

$$\delta M_1 = 9.82 \times 10^{30} \text{ g}$$

$$\delta M_2 = 9.83 \times 10^{32} \text{ g} \times \sqrt{(\frac{0.005}{8.83})^2 + (\frac{1.414}{144+156})^2 + (\frac{1.414}{144+156})^2 + (\frac{1}{144})^2}$$

$$\delta M_2 = 9.48 \times 10^{30} \text{ g}$$

$$M_1 = 1.06 \times 10^{33} \pm 9.82 \times 10^{30} \text{ g}$$

$$M_2 = 9.83 \times 10^{32} \pm 9.48 \times 10^{30} \text{ g}$$

Separation

The calculation for separation distance is comprised of a multiplication of the period by the sum of the two velocities.

$v_1 = 144  \pm 1 \text{ km/s}$
$v_2 = 156 \pm 1 \text{ km/s}$
$P = 8.83 \text{ hours} \pm 0.005 \text{ hours}$

$$\delta v = \sqrt{(1)^2 + (1)^2} = 1.414 \text{ km/s}$$
$$\delta a = 1.51 \times 10^{11} \text{ cm} \times \sqrt{(\frac{0.005}{8.83})^2 + (\frac{1.414}{144+156})^2} = 7.17 \times 10^8 \text{ cm}$$

$$a = 1.51 \times 10^{11} \pm 7.17 \times 10^8 \text{ cm}$$


Possible sources of error include weather conditions on the various observation days. For example, there is a lot of light pollution in Boston, and clouds in the sky could cause that light to scatter. There could also be variation in the flat fielding, depending on at what point after sunset they were performed, how how much the telescope was slewed in between frames.

In conclusion, the determined values were:

$R_1 = 3.91 \times 10^{10} \text{ cm}$

$R_2 = 3.70 \times 10^{10} \text{ cm}$
$M_1 = 1.06 \times 10^{33} \text{ g}$
$M_2 = 9.83 \times 10^{32} \text{ g}$
$P = 8.83 \text{ hours}$
$a = 1.51 \times 10^{11} \text{ cm}$

These match up reasonably closely to the values from the paper, which lends support to the accuracy of our procedures and calculations.

References:

http://hosting.astro.cornell.edu/academics/courses/astro201/bin_vels.htm

https://www.astronomynotes.com/starprop/s10.htm

http://www.astro.keele.ac.uk/astrolab/manual/week11.pdf

https://imagine.gsfc.nasa.gov/features/yba/CygX1_mass/binary/equation_derive.html

https://imagine.gsfc.nasa.gov/science/toolbox/timing1.html

https://ay16-demery.blogspot.com/2015/04/eclipsing-binary-lab.html

https://vimeo.com/92354083

https://www.cfa.harvard.edu/~avanderb/tutorial/tutorial.html

https://sites.fas.harvard.edu/~astrolab/EB_Ay16.html

Week 10, Post 6

General Post 1:

How does the result of the first imaged Black Hole relate to our galaxy? Write more than a simple statement of fact.

https://www.nytimes.com/2019/04/10/science/black-hole-picture.html
https://eventhorizontelescope.org/

The image of the black hole came from Messier 87 (abbreviated as M87), a galaxy in the constellation Virgo, located 55 million light-years away from Earth. The data was collected by a large telescope virtually the size of Earth, the Event Horizon Telescope. The telescope is comprised of a network of eight radio observatories, precisely synchronized together using atomic clocks. This observation technique is called very-long-baseline interferometry (VLBI), which leverages the rotation of the Earth to form a huge telescope which observes at a wavelength of 1.3 mm, giving an angular resolution of 20 micro-arcseconds. VLBI can be used to observe other things as well; for example the Event Horizon Telescope also recorded data from our Milky Way galaxy, where there is source of radio noise called Sagittarius A*, likely the location of another black hole. However, it is much smaller than M87, making it harder to image. The Event Horizon Telescope is at work imaging the black hole in our galaxy.

Furthermore, the image of M87 also enables astronomers to measure the radius of the black hole, as it provides a clear, circular shadow. Using the radius, we can estimate the mass of the black hole. This measurement yielded an estimated mass of 6.5 billion solar masses, and can be generalized across other black hole. Since it is heavier than most other estimates, it suggests that the masses of other black holes may be underestimated.

Finally, the image confirms Einstein's theory of general relativity, which predicts that when too much matter or energy is crammed into a small enough volume, space-time would collapse, leading to an eternal gravity trap from which light cannot even escape. The black hole casts a shadow by creating a dark region when immersed in a bright region such as a disc of glowing gas. The circular shadow is caused by the gravitational bending and capture of light by the event horizon of the black hole.

General Post 2:

http://www.astronomy.com/news/2019/04/a-new-neutron-star-merger-is-caught-on-x-ray-camera

When neutron stars collide, they give off powerful signals, both on the electromagnetic spectrum and as gravitational waves. Back in October 2017, the first gravitational waves were detected from the merger of two neutron stars. This event was also the advent of multi-messenger astronomy, in which many telescopes observed a single event as different types of signals, including optical light, X-rays, and gamma rays. Now, a second neutron star merger has been detected in the form of an X-ray signal named XT2. The signal originated from a galaxy 6.6 billion light years away, in which two neutron stars combined into a single, heavier neutron star. This body is called a magnetar, because of its extremely strong magnetic field This is a completely new way to detect a neutron star merger, using X-rays instead of gravity waves. The X-ray signal matched astronomer’s predictions for this kind of event.

A neutron star is the final stage of some massive stars, formed after a supernova, when the star’s core collapses into dense, hot ball of neutrons about 20 meters in diameter. These neutron stars spin very quickly and possess very strong magnetic fields, trillions of times that of Earth. Magnetars, as their name suggests, are neutron stars with particularly strong magnetic fields, on the order of thousands of times that of regular neutron stars, or a quadrillion that of Earth. They are very rare, as only 30 of them are known. Their behavior can’t be replicated in a lab, so they need to be observed in order to learn more about them. For example, in the collision and merger of two neutron stars, a heavy neutron star emerges, which indicated that their structure is relatively resilient.

Week 9, Post 5

Worksheet 7, Problem 1

a)

The question asks us to set up the problem, which describes a small, cylindrical parcel of gas, oriented vertically in the Earth's atmosphere. The drawing is shown below.

$P_{up} = P(r)$
$P_{down} = P(r + \Delta r)$

b)

The other force acting on the parcel of gas is gravity. The Earth has a mass of $M_{\oplus}$, while the parcel has a mass of $m = \rho V = \rho(r) A \Delta r$.

$$F_G = \frac{GMm}{r^2}$$
$$F_G = \frac{G M_{\oplus} \rho(r) A \Delta r}{r^2}$$

c)

We can draw a free-body diagram to show all the forces acting on the parcel of air.

Note that pressure is force per unit area, so the forces exerted by $P_{up}$ and $P_{down}$ are therefore:

$F_{up} = A \cdot P_{up}$

$F_{down} = A \cdot P_{down}$

Since it is not moving, we know that the sum of the forces is 0. 

$$F_{net} = F_{up} - F_{down} - F_G = 0$$

$$A(P_{up}-P_{down}) = F_G$$

$$A(P(r) - P(r + \Delta r)) = \frac{G M_{\oplus} \rho(r) A \Delta r}{r^2}$$

d)

$$g = \frac{F_G}{m} = \frac{GM_{\oplus}}{r^2}$$

e)

First, we start with the result from part c)

$$A(P(r) - P(r + \Delta r)) = \frac{G M_{\oplus} \rho(r) A \Delta r}{r^2}$$

Substitute the result from d)

$$A(P(r) - P(r + \Delta r)) = \left(\frac{GM_{\oplus}}{r^2}\right) \rho(r) A \Delta r = g \rho(r) A \Delta r$$

Solve:

$$P(r + \Delta r) - P(r) = -g \rho(r) \Delta r$$

$$\lim_{x \rightarrow \infty} \frac{P(r + \Delta r) - P(r)}{\Delta r} = -g \rho(r)$$

$$\frac{dP(r)}{dr} = -g\rho(r)$$

f)

We refer to the ideal gas law given at the beginning of the worksheet: $P = nk_B T$.

We note that $n = \frac{\rho(r)}{m}$, where $n$ is the number density of particles (with units $cm^{-3}$) and $m$ is the mass of a particle of gas.

This means we can rewrite the ideal gas law as: 

$$P(r) = \frac{\rho(r) k_B T}{m}$$

$$\frac{dP(r)}{dr} = \frac{k_B T}{m} \frac{d\rho(r)}{dr} = -g\rho(r)$$

$$\frac{d\rho(r)}{dr} = \frac{-gm\rho(r)}{k_B T}$$

Integrate:

$$\int \frac{d\rho(r)}{\rho(r)} = \int\frac{-gm}{k_B T}dr$$

$$\ln(\rho(r)) = \frac{-gmr}{k_B T} + C$$

$$\rho(r) = \rho_0 e^{\frac{-gmr}{k_B T}}$$

g)

The scale height is given by: $H = \frac{kT}{\bar m g}$.

Let's check the units on this:

$$\frac{[erg][K]^{-1}[K]}{[g][cm][s]^{-2}} = \frac{[erg][s]^2}{[g][cm]}$$

$$= \frac{[g][cm]^2}{[s]^2}\frac{[s]^2}{[g][cm]} = [cm]$$

The scale height is the height $H$ added to $r$ that would make the density fall off by a factor of $1/e$. We can express this mathematically as:

$$\frac{\rho(r + H)}{\rho(r)} = \frac{1}{e}$$

$${\displaystyle \frac{\rho_0 e^{\frac{-gm(r+H)}{k_B T}}}{\rho_0 e^{\frac{-gmr}{k_B T}}} = \frac{1}{e}}$$

$${\displaystyle \exp(\frac{-gm(r+H)}{k_B T} + \frac{gmr}{k_B T}) = \exp(-1)}$$

$$1 = \frac{gmH}{k_B T}$$

$$H = \frac{k_B T}{mg}$$

h)

We are told that the Earth's atmosphere is mostly molecular nitrogen, $\text{N}_2$. Each molecule contains a total of 14 protons and 14 neutrons, for a total mass of 28 protons.

$$H_{\oplus} = \frac{kT}{\bar m g}$$

$$ = \frac{1.4 \times 10^{-16} \text{ erg K}^{-1} \cdot 300 \text{ K}}{28 \cdot 1.7 \times 10^{-24} \text{ g} \cdot 980 \text{ cm/s}^{2}} = 900360 \text{ cm} = 9 \text{ km}$$

Worksheet 8, Problem 2

a)

To solve this problem, we start with a diagram of the situation:

Depicted is a star, with a radius of $r$ and a mass shell of width $\Delta r$. The mass shell has a density of $u + \Delta u$, while the next mass shell out, at radius $r + \Delta r$, has an energy density of $u$. Notably, as radius increases, energy density decreases.

We are told that the net outwards flow of energy, $L(r)$, must equal the total excess energy in the inner shell divided by the amount of time needed to cross the shell's width $\Delta r$.

$$L(r) = \frac{\text{excess energy}}{\text{time to cross} \Delta r}$$

Since we have the relation in which energy density decreases as radius increases, we need to be careful of the sign.

$$\frac{du}{dr} = \lim_{\Delta r, \Delta u \rightarrow 0} -\frac{\Delta u}{\Delta r}$$

First, we need to solve for change in total energy. We know that change in volume is $\Delta V = \Delta r \cdot SA = 4 \pi r^2 \Delta r$. Then we can multiply this change in volume by energy density.

$$\text{energy} = u \cdot 4\pi r^2 \Delta r$$

Next, we need to solve for time to cross $\Delta r$, which is $\Delta t$.

$$\Delta t = \frac{\Delta r}{v_{diff}}$$

In question 1, we found that $v_{diff} = \frac{cD \rho \kappa}{N}$ and $N = \left(\frac{R_{\odot}}{l}\right)^2$

In this question, we know that $D = R_{\odot} = \Delta r$, so we can rewrite $v_{diff}$:

$$v_{diff} = \frac{cD\rho\kappa}{N} = \frac{c\Delta r \rho \kappa l^2}{(\Delta r)^2} = \frac{c\rho\kappa}{\Delta r \sigma^2 n_e^2}$$

Now, we can plug in and solve:

$$L(r) = \frac{energy}{\Delta t} = \frac{u \cdot 4\pi r^2 \Delta r}{\frac{\Delta r}{v_{diff}}}$$

$$= \Delta u 4 \pi r^2 v_{diff} = \frac{4\pi r^2 c \rho \kappa}{\sigma^2 n_e^2} \frac{\Delta u}{\Delta r}$$

$$L(r) = -\frac{4\pi r^2 c \rho \kappa}{\sigma^2 n_e^2}\frac{du}{dr}$$

b)

From the diffusion equation, we use the fact that the energy density of a blackbody is $u(T(r)) = aT^4$ to derive the differential equation:

$$\frac{dT(r)}{dr} \propto -\frac{L(r)\kappa \rho(r)}{\pi r^2 acT^3}$$

where $a$ is the radiation constant.

Starting with $u(T(r)) = aT^4$, we take the derivative with respect to r:

$$\frac{du(T(r))}{dT} \cdot \frac{dT(r)}{dr} = 4aT^3 \cdot \frac{dT}{dr}$$

$$\frac{du(T(r))}{dr} = 4aT^3 \frac{dT}{dr}$$

$$\frac{dT}{dr} = \frac{1}{4aT^3}\frac{du(T(r))}{dr}$$

Using $L(r)$ from part a), we can solve for $\frac{du}{dr}$:

$$\frac{du}{dr} = \frac{\sigma^2 n_e^2}{4\pi r^2 c \rho \kappa} L(r)$$

Plug this in and solve:

$$\frac{dT}{dr} = \frac{1}{4aT^3}\frac{\sigma^2 n_e^2}{4\pi r^2 c \rho \kappa} L(r)$$

$$= \frac{\sigma^2 n_e^2}{16 \pi r^2 a T^3 c \rho \kappa} L(r)$$

Then, substitute in $\sigma = \frac{\rho \kappa}{n_e}$:

$$\frac{dT}{dr} = \frac{\rho\kappa}{16\pi r^2 aT^3c} L(r)$$

This matches the form given at the beginning of the problem.

Acknowledgement: Worked with Awnit, Elton, and Simon on these problems.

General Post

https://www.nytimes.com/2019/04/10/science/black-hole-picture.html

Today’s exciting news item comes from one of our very own! A team of astronomers led by Shep Doeleman, an astronomer at the Harvard-Smithsonian Center for Astrophysics, announced that they had finally captured an image of a black hole, once thought to be unobservable because it is so dense that light cannot escape from it.

The image captured was of Messier 87 (M87), a galaxy in the constellation Virgo. In this galaxy is a black hole billions of times more massive than the sun, which unleashes a jet of energy into space about 5000 light-years away. The image of the black hole provides confirmation of a strange conclusion from Einstein’s theory of general of relativity that even he was reluctant to accept. When too much matter is forced into a small volume, the density becomes so extreme that space-time collapses. As Einstein’s theory predicts, the shadow formed by the black hole is circular in shape.

Many black holes are the remains of stars that burned through their fuel and collapsed into themselves. However, there are also supermassive black holes millions or billions of times larger than the sun that reside in the center of almost every galaxy. Scientists don’t know the origins of these black holes. They also don’t know exactly what would happen to something to an object that falls into a black hole.

The image was created from two years of computer analysis of observations from a large telescope the size of Earth, the Event Horizon Telescope. In reality, this telescope is actually a network of radio antennas comprising eight radio observatories on six mountains across four continents. These observatories recorded data on M87 for 10 days in April 2017. In order to prove that the black hole really was a black hole, scientists needed to measure the size of its shadow. This is difficult because of their distance, which makes it very hard to resolve them, even with a very large telescope. However, by combining data from multiple radio telescopes using a technique called very long baseline interferometry, the team was able to create a very large “virtual” telescope, the Event Horizon Telescope.

The accretion disk of a black hole is where matter swirls around a black hole, just at its edge. The ring of light in the image is the innermost orbit of light in the accretion disk. This provides a way for astronomers to measure the size of a black hole. This size measurement also provides an estimate for the mass of M87, which is 6.5 billion solar masses.

Day Lab Report

The Measurement of the Astronomical Unit

The purpose of this lab is to determine the length of the Astronomical Unit (AU), which is the distance between the Earth and the Sun. We will use geometry principles and an array of tools from the lab. The AU is important because it's the fundamental unit used in astronomy - if we have error, it would affect other measurements such as parallax measurements to nearby stars. This lab consists of three steps to gather the three pieces of information we need: angular size of the sun, rotational speed, and rotational period.

Step 1: Determine the Angular Size of the Sun

Procedure

The angular size of the Sun is how much of the sky the Sun covers, in degrees. From our perspective on the Earth, it looks like the sun moves around the entire Earth (360 degrees) in a day (24 hours). By measuring the amount of time it takes for the Sun to move its diameter, we can determine the angle that the Sun occupies in the sky.

We focus the Sun's light through a lens on a window onto an easel, producing a bright projected spot on the paper. First, we mark the right edge of the bright spot on the paper, and at the same time start timing on a stopwatch. Then, we wait for the Sun to move across the paper until the left edge crosses the first mark made on the paper, and stop the stopwatch. We repeat this to get 3 measurements to take an average and estimate the uncertainty.

Analysis

Because the sky was cloudy on our lab day, we used the Friday lab's measurements:

Trial 1: 02:12.65 min = 2.2108 min

Trial 2: 02:12.41 min = 2.2068 min

Trial 3: 02:10.99 min = 2.1832 min

We get a mean time of 2.200 min. A day is 24 hours (1440 minutes) and it's 360 degrees all around a circle. We can use the ratio of this time over a day to find the angular diameter.

$$\frac{\text{movement time}}{\text{length of day}} = \frac{\theta}{360^{\circ}}$$

$$\theta = \frac{360^{\circ} \times \text{movement time}}{\text{length of day}}$$

$$\theta = \frac{360^{\circ} \times 2.200 \text{ min}}{1440 \text{ min}}$$

$$\theta = 0.55^{\circ}$$

Error

The error comes from the reaction time of marking the line and pressing the stopwatch. We're using this neat reference for propagation of error by multiplication of measured quantities. According to Human Benchmark, the average reaction time is 284 ms.

$$\delta \theta = |\theta| \cdot \sqrt{\left(\frac{\delta t}{t}\right)^2}$$

$$\delta \theta = 0.55^{\circ} \cdot \sqrt{\left(\frac{0.284 \text{ s}}{2.200 \text{ s} \times 60 \text{ s/min}}\right)^2}$$

$$\delta \theta = 0.0255^{\circ}$$

Step 2: Determine the Rotational Speed of the Sun

Procedure

We can use Doppler shift to measure radial velocity. We measure NaD emission lines from the East limb and then the West limb of the Sun. Using the CCD (charge-coupled device, a digital imaging instrument) we take measurements of the Doppler shifts of these two emission lines relative to that of a Telluric (meaning terrestrial, not the element) absorption line which comes from water in the Earth's atmosphere. Since we don't know how the Sun is oriented, we take 8 measurements around the edge of the Sun in pairs: left, right; top, bottom; top right, bottom left; top left, bottom right.

The setup is tedious because we need to align the CCD with the NaD lines. We use the Sodium lamp to adjust the slit width to get the sharpest lines possible without losing out on brightness. Once everything is aligned, we remove the lamp and allow sunlight into the slit. We take images at 8 locations of the Sun's limb because we don't know the orientation of the Sun - the Sun is tilted, the Earth is tilted, and the image is reversed after being reflected off all the mirrors. The two opposite locations with the biggest shift are the ones closest to the equator of the Sun.

Analysis

We import the image files using MaximDL. Then, we draw a box from the very left edge to the very right edge, to ensure that the data will be 765 pixels wide. We have 4 pairs of datasets to compare because we want to choose the pair that will be closest-aligned to the solar equator. We can find the index of the minimum value for each dataset, and determine the pair with the biggest difference between the indices of their respective minimum values. We ended up choosing the Bottom Left/Top Right pair.

Using this pair of datasets, we want to derive the most accurate minimum positions of the NaD lines. We achieve this by fitting the curves to a quadratic $ax^2 + bx + c$ and then use the $a$ and $b$ values to calculate the minimum pixel value, using the equation $min = \frac{-b}{2a}$.

Between the left and right NaD lines, there are 343.0319 pixels separating them. We know in fact that there are 5.97 angstroms separating the two lines, which means that we have a conversion factor of $5.97 \text{ ang}/343.0319 \text{ pixel} = 0.017404 \text{ ang/pixel}$.

For the left emission line, we have a 2.478508 pixel separation, which converts to 0.043135 angstroms. For the right emission line, we have a 2.083325 pixel separation, which converts to 0.036257 angstroms.

We can take these distances and convert them into a velocity, using the equation:

$$\Delta v = \frac{\Delta \lambda \cdot c}{\lambda}$$

For the left side, we get:

$$\Delta v_{left} = \frac{\Delta \lambda \cdot c}{\lambda} = \frac{0.043135 \text{ ang} \cdot c}{5895.94 \text{ ang}} = 2.19 \text{ km/s}$$

For the right side, we get:

$$\Delta v_{right} = \frac{\Delta \lambda \cdot c}{\lambda} = \frac{0.036257 \text{ ang} \cdot c}{5889.96 \text{ ang}} = 1.85 \text{ km/s}$$

We have to be careful not to double count the rotational velocity, since one side is going away from us, while the other is going towards us. To account for this, we first take the average of the two sides, and then halve the result.

$$v_{rot} = \frac{\Delta v_{left} + \Delta v_{right} }{2} \cdot \frac{1}{2} = \frac{2.19 \text{ km/s} + 1.85 \text{ km/s}}{4} = 1.01 \text{ km/s}$$

Error

To account for error, we use the Telluric absorption line, which, as its name suggests, is purely terrestrial. It arises from water in the Earth's atmosphere. Again, we fit the curves to a quadratic to find the pixels with the minimum values.

We determine the offset to be 0.2656 pixels, so we can convert to angstroms using our conversion factor:

$$0.2656 \text{ pixels} \times 0.017404 \text{ ang/pixel} = 0.004622 \text{ angstroms}$$

Meanwhile, the average separation between left and right is as follows:

$$\Delta \lambda_{avg} = \frac{0.043135 \text{ ang} + 0.036257 \text{ ang}}{2} = 0.039696 \text{ ang}$$

Let's use the propagation of error from multiplied quantities:

$$\delta v = |v| \sqrt{\left(\frac{\delta \lambda}{\Delta \lambda_{avg}}\right)^2}$$
$$\delta v = 1.01 \text{ km/s} \cdot \sqrt{\left(\frac{0.004622 \text{ ang}}{0.039696 \text{ ang}}\right)^2}$$
$$\delta v = 0.118 \text{ km/s}$$

Step 3: Determine the Rotational Period of the Sun

Procedure

We want to determine how long it takes for the Sun to rotate about its axis. We can use sunspots as markers of the sun's surface to measure how long it takes for the sun to rotate a certain number of degrees. This relies on the assumption that the spots are carried across the Sun by its rotation.

Ideally, we would track the movement of sunspots ourselves using a heliostat. However, since the weather was cloudy, we ended up using archival data. We opened the sunspot movie on the computer and resized the window to fit in the transparency grid of the Sun. We marked the transparency as the sunspot moved across the Sun, noting the angle the sunspot traveled and the amount of time it took. We measured 3 different sunspots to take an average.

Analysis

The following are the measurements for the 3 sunspots we tracked.

Sunspot 1	Time	Degree
Start	2013-11-28_133000	30°
End	2013-12-04_210000	110°
Delta	151.5 hours (6 days, 7 hours, 30 mins)	80°
Sunspot 2	Time	Degree
Start	2013-11-27_090000	40°
End	2013-12-01_163000	90°
Delta	103.5 hours (4 days, 7 hours, 30 mins)	50°
Sunspot 3	Time	Degree
Start	2013-11-29_193000	40°
End	2013-12-03_000000	80°
Delta	76.5 hours (3 days, 4 hours, 30 mins)	40°

Assuming that time and angle scale linearly with each other, we can take the average of the times and angles to get the average rate of rotation.

$$\bar{t} = \frac{151.5 \text{ h} + 103.5 \text{ h} + 76.5 \text{ h}}{3} = 110.5 \text{ h}$$
$$\bar{\theta} = \frac{80^{\circ} + 50^{\circ} + 40^{\circ}}{3} = 56.67^{\circ}$$

We want to use the rate of rotation to determine the rotational period $P$, knowing that there are 360° in a circle.

$$\frac{P}{360^{\circ}} = \frac{t }{\theta}$$
$$P = \frac{110.5 \text{ h} \cdot 360^{\circ}}{56.67^{\circ}}$$
$$P = 702 \text{ h}$$

Therefore, it takes 702 hours for the Sun to rotate about its axis.

Error

A source of error comes from the measuring tool we used to measure angle, the transparency grid we placed on the computer monitor. The resolution we get comes from the width of the line marking 10° intervals, which are about 1° wide. Let's use the propagation of error from multiplied quantities:

$$\delta P = |P| \sqrt{\left(\frac{\delta \theta}{\theta}\right)^2}$$
$$\delta P = 702 \text{ h} \cdot \sqrt{\left(\frac{1^{\circ}}{55.67^{\circ}}\right)^2}$$
$$\delta P = 12.61 \text{ h}$$

Putting it all together

Calculation

Now we have everything we need to solve for the AU.

First, solve for the solar radius:

$$v = \frac{2 \pi R}{P}$$

$$R = \frac{vP}{2\pi} =  \frac{1.01 \text{ km/s} \times 702 \text{ h} \times 3600 \text{ s/h}}{2\pi}$$
$$R = 406239 \text{ km}$$

Then, using the skinny triangle approximation, we can determine the AU:

$$\sin(\theta/2) = \frac{R}{\sqrt{AU^2 + R^2}}$$

By approximation:

$$\theta/2 = \frac{R}{AU}$$
$$AU = \frac{R}{\theta/2} = \frac{406239 \text{ km}}{\frac{0.55^{\circ}}{2} \cdot \frac{\pi}{180^{\circ}}}$$
$$AU = 8.46392 \times 10^7 \text{ km} = 8.46392 \times 10^{12} \text{ cm}$$

Error

Let's use the propagation of error from multiplied quantities:

$$\delta AU = |AU| \cdot \sqrt{\left(\frac{\delta \theta}{\theta}\right)^2 + \left(\frac{\delta v}{v}\right)^2 + \left(\frac{\delta P}{P}\right)^@}$$
$$\delta AU = 8.46392 \times 10^{12} \text{ cm} \cdot \sqrt{\left(\frac{0.0255^{\circ}}{0.55^{\circ}}\right)^2 + \left(\frac{0.118 \text{ km/s}}{1.01 \text{ km/s}}\right)^2 + \left(\frac{12.61 \text{ h}}{702 \text{ h}}\right)^2}$$
$$\delta AU = 1.07468 \times 10^{12} \text{ cm}$$

Our final value of AU, including error, is $8.46392 \times 10^{12} \pm 1.07468 \times 10^{12} \text{ cm}$

Our calculated value of AU is within a factor of 2 of the actual measurement, $AU = 1.49597870 \times 10^{13} \text{ cm}$, which is pretty good. If we want to feel bad about ourselves, we notice that we're about 43% off.

There are many sources of error in the various procedures. In particular, the measurement of the rotational period of the Sun seemed particularly crude, because the transparency was very imprecise, and there is ambiguity as to the start or end of a sunspot. Overall, the lab managed to get a reasonably close estimate of the AU.