The Hubble Law

Goals of this project

• To learn how read and display FITS images
• To learn how to analyze imaging data
• To learn how to extract a spectrum from a 2-dimensional image
• To verify that the universe is expanding

You must also have read the IDL primer. This is a data intensive lab and you will not be able to waste much time the first period reading the notes.

Grading Policy

Grading is based on the quality of the data analysis, the discussion of the uncertainties, and the answers to the questions at the end of this page.

In 1920 the National Academy of Sciences held a great debate, between Heber Curtis and Harlow Shapley about the scale of the Milky Way. Curtis argued for a large scale of the Milky Way, with the spiral nebulae as gas clouds within the Milky Way. Shapley argued for a much smaller Milky Way, with the spiral nebulae as "island universes", similar in size to the Milky Way, but very far away. Curtis was correct about the size of the Milky Way, but Shapley's view about the spiral nebulae was ultimately proven correct when in 1925 Edwin Hubble, using the Mt. Wilson 100 inch telescope, resolved the outer parts of M31, the Andromeda Nebula, into stars. The spiral nebulae are galaxies much like our Milky Way.

In the 1930's Edwin Hubble discovered and quantified a relationship between the distance of galaxies and their recessional velocities.
Left: Hubble guiding the Palomar Observatory 48" Schmidt telescope. This is a publicity shot: it is generally not recommended that one smoke while observing.

Solutions of the equations of general relativity show that a static universe is stable provided it is empty, but a universe filled with matter cannot be static. This led Albert Einstein to formulate the Cosmological Constant. Hubble's observations provided the first proof that the universe was not static, but was expanding, and led Einstein to renounce the Cosmological Constant (it is now back, but for other reasons).

The inverse of H₀is proportional to the age of the universe. Measuring the value of the Hubble Constant has not proved simple; the currently accepted value is about a factor of 7 smaller than Hubble's first estimate.

Your task in this lab is to verify that Hubble got it right. We have images and spectra of a sample of galaxies. Your task is to determine the redshifts, and to estimate the relative distances.

• Hubble, E. The Realm of the Nebulae, 1936, Yale University Press, reprinted 1958, Dover Publications

Notes:
S: Seyfert 1 galaxy. Nucleus may be abnormally bright.
X: Possibly tidally truncated. May give misleading results.

Coordinates are epoch 2000. All have nominal Hubble types of Sc.

The data are stored on two CDs which should be easy to locate in the lab room. The data for each galaxy are in separate directories. Refer to the Table above to locate a specific data set.

The archival data include:

• Images of the galaxies (*_icube*.fits). You will use these to determine the relative distances.
• Extracted spectra (*_s_gx*.fits, where x is 26 or 47.). Use these to determine the redshifts, unless you prefer to do your own extractions.
• Two dimensional spectra (*_sp_x.fits, where x is r or b). You may use these to extract your own spectra.
• A list of stars in the field (stars.dat).

All the images and spectra are stored in FITS format.

II A. The Images

The images were obtained using the SMARTS 0.9m CCD imager at the Cerro Tololo Interamerican Observatory during the spring of 2003. You should read up on CCDs. The data have been bias-corrected, trimmed, and flat-fielded.

If there is a single image cube, the name is GAL_icube.fits (where GAL is the name of the galaxy); if there are two images cubes, the names are GAL_icube*.fits, where * is either 1 or 2. The two images were taken on different days, and one is generally much longer (deeper) than the other. In general, use the deeper exposure.

Use IDL to read in the images. For example, to read the I band image for NGC2525, use the following sequence of commands:

cd,'NGC2525'
d=readfits('NGC2525_icube.fits',h)
print,sxpar(h,'filter*')     ; print header words filter1, filter 2, etc.
                             ; result is "B        V        R        I"
image=d(*,*,3)

To automate the image analysis, you could find and extract the I band image automatically with the following two lines:

k=where(strtrim(sxpar(h,'filter*'),2) eq 'I')
image=d(*,*,k)

Display the image using the TV command. The TV command is primitive: it does not scale the data, or resize the TV window. If the image is larger than the size of the monitor, you should probably rebin it (see this section) to a manageable size. Use the TVPLOT command or, for more control, size the plot window with the XSIZE and YSIZE keywords and do the flux scaling manually (see the "Using the Cursor to Measure Positions" section in the IDL primer). The images contain some bad columns, which can take very high or very low values. Some images contain very bright stars, which can make it hard to find a low surface brightness galaxy. I suggest you find the sky level, which to a good approximation is the median level in the image. In IDL, type m=median(image) to determine the median value m.

Plotting image>m will set the floor of the image to the median, getting rid of all negative points. To truncate high points, try plotting (image>m)<(m+n), where n is some threshhold level of counts above the median. You should also try plotting the log of the image (make sure there are no negative of zero values). Experiment until you find a good visualization of the galaxy. You could also experiment with plotting the logarithm of the image to bring up the fainter details while retaining dynamic range.

