Written by Elvi Haaramaki and Sof Hachem

Star clusters

In this project, we are looking at objects called star clusters, which are groups of stars that are gravitationally bound together. There are two types – globular clusters and open clusters.

Globular clusters are large collections of thousands of millions of red stars, most of which are population II stars formed a few million years after the Big Bang. The stars are densely packed together, and the shape of the cluster is roughly spherical¹. 

Open clusters are made up of a few hundred stars which are much less densely packed together than globular clusters, and are without a distinct shape. These stars are younger blue stars, tens of millions of years old, and are weakly gravitationally bound². For this project we took pictures of M35, M36, and M37, which are all open clusters. 

Open clusters are ideal for the purposes of the project, since the stars are further apart from each other, making them easier to resolve individually. This is necessary to be able to find the apparent magnitudes of the stars. Additionally, it is assumed that all the stars in the cluster are the same distance from us, allowing us to calculate the absolute magnitudes.

HR diagrams

HR diagrams are a useful tool in astronomy, as plotting stars within a cluster based on their magnitude and colour-temperature allows the study of the stars’ relative stellar evolutionary stages and provides an idea about the age and the relative distances of given clusters³, which was also the objective during the project. 

For the data here, green and blue filters⁴ were used to attain the difference in the magnitudes corresponding to a colour index in the observational HR diagram, being plotted against the apparent (green) magnitude.  

The upper left region in the diagrams corresponds to the bluer and brighter stars, while the redder and dimmer stars are located at the bottom right. The stars in a HR diagram for a relatively young open cluster, which mostly contains stars burning hydrogen into helium, are expected to lie along the main sequence branch – this is seen in the HR diagram for M36.  In older clusters, the main sequence cut-off becomes more visible as the more massive stars burn hydrogen at faster rates and are consequently spending less time in the main sequence, leading older clusters to have more prominent red-giant branches. In globular clusters, more stellar evolutionary stages are present due to the clusters age, which can be seen in the HR diagram for M3.   

Identifying the turnoff point from the diagram and comparing the mass and brightness of the most massive remaining star in the main sequence, allows us to further estimate the age of the star and – as stars within the cluster are estimated to have formed around the same time – the age of the cluster⁵. 

References: 

1: https://www.britannica.com/science/globular-cluster 

2: https://www.schoolsobservatory.org/learn/astro/stars/clusters/openclus 

4: https://www.astronomik.com/en/photographic-filters/deep-sky-rgb-colour-filters.html

3: https://en.wikipedia.org/wiki/Hertzsprung%E2%80%93Russell_diagram 

5: https://www.scientificamerican.com/article/how-do-scientists-determi/ 

6: https://jila.colorado.edu/~ajsh/courses/astr1120_03/text/chapter4/l4S8.htm (Figure 6)

By Ceinwen Cheng

The most typical way to classify stars is based on their spectral profiles. Electromagnetic radiation is passed through filters displaying lines on a spectrum. The intensity of each spectral line provides information on the abundance of the element and the temperature of the star’s photosphere.

The system for classification is named the Morgan-Keenan (MK) system, using letters O, B, A, F, G, K and M, where O-type is the hottest and M-type is the coolest. Each letter class has a further numerical subclass from 0 to 9 where 0 is the hottest and 9 is the coolest. In addition to this, a luminosity class (Yerkes Spectral Classification) can be added using roman numerals, based on the width of the absorption lines in the spectrum which is affected but the density of the star’s surface.

The aim of  our very last imaging session was to image spectra, including stars from all the classes in the MK-system:

• 0-Class: Alnitak
• B-Class: Regulus and Alnilam
• A-Class: Alhena and Sirius
• F-Class: Procyon
• G-Class: Capella
• K-Class: Alphard
• M-Class: Betelgeuse

For each star, two sets of spectra were recorded at a 60 s exposure set on 2×2 binning. As the exposure time is relatively long, it was important to manually keep the star between the cross-hairs of the eyepiece while imaging, as there is a tendency for the star to stray from the center of the telescope as the earth rotates.

There is a trend of worsening spectra as you move down the spectral types. As luminosity is proportional to the fourth power of temperature, cooler stars in the lower end of the spectral classes such as K and M-classes are dimmer than an O-class. Displayed below are the spectra we collected for each star, excluding Alphard and Betelgeuse as no identifiable spectral lines were seen for both stars.

