{ "cells": [ { "cell_type": "markdown", "metadata": {}, "source": [ "# ZPEM2509 Astrophysics Lab 101: Eclipsing Binaries \n", "\n", "Developed by Dr. Simon Murphy, UNSW Canberra ([s.murphy@adfa.edu.au](mailto:s.murphy@adfa.edu.au))\n", "\n", "Please read the background material in the lab script and complete the light curve simulations in Part II prior to starting this exercise.\n", "\n", "
\n", "\n", "# Part III: Fundamental properties of the young eclipsing binary 'THOR 42'\n", "\n", "
\n", "\n", "Now that we have an understanding of how the physical and orbital parameters of an eclipsing binary affect its light curve, we can try to solve the inverse problem and derive the fundamental properties (in particular the masses, radii and semi-major axis) of a binary using its light and velocity curves. \n", "\n", "Our target of interest in this lab is the low-mass eclipsing binary __THOR 42__ (otherwise known as CRTS J055255.7-004426), which is a member of a young group of stars called the [32 Orionis Moving Group](https://ui.adsabs.harvard.edu/abs/2017MNRAS.468.1198B/abstract). THOR 42 comprises two M dwarf stars which are smaller, less massive, fainter and cooler than the Sun. You can find more information on the system in the [Simbad database](http://simbad.u-strasbg.fr/simbad/sim-basic?Ident=THOR+42&submit=SIMBAD+search). \n", "\n", "
\n", " \n", " \n", " \n", "
\n", " Location of THOR 42 in the constellation of Orion, to the east of Orion's belt
(remember that east and west are swapped when looking out at the sky from the Earth).\n", "
\n", "\n", "THOR 42 was observed by the NASA *TESS* ([*Transiting Exoplanet Survey Satellite*](https://www.nasa.gov/tess-transiting-exoplanet-survey-satellite)) mission between 2018 December 15 and 2019 January 6 as part of its all-sky survey. We will use *TESS* data to construct a light curve for THOR 42, then analyse radial velocity measurements from the [ANU 2.3-m telescope](https://rsaa.anu.edu.au/observatories/telescopes/anu-23m-telescope) at [Siding Spring Observatory](https://goo.gl/maps/UBrdy1XS57sKuqrR9) near Coonabarabran, NSW. \n", "\n", "------\n", "\n", "First, let's load some useful Python modules to read the data and make plots. \n", "\n", "Click in the code cell so it turns green and press `SHIFT` + `RETURN` to run the code. If everything loads without error the focus will move to the next cell.\n", "\n", "
\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "scrolled": true }, "outputs": [], "source": [ "# Import various libraries\n", "%matplotlib widget\n", "from astropy.table import Table\n", "import numpy as np\n", "import matplotlib.pyplot as plt\n", "from interactive_figures import RVCurve, LightCurve, string_length" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## A. Light curve analysis\n", "\n", "
We begin by loading the *TESS* photometry from a text file. This file has 4 columns - a shortened version of the Barycentric Julian Date (BJD - 2450000), the flux (brightness) from the *TESS* pixel containing THOR 42, a flux error and a flag (0 or 1) indicating whether the measurement is good or bad. The first few rows of the file look like this:" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "!head TESS.txt # Run the Unix 'head' command to show the first 10 lines" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "\n", "Rather than Earth-centred Julian Dates (JD; the decimal number of days since noon on January 1, 4713 BC), we use *Barycentric* Julian Dates (BJD) calculated at the centre-of-mass of the solar system (approximately the centre of the Sun). This correction is necessary because otherwise the times recorded by *TESS* would be up to $\\pm$8.3 min different depending on which side of the Sun the Earth was at the time of observation.\n", "\n", "First, load the file and select the 'good' measurements:\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "scrolled": true }, "outputs": [], "source": [ "# Load the TESS photometry\n", "tess = Table.read('TESS.txt', format='ascii')\n", "# Only keep measurements with no issues (FLAG = 0)\n", "good = tess['FLAG'] == 0 # this is an array of True or False values which we can apply to the table\n", "bjd = tess['BJD'][good]\n", "flux = tess['FLUX'][good]" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "And now we plot the time-series light curve. Once the figure has been created you can interact with it using the toolbar on the left:\n", "\n", "
\n", " Zoom to rectangle \n", " Home (return plot to initial state) \n", " Pan (click and drag) \n", " Save the current view \n", "
\n", "\n", "The $x$ and $y$ coordinates of the cursor are displayed at the bottom of the plot." ] }, { "cell_type": "code", "execution_count": null, "metadata": { "scrolled": false }, "outputs": [], "source": [ "# Plot the TESS time series\n", "fig, ax = plt.subplots(constrained_layout=True, figsize=(6, 4))\n", "plt.plot(bjd, flux, 'o-', ms=2, lw=0.5, color='black', label='Good')\n", "# Also plot the rejected epochs in red. \n", "# The '~' operator does the inverse and selects the 'not good' (i.e. bad) points\n", "plt.plot(tess['BJD'][~good], tess['FLUX'][~good], 'o', ms=2, color='tab:red', label='Bad')\n", "plt.xlabel('BJD $-$ 2450000 (days)')\n", "plt.ylabel('Flux (normalised)')\n", "plt.title('$TESS$ light curve')\n", "plt.legend()\n", "plt.show()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "The alternating deep primary and shallower secondary eclipses are immediately visible. Observe that unlike the light curves in the simulator, THOR 42 does not return to a constant brightness outside of the eclipses. We will return to this point later.\n", "\n", "The red data points were rejected for failing quality checks (`FLAG = 1`). The large group of bad points between $8475 < \\mathrm{BJD} - 2450000 < 8479$ are contaminated by Earth-shine reflecting off the spacecraft's sun shade and the bad points every 3 days are due to thruster firings to dump angular momentum from gyroscopes used to keep the spacecraft pointing accurately.\n", "\n", "
\n", "\n", "__Q1. Zoom in and explore the plot. By measuring the time between successive primary or secondary eclipses, estimate an approximate period ($P$, in days). Also estimate the BJD for the middle of a single primary eclipse ($t_{0}$, BJD - 2450000). Save a copy of the plot illustrating your method.__ \n", "\n", "
\n", "\n", "We can use these values to create a *phased* (sometimes called *folded*) light curve by plotting the flux versus the orbital phase $\\phi$. This is the fraction of the orbital period that has elapsed since some reference time $t_0$ and is defined as:\n", "\n", "\\begin{align}\n", "\\phi = (\\mathrm{BJD} - t_{0})/P - \\mathrm{floor}([\\mathrm{BJD} - t_{0}]/P)\n", "\\end{align}\n", "\n", "where the $\\mathrm{floor()}$ function returns the rounded down integer value. The phase varies between 0 and 1 as the stars move in their orbits and our choice of $t_{0}$ places the primary eclipse at $\\phi = 0$ (also $\\phi = 1$). \n", "\n", "
\n", "\n", "__Q2 a. Use the figure below with your $P$ and $t_{0}$ values to create a phased light curve for THOR 42.__\n", "\n", "Remember to click inside the code cell and create the plot by pressing `SHIFT` + `RETURN`.\n", "\n", "
\n", "\n", "\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "scrolled": false }, "outputs": [], "source": [ "# Plot the interactive figure\n", "LightCurve(bjd, flux)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "\n", "Does the light curve resemble that of an eclipsing binary? (perhaps not...)\n", "\n", "Given the short (~2 hr) duration of the eclipses and the 30 min cadence of *TESS* observations, it is likely your initial choices of $P$ and $t_{0}$ are not optimal. We can assess how well the light curve is phased by computing its ['string length'](https://academic.oup.com/mnras/article/203/4/917/1029604). To do this we simply connect up all the points in order of phase and compute the total length. The best period should have the shortest string length (within some uncertainty).\n", "\n", "
\n", "\n", "__Q2 b. Turn on the string length option in the plot and explore how changing the period affects the string length. Can you find a better period with a string length less than 5? (do not spend more than a few minutes trying, we will do this more methodically below).__\n", "\n", "
\n", "\n", "To help refine the period, we can construct a crude *periodogram* by computing the string length for a range of trial periods and selecting the period with the smallest string length. \n", "\n", "Based on the unphased *TESS* light curve we can be pretty sure that the period is between 0.8 and 0.9 days. Let's compute the string length for 0.0001 day increments and plot the results versus period. The function `string_length(phase, flux)` below does the actual string length calculation:\n", "\n", "
\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "# Create a list of trial periods every 0.0001 days\n", "periods = np.arange(0.8, 0.9, 0.0001)\n", "string_lengths = [] # empty list to hold the string lengths\n", "\n", "# Loop through the trial periods, computing a string length for each one\n", "for p in periods:\n", " phase = bjd/p - np.floor(bjd/p) # don't worry about t0, it just shifts the light curve\n", " string_lengths.append(string_length(phase, flux))\n", "\n", "# Plot log(string length) versus the period to emphasize the smallest values\n", "fig, ax = plt.subplots(constrained_layout=True, figsize=(6, 4))\n", "plt.plot(periods, np.log10(string_lengths), '-', lw=1.0, color='tab:blue')\n", "plt.xlabel('Period (days)')\n", "plt.ylabel('$\\log_{10}$(string length)')\n", "plt.title('String length periodogram')\n", "plt.show()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "\n", "\n", "\n", "__Q2 c. Zoom in and record the best period, including a copy of the plot in your report. Does this new period improve your phased light curve?__\n", "\n", "
\n", "\n", "A more important problem is the natural degeneracy between $P$ and $t_0$; different values of each parameter can yield equally well-phased light curves, especially for data taken over a short time period like *TESS*. If the period is even slightly incorrect, this error is amplified when trying to phase observations taken at much earlier or later times (e.g. our radial velocity measurements).\n", "\n", "We can break the degeneracy between $P$ and $t_0$ by using a well-sampled light curve of a primary eclipse to accurately determine $t_0$ without knowing the exact period. Unfortunately our *TESS* observations are not suitable for this task. Luckily, we have 50-sec cadence observations from the ANU [SkyMapper telescope](http://skymapper.anu.edu.au/about-skymapper/) at Siding Spring Observatory, which captured a single primary eclipse on 2017 December 13:\n", "\n", "
\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "scrolled": false }, "outputs": [], "source": [ "# Load the SkyMapper photometry\n", "sm = Table.read('SkyMapper.txt', format='ascii')\n", "# Only consider the fraction of a day around the eclipse\n", "# Remember to add 8101 back to the x-axis value.\n", "sm['BJD'] = sm['BJD'] - 8101\n", "\n", "# Plot the eclipse\n", "fig, ax = plt.subplots(constrained_layout=True, figsize=(6, 4))\n", "plt.plot(sm['BJD'], sm['dmag'], 'o', ms=2, color='black')\n", "plt.xlabel('BJD $-$ 2458101.0 (days)')\n", "plt.ylabel('Differential $i$ magnitude')\n", "plt.title('SkyMapper primary eclipse (2017 Dec 13)')\n", "plt.xlim(sm['BJD'][0],sm['BJD'][-1])\n", "plt.ylim(plt.ylim()[::-1]) # Reverse the y-axis\n", "plt.show()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "This light curve has __not__ been phased, so we have no information about the period (aside from a reasonable prior estimate so we could schedule the observations!). The brightness is now also measured in *magnitudes*, not flux. Recall that a larger magnitude indicates a *fainter* measurement. Also note the different $x$-axis, which differs by 8101 days from the first *TESS* light curve.\n", "\n", "
\n", "\n", "__Q3. Use the SkyMapper observations to accurately determine the point of mid-eclipse, $t_0$. You may need to zoom in to obtain the best results. Include an annotated copy of the plot in your report. Remember to add 8101.0 days to the value you read off the plot to return it to $\\mathrm{BJD}-2450000$ like the other data sets.__\n", "\n", "\n", "
\n", "\n", "__Q4. Update the *TESS* light curve below, fixing $t_0$ to the value you just obtained and adjusting the period to align the primary eclipse with $\\phi=0$. You may need to zoom in to obtain the best results. When you are happy, record the final values of $P$ and $t_0$, and save a copy of the light curve for your report. How does the new period compare to your earlier estimates?__\n", "\n", "It should be possible to obtain a string length of <5 with the right values. Check with the demonstrator that your $P$ and $t_0$ are correct before moving on.\n", "\n", "
\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "# Plot the interactive figure\n", "LightCurve(bjd, flux, P=0.8589, t0=8101.0, dP=1e-6)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "\n", "\n", "__Q5. Use your phased light curve to answer the following questions. You may wish to refer back to the light curve simulator for guidance.__ \n", "\n", "- __Assuming the stars have different radii, what does the shape of the primary eclipse tell us about the inclination of the system?__\n", "- __What do the different eclipse depths tell us about the temperatures of the stars?__\n", "- __What does the phase of the secondary eclipse tell us about the eccentricity of the orbits?__ \n", "\n", "
\n", "\n", "As we mentioned earlier, comparing the light curve to those from the simulator you will notice THOR 42 has a broad dip in brightness around primary eclipse. This is probably due to a large, cool (therefore darker) group of star spots on the surface of the primary star which rotates into view just before the eclipse. Since THOR 42 is such a compact binary we expect the rotation rate of each star to be the same as the orbital period (this is called *tidal synchronisation*) so the spots will never fall out of phase with the eclipse. The same effect is seen with the Moon keeping the same face towards the Earth during its orbit.\n", "\n", "\n", "------\n" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## B. Radial velocity curves\n", "\n", "
Now that we have a good idea of the period and eclipse timing we can analyse the radial velocity data. We again read it in from a text file containing the time of observation (BJD - 2450000), the radial velocities for each star (in km s$^{-1}$), their uncertainties and a quality flag:\n", "\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "!head RV.txt # Display the first 10 rows of the file" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "First, load the data and select the 'good' measurements:" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "scrolled": true }, "outputs": [], "source": [ "# Load the RV data \n", "rv = Table.read('RV.txt', format='ascii')\n", "t1 = rv['BJD']\n", "rv1 = rv['RV1']\n", "\n", "# Inflate the velocity errors to force the reduced chi-squared\n", "# to be approximately 1 for the best fitting model.\n", "# This is not best practice, see https://arxiv.org/pdf/1012.3754.pdf\n", "e_rv1 = np.sqrt(rv['E_RV1']**2 + 3.5**2)\n", "\n", "# Use only those secondary epochs with good velocities (FLAG = 0)\n", "good = rv['FLAG'] == 0 \n", "t2 = rv['BJD'][good]\n", "rv2 = rv['RV2'][good]\n", "# Again, inflate the errors\n", "e_rv2 = np.sqrt(rv['E_RV2'][good]**2 + 5.0**2)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "\n", "We reject those measurements that were taken too close to the eclipses and have erroneous secondary velocities (`FLAG > 0`). This leaves 65 primary and 52 secondary velocities.\n", "\n", "\n", "Next, plot the radial velocities of both stars versus time:\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "scrolled": false }, "outputs": [], "source": [ "fig, ax = plt.subplots(constrained_layout=True, figsize=(6, 4))\n", "plt.axhline(0, color='tab:grey', lw=1)\n", "plt.plot(t1, rv1, 'o', ms=4, color='tab:red', label='Primary')\n", "plt.plot(t2, rv2, 'o', ms=4, color='tab:blue', label='Secondary')\n", "plt.xlabel('BJD $-$ 2450000 (days)')\n", "plt.ylabel('Radial velocity (km s$^{-1}$)')\n", "plt.title('WiFeS radial velocities')\n", "plt.legend()\n", "fig.show()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "As you can see, the first measurement was made at $\\mathrm{BJD}\\approx2457319$ (2015 October 23) and the velocities had swapped signs at the time the second measurement was made almost a year later. This immediately identified THOR 42 as a spectroscopic binary and it was observed frequently over the following 18 months, often multiple times per night. \n", "\n", "
\n", "\n", "Both stars will obey Kepler's Laws as they orbit each other. For stars moving in elliptical orbits with period $P$, eccentricity $e$, inclination to the line of sight $i$ and semi-major axes $a_1$ and $a_2$, it can be shown that the *observed* radial velocities $v_{r}(t)$ for each star should obey the following relations:\n", "\n", "\\begin{align}\n", "v_{r}(t) &= v_{\\mathrm{sys}} + K[e\\cos\\omega + \\cos(\\omega + \\phi)] \\\\\\\\\n", "\\mathrm{with}\\ K_1 &= \\frac{2\\pi a_1 \\sin i}{P \\sqrt{1-e^{2}}} \\ \\ \\mathrm{and}\\ \\ K_2 = \\frac{2\\pi a_2 \\sin i}{P \\sqrt{1-e^{2}}}\\\\\n", "\\end{align}\n", "\n", "$v_{\\mathrm{sys}}$ is the constant *systemic* or centre-of-mass velocity and $\\omega$ relates to the angle between periastron (closest approach between the stars) and the line of sight. The angles $\\omega$ and $i$ orient the orbits relative to the plane of the sky. For circular orbits ($e=0$) the expression simplifies to a simple cosine and $\\omega$ is undefined (as the stellar separation is constant). $K_1$ and $K_2$ are called the radial velocity *semi-amplitudes*. We cannot be sure of the *true* semi-amplitudes because the orbital inclination $\\sin i$ is still unknown, but for an eclipsing system we know it must be close to $i=90^{\\circ}$.\n", "\n", "\n", "
\n", "\n", "\n", "__Q6 a. Use the figure below to create phased velocity curves for THOR 42 with your improved $P$ and $t_0$ values. Then adjust the sliders to fit radial velocity models $v_{r}(t)$ to the measurements (you can also click on the model values and adjust them manually). The values of $K_1$, $K_2$ and $v_{\\mathrm{sys}}$ are in km s$^{-1}$ and $\\omega$ is measured in degrees.__\n", "\n", "
\n", "\n", "The numbers overlaid on the plot are the [reduced $\\chi^2$ (chi-squared) statistic](https://en.wikipedia.org/wiki/Reduced_chi-squared_statistic), which is a measure of how well a particular model (e.g. the radial velocity equation) fits the observations and their uncertainties. A $\\chi_{r}^2$ value close to 1.0 generally indicates a good fit to the observations. It should be possible to obtain $\\chi_{r}^2$ values of <1.5 for both sets of velocities.\n", "\n", "
\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "scrolled": false }, "outputs": [], "source": [ "# Plot the interactive figure\n", "RVCurve(t1, rv1, e_rv1, t2, rv2, e_rv2)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "\n", "__Q6 b. When you are happy with the results, record the values of $K_1$, $K_2$, $v_{\\mathrm{sys}}$ (all in km s$^{-1})$, $e$ (dimensionless) and $\\omega$ (in degrees), and save a copy of the velocity curves for your report. Is the eccentricity consistent with the estimate you made using the light curve in Q5?__\n", " \n", "
\n", "\n", "We now (finally!) have everything we need to determine the masses and radii of both stars.\n", "\n", "------\n", "\n", "
" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## C. Masses and orbital semi-major axis of THOR 42\n", "\n", "
From centre-of-mass arguments, it can be shown that for two stars of masses $M_1$ and $M_2$, moving in elliptical orbits with semi-major axes $a_1$ and $a_2$;\n", "\n", "\\begin{align}\n", "a_1 M_1 = a_2 M_2\\quad\\mathrm{and}\\quad M_1 K_1 = M_2 K_2\n", "\\end{align}\n", "\n", "for radial velocity semi-amplitudes $K_1$ and $K_2$. We can therefore immediately calculate the *mass ratio* of the two stars:\n", "\n", "\\begin{align}\n", "\\frac{M_2}{M_1} = \\frac{K_1}{K_2}\n", "\\end{align}\n", "\n", "
\n", "\n", "__Q7. Using your fitted values for $K_1$ and $K_2$, calculate the mass ratio $M_2/M_1$ of THOR 42.__\n", "\n", "
\n", "\n", "We can also calculate the sum of the masses by considering Kepler's third law, namely:\n", "\n", "\\begin{align}\n", "\\frac{a^3}{P^2} &= \\frac{G(M_1 + M_2)}{4\\pi^2}\n", "\\end{align}\n", "\n", "where the secondary star appears to orbit around the primary star with semi-major axis $a=a_1 + a_2$ (or vice versa).\n", "\n", "By rearranging the expressions for $K_1$ and $K_2$ from Part B, we see that:\n", "\n", "\\begin{align}\n", "a_1 = \\frac{K_1 P \\sqrt{1-e^2}}{2\\pi \\sin i} \\quad\\mathrm{and}\\quad a_2 = \\frac{K_2 P \\sqrt{1-e^2}}{2\\pi \\sin i} \n", "\\end{align}\n", "\n", "
\n", "\n", "__Q8. Derive an expression for $a=a_1+a_2$ in terms of $(K_1,K_2,e,i,P)$ and use this to calculate the separation of the stars in solar radii $(1 R_{\\odot} = 6.957\\times10^{8}\\ \\mathrm{m})$ and astronomical units $(1\\ \\mathrm{au} = 1.4959\\times10^{11}\\ \\mathrm{m})$, assuming circular orbits and an inclination of $i=85^{\\circ}$. Be aware of units in your calculation, noting that $K_1$ and $K_2$ are in km s$^{-1}$ and $P$ is in days.__\n", "\n", "
\n", "\n", "__Q9 a. Using Kepler's third law and your expression for $a$, show that the sum of the masses depends on the inclination and is given by the following:__\n", "\n", "\\begin{align}\n", "(M_1 + M_2)\\sin^3 i = \\frac{P}{2\\pi G}(1-e^2)^{3/2}(K_1 + K_2)^3\n", "\\end{align}\n", "\n", "
\n", "\n", "__Q9 b. Calculate the sum of the masses in solar units $(1 M_{\\odot} = 1.989\\times10^{30}\\ \\mathrm{kg})$, assuming an inclination of $i=85^{\\circ}$ and circular orbits.__\n", "\n", "
\n", "\n", "__Q10. Use the sum and ratio of the masses to calculate the individual masses of THOR 42 in solar units.__\n", "\n", "------\n", "\n", "
" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## D. Radii of THOR 42\n", "\n", "
The *relative* radii of both stars (i.e. $R/a$, for semi-major axis $a$) can be extracted from the light curve and then multiplied by $a$ to calculate the *absolute* radii. The most accurate values come from a full physical model of the binary, taking into account things like limb darkening, non-spherical stars, star spots etc. Nevertheless, we can still extract reasonable values using nothing more than geometry and some basic assumptions.\n", "\n", "The duration of the eclipses tells us about the *sum* of the radii, as the smaller, fainter star must transit across the disk of the brighter star during primary eclipse, and is occulted by the primary star during secondary eclipse. The following discussion is taken from [work on transiting exoplanets](https://www.paulanthonywilson.com/exoplanets/exoplanet-detection-techniques/the-exoplanet-transit-method/) but is also applicable to eclipsing binaries. Consider an eclipsing binary as seen from Earth the moment the primary eclipse starts:\n", "\n", "
\n", "\n", "
\n", "\n", "In our case $R_{\\star} = R_1$ and $R_{p}= R_2$. The secondary star is moving right to left. Using Pythagoras, the straight-line distance the secondary star has to travel across the disk of the primary, as projected onto the plane of sky, can be expressed as\n", "\n", "\\begin{align}\n", "2\\ell = 2\\sqrt{(R_1 + R_2)^2 - (b R_1)^2}\n", "\\end{align}\n", "\n", "where the distance between stellar centres, $b=a\\cos (i)/R_1$, is called the *impact parameter* of the eclipse and is 0 at the centre of the disk of the primary and 1 at the limb. For a completely edge-on orbit ($i=90^{\\circ}$) the impact parameter is zero. \n", "\n", "Now consider how the situation looks in three dimensions:\n", "\n", "
\n", "\n", "
\n", "\n", "During the primary eclipse the secondary moves from $A$ to $B$ around its orbit, creating an angle $\\alpha$ (measured in radians) with respect to the centre of the primary star. If we assume the orbit is circular (a reasonable assumption for THOR 42), the distance around an entire orbit is $2\\pi a$, where $a$ is now the radius of the orbit. The arc length between $A$ and $B$ is $\\alpha a$ and the distance along a straight line between $A$ and $B$ is $2\\ell$, as calcuated above. \n", "\n", "From the triangle formed by $A$, $B$ and the centre of the primary star: \n", "\n", "\\begin{align}\n", "\\sin(\\frac{\\alpha}{2}) = \\frac{\\ell}{a}\n", "\\end{align}\n", "\n", "and therefore the transit duration $\\Delta t$ (in time units) is given by:\n", "\n", "\\begin{align}\n", "\\Delta t & = P \\frac{\\alpha}{2\\pi} \\\\\n", " & = \\frac{P}{\\pi}\\sin^{-1}(\\frac{\\ell}{a}) \\\\\n", " & = \\frac{P}{\\pi}\\sin^{-1}(\\frac{1}{a}\\sqrt{(R_1 + R_2)^2 - (b R_1)^2}) \\\\\n", " & = \\frac{P}{\\pi}\\sin^{-1}(\\frac{1}{a}\\sqrt{(R_1 + R_2)^2 - (a \\cos i)^2}) \n", "\\end{align}\n", "\n", "where $\\sin^{-1}$ is the inverse sine function (sometimes called asin or arcsin) and $\\sin(\\sin^{-1}(x)) = x$.\n", "\n", "
\n", "\n", "__Q11. Rearrange the above equation for $\\Delta t$ to show that the sum of the relative radii is given by__\n", "\n", "\\begin{align}\n", "\\frac{R_1}{a} + \\frac{R_2}{a} & = \\sqrt{\\cos^2 i + \\sin^2(\\pi \\Delta\\phi)}\n", "\\end{align}\n", "\n", "__for eclipse duration $\\Delta\\phi = \\Delta t\\ /\\ P$ ($\\Delta \\phi$ now in dimensionless phase units).__ \n", "\n", "
\n", "\n", "__Q12. From your phased light curve in Q4, measure the primary and secondary eclipse durations $\\Delta \\phi$ and calculate an average value. Remember to be consistent in how you define the start and end of each eclipse. You may wish to include in your report annotated copies of the eclipse light curves to illustrate your method.__\n", "\n", "
\n", "\n", "__Q13 a. Use this average $\\Delta \\phi$ and the equation from Q11 to calculate the sum of the *relative* radii, $\\frac{R_1}{a} + \\frac{R_2}{a}$, assuming an inclination of $i=85^{\\circ}$. Remember to convert the inclination to radians (or convert $\\pi \\Delta\\phi$ to degrees) in order to correctly apply the equation.__\n", "\n", "
\n", "\n", "__Q13 b. Then, using the value for the semi-major axis $a$ you found in Q8, convert this to the sum of the *absolute* radii, $R_1+R_2$, in solar units.__\n", "\n", "
\n", "\n", "The last remaining parameter to determine is the ratio of the radii, $R_2/R_1$. We can obtain this from the eclipse *depths*, assuming that the eclipses are total. Since our eclipses are not flat-bottomed we know THOR 42 cannot be exactly at $i=90^{\\circ}$. However, if $R_2$ and the impact parameter $b$ are small enough they can still cause a near-total eclipse at $i\\lesssim90^{\\circ}$. \n", "\n", "Consider the following toy light curve:\n", "\n", "
\n", "\n", "
\n", "\n", "Assuming the (circular) stellar disks of area $\\pi R^2$ emit uniformly (i.e. no limb darkening), then the luminosity $L$ (in W or erg s$^{-1}$) from each star directed towards the Earth is $L= \\pi R^2 F$ for the surface flux $F$ (in W m$^{-2}$ or erg s$^{-1}$cm$^{-2}$). Outside of the eclipses we see the full luminosity from both stars, $L_1 + L_2$, and the observed brightness $B_0$ recorded at the telescope is:\n", "\n", "\\begin{align}\n", "B_0 & = c\\ (L_1 + L_2) \\\\\n", " & = c\\ \\pi(R_1^2 F_1 + R_2^2 F_2)\n", "\\end{align}\n", "\n", "where the constant $c$ relates the flux leaving the star to the signal measured by the detector and depends on things like the inverse square of the distance, the amount of interstellar absorption, atmospheric extinction and the workings of the detector. During secondary eclipse ($\\phi=0.5$) when the secondary star is completely blocked by the primary star, the brightness recorded is just that of the primary star:\n", "\n", "\\begin{align}\n", "B_2 = c\\ \\pi R_1^2 F_1\n", "\\end{align}\n", "\n", "and during primary eclipse ($\\phi=0$) we add the stars together but subtract the light lost from the primary star as the secondary passes across it:\n", "\n", "\\begin{align}\n", "B_1 = c\\ \\pi (R_1^2 F_1 + R_2^2 F_2 - R_2^2 F_1)\n", "\\end{align}\n", "\n", "\n", "
\n", "\n", "__Q14 a. Use the equations for $B_0$, $B_1$ and $B_2$ to show that the ratio of eclipse depths does not depend on $c$ and is equal to the surface flux ratio:__\n", "\n", "\\begin{align}\n", "\\frac{B_0-B_1}{B_0-B_2} & = \\frac{F_1}{F_2} \\\\\n", "\\end{align}\n", "\n", "__Q14 b. Then use the definition of luminosity ($L\\propto R^2 F$) to show that__\n", "\n", "\\begin{align}\n", "\\frac{R_2}{R_1} = \\sqrt{\\frac{L_2}{L_1}\\frac{F_1}{F_2}}\n", "\\end{align}\n", "\n", "
\n", "\n", "__Q15. Carefully measure both eclipse depths from the phased light curve and calculate the ratio of the radii $R_2/R_1$. As the *TESS* fluxes have already been normalised, you can assume $B_0 = 1.0$ for both eclipses (i.e. ignore the dip due to the spots).__ This calculation also requires knowing the *luminosity* ratio of the stars, $L_2/L_1$, of which we have some prior knowledge from spectroscopy. As seen from Earth, the secondary appears roughly 1/5 as bright as the primary at visible wavelengths, or more accurately, $L_2/L_1 = 0.22$.\n", "\n", "
\n", "\n", "__Q16. Combine the sum and ratio of the radii to derive the individual radii of THOR 42 in solar units.__\n", "\n", "------\n" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## E. Conclusions\n", "\n", "
\n", "\n", "__Q17. Construct a summary table of results for your report. Include the period, eccentricity and semi-major axis of the system, as well as the individual masses and radii, with appropriate units and significant figures.__\n", "\n", "
\n", "\n", "Input your mass, radius and separation measurements into the light curve simulator. You may assume circular orbits, $i=85^{\\circ}$ and temperatures of $T_1 = 3100\\ \\mathrm{K}$ and $T_2 = 3000\\ \\mathrm{K}$. The light curve should look broadly similar to the *TESS* data (with the exception of the dip due to the spots). \n", "\n", "
\n", "\n", "__Q18. Take screenshots to include in your report showing the configuration of THOR 42 as viewed from Earth at phases $\\phi=0$ (primary eclipse) and $\\phi=0.5$ (secondary eclipse). Also include a plot showing the plane of the orbits as viewed from above.__\n", "\n", "
\n", "\n", "__Q19. Was our assertion of total eclipses valid at the assumed inclination of $i=85^{\\circ}$? Why or why not?__\n", "\n", "
\n", "\n", "__Q20. In a short paragraph, briefly summarise your analysis and findings for THOR 42, stating your main results and including answers to the following questions:__\n", "\n", "- Which physical quantities can be derived from the *light curve* of an eclipsing binary? What about the *radial velocity curves*? \n", "- What were the major sources of uncertainty in your analysis? This can include limitations of the data or techniques, as well as any assumptions made.\n", "- How might you improve your parameter estimates using different techniques and/or observations?\n", "\n", "
\n", "\n", "__Submit a PDF of your report, including the answers to Part II, via Moodle before the agreed deadline.__\n", "\n", "*Extra for those who are interested:* You may wish to compare your findings to recent published values for THOR 42 made using much the same observational data. Table 6 in [this paper in the Monthly Notices of the Royal Astronomical Society](https://arxiv.org/pdf/1911.05925.pdf) contains estimates of the system parameters and their uncertainties. Some of the authors may be familiar :) How do your values and analysis methods compare?\n", "\n", "


\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "scrolled": true }, "outputs": [], "source": [] } ], "metadata": { "kernelspec": { "display_name": "Python 3", "language": "python", "name": "python3" }, "language_info": { "codemirror_mode": { "name": "ipython", "version": 3 }, "file_extension": ".py", "mimetype": "text/x-python", "name": "python", "nbconvert_exporter": "python", "pygments_lexer": "ipython3", "version": "3.7.2" } }, "nbformat": 4, "nbformat_minor": 4 }