Note that some of the images are rather large (the original images are 2048 x2048 pixels). Feel free to trim the images as appropriate. For example, to make an 800 x 800 image out of a 1200x1200 image, retaining the same image center, you would type IMAGE=IMAGE(200:999,200:999)

Note that the keywords RA and DEC correspond to the actual pointing position of the telescope. You need to precess this from the epoch of the observation to epoch 2000, and then figure out where the central pixel is.

When you convert seconds of time to seconds of arc, remember that there is a cos(dec) factor in the conversion.

II B. The Spectra

The images were obtained using the SMARTS 1.5m RC spectrograph at the Cerro Tololo Interamerican Observatory, mostly during the spring of 2003.

The two dimensional spectral images have been bias-corrected, trimmed, and flat-fielded. Three spectra were obtained of each object. These three spectra were median-filtered to remove cosmic rays. A wavelength comparison spectrum of a hollow-cathode arc lamp was obtained along with each spectrum. These have been reduced to produce the wavelength solution. The solution in a 4^th or 5^th order polynomial; the coefficients are in the FITS header.

The slit runs east-west; the slit width is 110 microns (1.5 arcsec).

• 0: the net extracted spectrum from a boxcar extraction
• 1: the background corresponding to the extracted spectrum
• 2: the uncertainty on the net extracted spectrum
• 3: the uncertainty on the background spectrum
• 4: the net extracted spectrum from a gaussian extraction
• 5: the background corresponding to the gaussian spectrum
• 6: the counts to flux conversion for spectra 0 - 3
• 7: the counts to flux conversion for spectra 4 and 5

Extract the spectra as spectra = readfits(file,header)
where file is the name of the fits file, and spectra and header are the output data and header vectors. spectra will be an 1199 x 8 array.

The spectra have been linearized in wavelength; the wavelength of the first pixel and the number of Angstroms per pixel are given in the header keywords wave_0 and delta_w. Construct the wavelenth vector as w=sxpar(header,'wave_0')+sxpar(header,'delta_w')*findgen(1199)

Typical spatially-resolved spectral images (actually of the brightest of the galaxies, M83), are shown and described on this page.

The spectral resolution is 3.1 Angstroms.

The spectral resolution in the blue spectra is 4.3 Angstroms.

II C. The stars.dat File

The stars.dat file is the output of a search of the USNO astrometric catalog, centered on the nominal position of the galaxy nucleus. Note that the positions given in the headers of the fits images are the positions of the telescope at the epoch of observation (usually in 2003), and may be off due to telescope drift.

The file provides the coordinates and magnitudes for stars within a square of size 7 arcmin centered on the galaxy. The columns headed Bn, Rn, I2 are, respectively, the B, R, and I magnitudes measured at two epochs (n=1, 2). An entry of 0.0 means no magnitude was measured. The B, R, and I magnitudes correspond approximately (but not exactly) to the B, R, and I bands in our images. See the USNO astrometric catalog help page more detailed descriptions. I have retained only stars with I2 magnitudes <18.

The last entry in each line is the distance from the center. Be leery of anything within a few arcsec of the galaxy: the nucleus is sometimes picked up as a star. Also, nebulosities and star clusters in the galaxy may be listed here as stars.

Because the image sizes vary, some stars may be outside the boundaries of the image.

II D. Which Data Should You Use?

We recommend the following:

• If there are two image cubes, use the deeper images.
• Select a single filter band to examine for all the galaxies. There is photometry available in the B, R, and I bands. HII regions and the spiral arms stand out best at the shorter wavelengths (B,V), while the smoother old disk population dominates at I.
• If you want a challenge, extract your own spectra from the two-dimensional spectral images. You will get brownie points for doing this (but only if you do it correctly). But since you only have 3-4 weeks for this lab, most of you are probably better advised to use the extracted spectra for the redshift determinations.
• The red spectra are easier to analyze than the blue spectra, but the blue spectra are more akin to the type of data Slipher and Hubble had to work with (photographic plates are blue-sensitive). Use the blue spectra if you want a challenge.

II A. Images

You want to measure some distance-dependent quantity. For example, you could measure the angular size of the galaxy, either in pixels or arcsec, using the XPIXSIZE and YPIXSIZE header words to convert. The angular size should scale linearly with distance.

You could measure the brightness of the nucleus. For this you will need to do an aperture-extraction to determine the total number of counts in the nucleus. Be careful with the background subtraction.

You could measure the brightness of the brightest nebulosities in the B or V band images. As in the Olympics, throw away the brightest few (they may be pathological), and compare the 4^th or 5^th brightest. Again, be careful with the background subtraction.

You can select other criteria, but you must justify them. You should discuss them with the instructor first.

Note that you will not be doing rigorous photometry here. A brief introduction to everything that photometry entails is presented in this AST 443 photometry page.

The simplest way is to extract the counts in a square region centered on some object of interest at coordinates xc, yc For a square extraction region of radius r, the total counts are cts = image(xc-r:xc+r,yc-r:yc+r). You extract a background region in a similar way, avoiding stars and galaxy light. The background region need not be the same are as the extraction area, but if it is not you must scale the background to the area in the extraction aperture. Then subtract this background off to get the net counts.