Alnilam B-class: For this class we expect medium amounts of Hydrogen Balmer lines, and some neutral Helium lines.  We observe spectral lines of    H-γ, H-δ,  He I, and C III in our experimental spectra analysis.

Alnitak O-class:  In the hottest class of stars, we expect to see ionized helium features. In this graph, we see prominent absorption lines of H-β, H–γ, H–δ,  as well as He I, He II and He III, which are in correlation with our expectations.

Capella G-class: We expect heavier elements such as calcium and for the Hydrogen Balmer lines to be less prominent. Compared to stars in A or F-classes, our experimental Capella spectra have less defined H-ε, H–γ, H–δ absorption lines, but a more obvious Ca II absorption line.

Sirius A-class: This class has the strongest features of the Hydrogen Balmer series. We observe H-β, H–γ, H–δ, H–ε strong spectral lines, and slight H-ζ and H-η absorption lines.  Additionally, characteristically of A-class stars, Sirius displays spectral lines of heavier elements such as Mg II and Ca II.

Alhena A-class: Alhena’s spectra are very similar to Sirius’s, being in the same class.

Regulus B-class: Displaying prominent H-β, H–γ, and H–δ Balmer lines, it is in the lower temperature end of B-class, almost an A-class. It has strong Balmer features but still a characteristic He I absorption line that places it in the B-class.

Procyon F-class: We expect to see ionized metals and weaker hydrogen lines than A-classes. Our data shows weaker H-β, H–γ, and H–δ Balmer lines than Alhena or Sirius, but more prominent absorption lines for ionized elements such as Ca II .

Stars can be assumed to be black bodies, and its absorption spectra will overall take the shape of its black body radiation. After extracting the absorption spectra, we attempted to normalize it by dividing the spectra data by its polynomial interpolation to the fourth degree. The polynomial interpolation is our estimate of the black body curve that would be displayed by the star.

Characteristically, an ideal normalize spectrum will not deviate from 1.0 on the y-axis, other than peaks where the absorption lines are. We can directly compare that the peaks in the normalized graphs are where the troughs in the initial spectra data are.

Normalized spectra, Alhena A-class

Normalized spectra, Alnilam B-class

Normalized spectra, Alnitak O-class

Normalized spectra, Capella G-class

Normalized spectra, Procyon F-class

Normalized spectra, Regulus B-class

Normalized spectra, Sirius A-class

Code: Extracting spectra and normalizing data

From here, we will go through how we extracted the spectra from the .fit files using Python, as well as the normalization of the spectra through a polynomial interpolation. This is for 1 star, specifically Alhena, where the code is repeated for the rest of the stars as well, where the only deviation is the addition of absorption lines manually.

Step 1:

#importing relevant modules
import numpy as np
import os
import math
from PIL import Image
from PIL import ImageOps
import astropy
%config InlineBackend.rc = {}
import matplotlib
import matplotlib.pyplot as plt
%matplotlib inline
from astropy.io import fits
from astropy.nddata import Cutout2D
import os
from numpy import asarray
import matplotlib.image as mpimg
from specutils.spectra import Spectrum1D, SpectralRegion
from specutils.fitting import fit_generic_continuum
from numpy import *
import numpy.polynomial.polynomial as poly
from scipy import interpolate

Step 2:

#opening both files to check stats, and positioning a cut out around both spectra
hdulista1 = fits.open(r"C:\Users\Ceinwen\Desktop\ProjectPics\Spectra\alhena_953.fit")
hdulista1.info()
print(repr(hdulista1[0].header))
dataa1 = ((hdulista1[0].data)/256)

positiona1 = (838.5, 660)
sizea1 = (150,1677)
cutouta1 = Cutout2D(np.flipud(dataa1), positiona1, sizea1)
plt.imshow(np.flipud(dataa1), origin='lower', cmap='gray')
cutouta1.plot_on_original(color='white')
plt.colorbar()



hdulista2 = fits.open(r"C:\Users\dhruv\Desktop\Project Pics\Spectra\alhena_954.fit")
dataa2 = ((hdulista2[0].data)/256)