The background should be taken near to, and preferably surrounding, the extraction region, in case there are gradients in the background. You can extract background in a larger area centered on the extraction box, or you can average a number of background regions. The most important thing is to explain carefully. what you did.

A more complex extraction would involve a circular aperture and concentric annuli for the background.

Always remember to propagate your errors correctly.

If you choose to measure the brightness of something, you will have to correct to absolute. To do this you need a reference of known brightness. This is where the stars.dat file comes in.

stars.dat is a printable text file. You can read the values into IDL variables using the rd_stars.pro procedure. The calling sequence is rd_stars,stars. The output variable stars is an IDL structure containing the data. Type rd_stars,/help for online help.

The structure contains the right ascension (stars.ra) and declination (stars.dec) of the stars. Identify a star in the image by taking the diffence in position between the star and the center of the galaxy in arcsec, and then convert to pixels. North is up in the image, East is left. Right ascension increases towards the east. The pixel scale is 0.401 arcsec per pixel. The position of the galaxy is given in the table above, or in the image header as keywords CRVAL2 and CRVAL2 (see the Astrometry section).

After you identify some stars in the image, you will perform aperture photometry on them to determine the number of detected counts. Convert the counts to the instrumental magnitude: mag = -2.5 log(counts). The difference betwen the instrumental and apparent magnitudes is the magnitude conversion factor. It is different for each filter, and varies with time due to weather condition, but should be the same everywhere within a given image. Use at least 3 stars to determine the conversion in order to reduce the scatter.

III B. Spectra

• First, extract the spectra, and then
• Measure the redshift

1. The data are in fits format. Use the READFITS function to extract the data.
2. Plot the image using, for example, the TVPLOT procedure. The nucleus of the galaxy is the bright horizontal line near the middle of the image
3. Locate the galaxy by performing a cross cut. Reset the plot scaling with SETXY and SVP. Extract a horizontal cut, e.g., CUT=IMAGE(X,*), where X is the X value at which the cut is made. To improve the S/N, sum over a range of X pixels with CUT=total(IMAGE(X1:X2,*),1).
   Then use the cursor to mark the peak, e.g., CURSOR,xc,yc & print,xc. Move the cursor to the peak position and click the left mouse button.
4. To verify this position, you can mark the position by typing OPLOT,[xc,xc],!y.crange Determine the width of the spectrum using the cursor. Note that the spectrum is not exactly horizontal. Verify this by overplotting a cut through the image near X=100 and near X=1100.
   Ideally you would write a program which would determine the center of the spectrum at each X pixel, and then extract the data centered on that central Y(X). However, it won't cost much in terms of S/N to simply extract a linear spectrum, so long as it is broad enough to contain all the light. In this case, the spectrum is given by SPECTRUM=TOTAL(D(*,Y1:Y2),2), where Y1 and Y2 are the lower and upper extent of the extraction box.
5. There is background from the night sky. In the red spectra much of this is in the form of emission lines; in the blue much is scattered continuum. Subtract the background. Define extraction boxes above and below the spectrum. Be sure that the area in the background region extracted is the same as the spectral region. You might want to average the bacgrounds above and below the nucleus. Then, the background-subtracted spectrum is SPECTRUM = SPECTRUM - BACKGROUND

   Note that the galaxies are extended on the sky, and there is galactic starlight in many cases. You may be able to see this in the vertical cuts. You may want to avoid this, if possible, but it is not crucial if all you intend to do is measure the radial velocities.

6. Establish the wavelength scale. The X axis (in pixels) is given by X = FINDGEN(1199). The coefficients of the polynomial wavelength solution are given by the header keywords WAVE_0, WAVE_1, and W_COEFFS. The wavelength scale is given by W=WAVE_0 + WAVE_1 * X + W_COEFFS(0) * X² + W_COEFFS(1) * X³ + ... .
7. You may use this wavelength scale, or you may linearize it using the INTERPOL function.
8. Then PLOT,W,SPECTRUM. Identify known features, perhaps using the CURSOR command.
9. Determine the redshift, Z = (Wmeas-W0)/W0), where Wmeas is the observd wavelength and W0 is the rest frame (laboratory) wavelength. and the radial valocity, V = cZ.
10. You should measure more than one feature, and average results to reduce the measurement errors.

The spectra are easier to measure and comprehend than the distances. To ensure you get a broad range of recession velocities, I recommend you first determine redshifts for ALL the galaxies. Then select a sample distributed over the full range of recessional velocities for the distance determinations.

You should also sample a number of galaxies with similar recessional velocities in order to investigate your scatter at a given distance.

Comment on the quality of your results. Why is there so much scatter among the closer galaxies? Do all galaxies show a redshift (positive velocity). How confident are you of your relative distances?

In order to determine the value for the Hubble constant, you will need to know the distance to one of the galaxies. Here are the actual distances (based on the Cepheid period-luminosity relation) in Mpc to three galaxies:

Question: Why is the mean surface brightness of the galaxy in the I band not a good distance indicator?

Return to top

Return to PHY 445/515 Astronomy page