positiona2 = (838.5, 690)
sizea2 = (150,1677)
cutouta2 = Cutout2D(np.flipud(dataa2), positiona2, sizea2)
plt.imshow(np.flipud(dataa2), origin='lower', cmap='gray')
cutouta2.plot_on_original(color='white')
plt.colorbar()

We open both .fit files to check the stats of the data, and cut out the part of the image with the spectra for both data sets, as you can see, the spectra are very faint.

Step 3:

#plotting the cut out

plt.imshow(cutouta1.data, origin='lower', cmap='gray')
plt.colorbar()

plt.imshow(cutouta2.data, origin='lower', cmap='gray')
plt.colorbar()

Plotting the cut-out, taking a closer look, the spectra lines are more visible now.

Step 4:

#converting the file's data into an array and plotting it

xa1 = cutouta1.data
a1 = np.trapz(xa1,axis=0)
plt.plot(a1)

xa2 = cutouta2.data
a2 = np.trapz(xa2,axis=0)
plt.plot(a2)

Step 5:

#finding the mean of both data and plotting it, alongside spectral lines.

a = np.mean([a1, a2], axis=0)
plt.plot(a,"r")
plt.axvline(x=206,color='b',ls=":")
plt.axvline(x=271,color='b',ls=":")
plt.axvline(x=325,color='b',ls=":")
plt.axvline(x=368,color='b',ls=":")
plt.axvline(x=525,color='b',ls=":")
plt.axvline(x=805,color='b',ls=":")
plt.axvline(x=967,color='b',ls=":")
plt.axvline(x=1401,color='b',ls=":")

Step 6:

#doing a line of best fit through the spectral lines's positions.

x1=np.array([206,271,325,368,525,805,967,1401])
y1=np.array([3835.384,3889.049,3933.66,3970.072,4101.74,4340.462,4481.325,4861.3615])
m1, b1 = np.polyfit(x1, y1, 1)
plt.scatter(x1,y1,color="red")
plt.plot(x1,m1*x1+b1)
print(m1)
print(b1)

Step 7:

#mapping the pixel number to wavelength
xo1 = np.arange(1,1678)
func1 = lambda t: (m1*t)+b1
xn1 = np.array([func1(xi) for xi in xo1])
plt.plot(xn1,a,"r")
plt.xlabel ('Wavelength ($\AA$)')
min(xn1),max(xn1)

Step 8:

#plotting the spectra in term of wavelength and including named spectral lines.
plt.plot(xn1,a,"r")
plt.xlabel ('Wavelength ($\AA$)')
plt.axvline(x=4861.3615,color='orange',label="H-\u03B2",ls=":")
plt.axvline(x=4481.325,color='magenta',label="Mg II",ls=":")
plt.axvline(x=4340.462,color='g',label="H-\u03B3",ls=":")
plt.axvline(x=4101.74,color='b',label="H-\u03B4",ls=":")
plt.axvline(x=3970.072,color='c',label="H-\u03B5",ls=":")
plt.axvline(x=3933.66,color='lime',label="Ca II",ls=":")
plt.axvline(x=3889.049,color='m',label="H-\u03B6",ls=":")
plt.axvline(x=3835.384,color='y',label="H-\u03B7",ls=":")
plt.legend()

Step 9:

#doing a polynomial fit to the 4th degree, plotting it on the same graph as the spectra
coefs1 = poly.polyfit(xn1, a, 4)
ffit1 = poly.Polynomial(coefs1)
plt.plot(xn1, ffit1(xn1))
plt.plot(xn1,a,"r")
plt.xlabel ('Wavelength ($\AA$)')

Step 10:

#diving the spectra by the polynomial fit to normalise, and plotting normalised graph

plt.plot(xn1,ffit1(xn1)/a)
plt.xlabel ('Wavelength ($\AA$)')

Finally, dividing the spectra by the polynomial interpolation to normalize the spectra.

Brief Note:

The end of the 2020 third year project came quicker than most years before us, unexpected situations such as strikes and social distancing due to coronavirus have unfortunately prevented us from collecting any more data on the variations in the apparent magnitude of Betelgeuse, detailed in the previous blog post. We have had to cancel our poster and presentation as well.

We are all saddened by the fact we can’t clamber about the dark rooftop that overlooks the Thames and half of London anymore, but I believe this project has inspired many in the group to go further in the field of Astronomy. We all want to thank our absolutely brilliant supervisor, Prof. Malcolm Fairbairn, for his guidance and leadership. He has been looking out for us in more than this project, being an amazing mentor.

Written by Ceinwen Cheng

Betelgeuse, what used to be the 10^th brightest star in the night sky, lies on the left shoulder of Orion’s constellation. This reddish star attribute’s its color to its classification of being a red supergiant around 15 times the mass of the Sun.

However, ever since October 2019, Betelgeuse’s apparent magnitude has had a drastic decrease that can be observed with the naked eye, dimming from an apparent magnitude of 0.6 in October 2019 to 1.8 in early February. The apparent magnitude scale is reverse logarithmic, i.e the lower the magnitude, the brighter the star. A star of apparent magnitude 1.0 is 2.512 times brighter than a star of apparent magnitude 2.0.

A graph published by Nature news article in February 26th c 2020 displaying a graphical representation of Betelgeuse’s apparent magnitude

Additionally, as recently as February 22^nd, Betelgeuse has started brightening again, an increase of almost 10% from the dimmest point. There are many proposed explanations, one being a change in extinction and therefore not a change in the actual luminosity of the star, which is supported by the fact that the infrared observations have of Betelgeuse has barely changed since the dimming phenomenon. In astronomy, extinction is the absorption and scattering of light particles by cosmic dust between the observer and the light emitter.

We acknowledge several problems measuring the apparent magnitude of Betelgeuse in Central London, one of the most light polluted cities in Europe. If we measure Betelgeuse directly, our data will be skewed by a range of obscure interference such as clouds, atmospheric seeing, thermal turbulence, air and light pollution. These factors are hard to mitigate and could alter or change the trend of our data. Therefore, we approach this project by also tracking a star close to Betelgeuse that apparent magnitude isn’t changing much, Betelgeuse’s sister star Bellatrix. We are interested in the ratio of luminosity between Bellatrix and Betelgeuse.

Imaging both Bellatrix and Betelgeuse in 5 different sessions in 2020, on the 28^th Jan, 3^rd Feb, 11^th Feb, 17^th Feb, 2^nd  March 2020,  we collected the images to be processed into a trend. Using Python 3, Dhruv Gandotra, our python expert, extracted the data from our images to produce the following trend in the graph below, showing the trend of the ratio of magnitudes between Betelgeuse and Bellatrix. Additionally, he also found a source online [twitter handle @Betelbot] that tracks the apparent magnitude of Betelgeuse over the same time frame as ours, allowing us to peer review our data.

@Betelbot ‘s measure of Betelgeuse apparent magnitude (in red), compared to our experimentally calculated difference in apparent magnitude between Betelgeuse and Bellatrix inverted.  Where Day 0 is the 28th Jan.

@Betelbot ‘s data and ours combined, with trend line calculated for both sets of data.

Brief Description of our code:

Importing our .fit files (a file format used in astronomy) into python, there are first and foremost several important points we need to take into account. All the images have noise, some images have hot pixels, and the position of the stars differ for each image. To solve these problems at the same time, Dhruv came up with an idea, to select a picture form a session, use an array to find the position of the maximum pixel brightness, then define a box around the star of fixed width, and only take data from that box. This way, any hot pixels can be cut out, the amount of noise will be relatively the same for the pictures in each session, and as the position of the stars do not move around too much, the star will still be in the box of data.

Note: We will be attaching parts of the code alongside its explanation for future references

After importing all our relevant modules into our code, we must first preview a single Betelgeuse image from the session, looking at the file’s information and then the data available in the file. From this we can extract the minimum, maximum, standard deviation, or mean of the pixel brightness.

Code:

#importing relevant modules
import numpy as np
import os
import math
from PIL import Image
import astropy
%config InlineBackend.rc = {}
import matplotlib
import matplotlib.pyplot as plt
%matplotlib inline
from astropy.io import fits

#opening single image and it's data
from astropy.nddata import Cutout2Dhdulist1a = fits.open(r"C:\Users\dhruv\Desktop\Project Pics\Betelgeuse\Session 1\betelgeuse_750.fit")
hdulist1a.info()
print(repr(hdulist1a[0].header))
data1a = ((hdulist1a[0].data)/256)

print('Min:', np.min(data1a))
print('Max:', np.max(data1a))
print('Mean:', np.mean(data1a))
print('Stdev:', np.std(data1a))

The next part of the code involves finding the position of the maximum pixel brightness i.e. where the star will be located, through an array. From this, we create a 50×50 cut out around the star, and graphed the cut out. Now that we have a preview on an image, we can repeat this easily for the 10-15 other images taken per session, add all their data and find the mean average of their data.

Code:

#locating position of star and creating a cut out
position1 = (223, 248)
size1 = (50,50)
cutout1a = Cutout2D(np.flipud(data1a), position1, size1)
plt.imshow(np.flipud(data1a), origin='lower', cmap='gray') #using np.flipud as the image is inverted
cutout1a.plot_on_original(color='white') #displacing the position of the cut out in the image
plt.colorbar()

#plotting cut out
plt.imshow(cutout1a.data, origin='lower', cmap='gray')
plt.colorbar()

#importing the rest of the images taken in the session
hdulist2a = fits.open(r"C:\Users\dhruv\Desktop\Project Pics\Betelgeuse\Session 1\betelgeuse_751.fit")
hdulist3a = fits.open(r"C:\Users\dhruv\Desktop\Project Pics\Betelgeuse\Session 1\betelgeuse_752.fit")
hdulist4a = fits.open(r"C:\Users\dhruv\Desktop\Project Pics\Betelgeuse\Session 1\betelgeuse_753.fit")
hdulist5a = fits.open(r"C:\Users\dhruv\Desktop\Project Pics\Betelgeuse\Session 1\betelgeuse_754.fit")
hdulist6a = fits.open(r"C:\Users\dhruv\Desktop\Project Pics\Betelgeuse\Session 1\betelgeuse_755.fit")
hdulist7a = fits.open(r"C:\Users\dhruv\Desktop\Project Pics\Betelgeuse\Session 1\betelgeuse_756.fit")
hdulist8a = fits.open(r"C:\Users\dhruv\Desktop\Project Pics\Betelgeuse\Session 1\betelgeuse_757.fit")
hdulist9a = fits.open(r"C:\Users\dhruv\Desktop\Project Pics\Betelgeuse\Session 1\betelgeuse_758.fit")
hdulist10a = fits.open(r"C:\Users\dhruv\Desktop\Project Pics\Betelgeuse\Session 1\betelgeuse_759.fit")
hdulist11a = fits.open(r"C:\Users\dhruv\Desktop\Project Pics\Betelgeuse\Session 1\betelgeuse_760.fit")
hdulist12a = fits.open(r"C:\Users\dhruv\Desktop\Project Pics\Betelgeuse\Session 1\betelgeuse_761.fit")
hdulist13a = fits.open(r"C:\Users\dhruv\Desktop\Project Pics\Betelgeuse\Session 1\betelgeuse_762.fit")
hdulist14a = fits.open(r"C:\Users\dhruv\Desktop\Project Pics\Betelgeuse\Session 1\betelgeuse_763.fit")

#extracting the data of each image
data2a = (hdulist2a[0].data)/256
data3a = (hdulist3a[0].data)/256
data4a = (hdulist4a[0].data)/256
data5a = (hdulist5a[0].data)/256
data6a = (hdulist6a[0].data)/256
data7a = (hdulist7a[0].data)/256
data8a = (hdulist8a[0].data)/256
data9a = (hdulist9a[0].data)/256
data10a = (hdulist10a[0].data)/256
data11a = (hdulist11a[0].data)/256
data12a = (hdulist12a[0].data)/256
data13a = (hdulist13a[0].data)/256
data14a = (hdulist14a[0].data)/256

#stacking the images and plotting
sum1 = data1a+data2a+data3a+data4a+data5a+data6a+data7a+data8a+data9a+data10a+data11a+data12a+data13a+data14a
average1 = sum1/14
ax1a = plt.plot(average1)

All 14 of the images of Betelgeuse taken in session 1 stacked, where the x axis is the number of pixels???, and the y-axis is the relative mean average brightness of?? the stacked images.

The graph above shows the cut out relative to the actual image taken of Betelgeuse in Session 1.

This graph displays the actual cut out, where we shall be focusing on taking our data from.

An interesting thing we can find in our stacked image is the amount of background noise there is, which can be seen if we reduce the y-axis stacked data to a limit of 1.10 to 1.50. Background noise is light registered by the CCD from a space in the sky where there are no light sources. Background noise is often caused by light diffusion or scattering in our atmosphere, and telescopes in space such as the Hubble Space Telescope do not encounter such a problem.

Code:

#changing the y-axis of the graph to show noise
ax1a = plt.plot(average1)
plt.ylim(1.1,1.5)

A zoomed in graph of the earlier stacked images to a y-axis limit of 1.10 to 1.50, displaying the noise present in the images.

This is the sum of the stacked cut out data after converting the image into an array. One thing to point out is that there are to be two peaks, i.e. the position of the star moved. This does not change anything as we are finding the summed maximum pixel brightness.

Code:

#repeating finding the position of the star and creating cut out for the stacked image
position2 = (224, 253)
size2 = (50,50)
cutoutf1 = Cutout2D(np.flipud(average1), position2, size2)
plt.imshow(np.flipud(average1), origin='lower', cmap='gray')
cutoutf1.plot_on_original(color='white')
plt.colorbar()
plt.imshow(cutoutf1.data, origin='lower', cmap='gray')
plt.colorbar()
xf1 = cutoutf1.data
ax1b = plt.plot(xf1)
mn1 = np.sum(xf1)
print(mn1)

print('Min:', np.min(xf1))
print('Max:', np.max(xf1))
print('Mean:', np.mean(xf1))
print('Stdev:', np.std(xf1))

Finally, we convert the stacked images into an array, and sum all the data, as seen in the graph above. The curve we produced is seen to be a nice Gaussian distribution. As we are looking for the ratio of apparent brightness between Bellatrix and Betelgeuse, we do not need the apparent magnitude. The whole process is repeated for Bellatrix, and then the sum of the respective data is placed into a ratio and normalised 2.512 (due to the reverse logarithmic scale of apparent magnitudes), to give the experimental difference in apparent magnitudes.

Code:

#where mn1 is the sum of brightness for Betelgeuse and mn2 for Belletrix 
mn1/mn2
#finding the difference in apparent magnitude
a = math.log(mn1/mn2,2.512)
print a
((a-1.22)/1.22)*100  #percentage error

Graphs of Interest

Discarded stacked data graph of failed session, with a hot pixel seen as a green line on the left.

The same graph with the y-axis reduced to 1.0-3.5 mean photons received in 0.1 s seconds. Still observing the hot pixel as well as an increase in noise past the x-axis 535 pixels mark.

When using the telescope, the CCD (Charged-Coupled device) needs to be cooled, up to around -31ºC. However, the graphs above shows the result of taking images when the CCD is not sufficiently cooled, in this image, the camera was cooled to around -7ºC. Hot pixels occur when the camera sensor heats up, an electrical charge is leaked, resulting in a pixel that is much brighter than expected. In our data, there is an obvious hot pixel, as there is a pixel that is very bright but not located at where the star is, additionally, there is no Gaussian spread in it neighbouring pixels. The hot pixel should go away after the CCD is sufficiently cooled.

Rapidly cooling the CCD may cause ice crystals to form in the detector, often dry ceramic pills placed in the detector to reduce the amount of moisture. On the far right of our data graphs, we see one of the stacked photos have a large amount of noise compared to the other photos. A few plausible explanations for such a phenomenon: ice crystals forming due to the CCD rapidly cooling, or a naughty undergrad using their phone near the telescope, increasing the background luminosity.

Summary

The span of our allocated project time is minuscule compare to many astronomical events. The dimming of Betelgeuse occurred over 4-5 months and should still be monitored in the future due to it sudden brightening. By the time this blog is published, we have only taken 5 session of images of Betelgeuse and Bellatrix, and to produce a trend or a timeline worth looking at, we need many more sessions and images. Clouds remain our worst enemy. Springtime in London brings rain and cloudy days, and we cannot image on a cloudy night.

However, even though Malcolm dislikes the company of optimists, there are still a few more weeks before the end of the allocated project time, where we will continue imaging Betelgeuse and Bellatrix in hopes of providing you with a more complete trend of Betelgeuse’s dimming and brightening.