ÐÏࡱá>þÿ /£¥þÿÿÿ“”•–—˜™š›œžŸ ¡¢¨©ª«¬­®¯°±²³´µëlínïpñrótõv÷eælÌÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿì¥Áq` ø¿<bjbjqPqP 3::Ñx;ÿÿÿÿÿÿ¤ÎÎÎÎê î ö "Äšî»î»î»8&¼Ôú¾¤æ¶L”hªÄ¢LÆ\Æ\Æ\Æö7÷T‹÷,ی݌݌݌݌݌݌$´•h˜‚ê ; óö;;ÎÎ\Æ\Æ­,”+1+1+1;4!Îö\Æê \ÆÛŒ+1;ÛŒ+1+1¾yrøÄ &ê Iv\ÆžÄ ÈdÇî»o#@ qt8€Ì ”0L”©t ž˜¯/`ž˜pIvIvªž˜î ów·÷Nû\+1aýäEÿö·÷·÷·÷1·÷·÷·÷L”;;;;æ¶æ¶æ¶D?*öÄÅæ¶æ¶æ¶*ö ä î0ÎÎÎÎÎÎÿÿÿÿ  Improving Radial Velocity Analysis at Lick Observatory A thesis submitted to the faculty of San Francisco State University in partial fulfillment of the requirements for the degree Master of Science in Physics by David A. Abouav San Francisco, California January, 2007 CERTIFICATION OF APPROVAL I certify that I have read Improving Radial Velocity Analysis at Lick Observatory by David A. Abouav, and that in my opinion this work meets the criteria for approving a thesis submitted in partial fulfillment of the requirements for the degree: Master of Science in Physics at San Francisco State University. _________________________________ Debra A. Fischer Professor of Astronomy _________________________________ Geoff W. Marcy Professor of Astronomy _________________________________ Ron Marzke Professor of Astronomy Improving Radial Velocity Analysis at Lick Observatory David A. Abouav San Francisco State University 2007 The mutual gravitational attraction between a star and an orbiting planet will cause the star to wobble about their common center of mass. By measuring the star's varying velocity over time, and applying Kepler's laws of motion and Newton's law of gravity, it is possible to infer the characteristics of the planet's orbit. Measuring these velocities is accomplished by modeling the Doppler shift in spectral absorption lines. This thesis describes modifications to the SFSU Doppler analysis code, and the resulting improvements in radial velocity analysis. Several years of existing data from Lick Observatory were reanalyzed and examined to search for previously undetected planets and to improve the Keplerian orbital parameters of known planets. I certify that the Abstract is a correct representation of the content of this thesis. ___________________________________________ _______________ Chair, Thesis Committee Date ACKNOWLEDGEMENTS Dr. Debra Fischer, my thesis advisor, for her continual support, guidance, and enthusiasm as I made my way through this project (especially during our energy-filled weekly meetings). Also, for giving me a chance to practice experimental science, which has benefited my education as much, or more, than learning theory and solving problems analytically. Jenny Abouav, my wife, for her love, support, and patience over the last three years as I pursued my dream of earning a degree in physics; a dream which she encouraged me to follow, despite my doubts. Mr. Art Fortgang, my high school physics teacher, for doing such an excellent job in teaching his Honors and AP classes. These classes gave me a solid foundation in physics that I still use when solving problems. Mr. Fortgang also helped me explore my interest in teaching, by finding me students to tutor. iv TABLE OF CONTENTS List of Tables………...…………..……………………………………………………...viii List of Figures……….…………………………………………………………………….x Chapter 1 – Background…………………………………………………………………..1 Section 1.1: “Looking” for Planets…………………………………………......…1 Section 1.2: Stellar Wobble……………………………………………………….4 Section 1.3: Taking Observations…………………………………………………6 Section 1.4: Studying Starlight……………………………………………………8 Section 1.5: Doppler Analysis – Theory…………………………………………10 Section 1.6: Doppler Analysis – Practice………………………………………..13 Section 1.7: PSF: The Spreading of Light……………………………………….14 Chapter 2 – Doppler Precision of 3 m s-1………………………………………………...18 Chapter 3 – The Doppler Analysis Computer Program………………………………….21 Section 3.1: Chunks and Observations…………………………………………..21 Section 3.2: Wavelength Solution………………………………………………..24 Section 3.3: IPCF and vdiods…………………………………………………....25 Section 3.4: The Stellar Template (DSST)………………………………………26 Section 3.5: The PSF…………………………………………………………….29 Section 3.6: (r2 Fits………………………………………………………………30 Section 3.7: The Marquardt Minimization Algorithm…………………………...32 Section 3.8: Frozen, Averaged PSF……………………………………………...35 Section 3.9: Velocity Analysis Code…………………………………………….38 v Section 3.10: Chunk-Sets and Weights…………………………………………..40 Section 3.11: Mean Velocity, and Errors………………………………………...49 Chapter 4 – New Iodine Atlas……………………………………………………………51 Section 4.1: A Baseline Test……………………………………………………..53 Section 4.2: First Test with the New Atlas……………………………………...55 Section 4.3: Fine Tuning the New Atlas…………………………………………58 Section 4.4: Further Tests with the New Atlas…………………………………..59 Chapter 5 – Polynomial Fitting, PSF Modifications, and IPCF Entries…………………66 Section 5.1: Polynomial Fitting………………………………………………….66 Section 5.2: PSF Modifications – First Pass……………………………………..68 Section 5.3: IPCF Entries………………………………………………………...74 Section 5.4: PSF Modifications – Second Pass………………………………….76 Chapter 6 – New Code and New Observations from UC Berkeley……………………...81 Section 6.1: EMU Observations…………………………………………………81 Section 6.2: Updated Code………………………………………………………82 Chapter 7 – Wavelength Continuity……………………………………………………..86 Section 7.1: Lack of Wavelength Continuity…………………………………….86 Section 7.2: Visualizing the Gaps……………………………………………..…88 Section 7.3: Freezing the Wavelength Solution………………………………….91 Section 7.4: New Baseline Tests…………………………………………………93 Section 7.5: Methods to Enforce Wavelength Continuity……………………….95 vi Section 7.6: Results of Applying Wavelength Continuity……………………….96 Section 7.7: Effect on the Wavelength Solution…………………………………98 Section 7.8: Investigation of Reduced Internal Error (rInt)…………………….101 Section 7.9: Investigation of Increased RMS Scatter (rBin)……………………107 Chapter 8 - Alternate Weighting Schemes……………………………………………...118 Chapter 9 – Marquardt Minimization Method& & & & & & & & & & & & & & & & .121 Chapter 10  IPCF Entries& & & & & & & & & & & & & & & & & & & & & & & ...124 Chapter 11  Results With Final Improvements& & & & & & & & & & & & & & & ..128 Section 11.1: Revisiting the 16 Ä Ceti Observations& & & & & & & & & & & 128 Section 11.2: A Large Set of Ä Ceti Observations& & & & & & & & & & & ...137 Section 11.3: Observations of 51 Pegasi& & & & & & & & & & & & & & & .138 Section 11.4: Analysis of Results with the Old Iodine Atlas& & & & & & & ...140 Chapter 12  Conclusion& & & & & & & & & & & & & & & & & & & & & & & & ..142 References& & & & & & & & & & & & & & & & & & & & …………………………146 vii LIST OF TABLES Table Page 2.1 – Radial velocity amplitudes induced on the Sun……………………………………18 4.1 – Comparison of rInt and rBin using the different iodine atlases……………………56 4.2 – Wave number and dispersion values for the ftslick05 iodine atlas………………...58 4.3 – Results of wave number and dispersion changes…………………………………..59 4.4 – HD217014 analyses done with different atlases (dewar 8)………………………..61 4.5 – HD217014 analyses done with different atlases (dewar 13)………………………62 4.6 – HD217014 analyses done with different atlases (dewar 6)………………………..62 5.1 – Results of changing pord from 3 to 4………………………………………………67 6.1 – Results of using newer Doppler analysis code obtained from UCB……………….84 7.1 – Comparison of wavelength solution in neighboring chunks……………………….86 7.2 – Gap information from baseline test ………………………………………………..90 7.3 – Gap information test with the new code (frozen dispersion)& & & & & & & & & 93 7.4  100 EMU observations analyzed with enforced wavelength continuity& & & & ..97 7.5  25 Ä Ceti observations analyzed with enforced wavelength continuity& & & & & 97 7.6  Comparison of rBin to S * rMed& & & & & & & & & & & & & & & & & & & .116 9.1  Results of altering with the Z step and (0 step for Ä Ceti observations& & & & ...122 9.2  Results of altering (0 step for 100 EMU observations& & & & & & & & & & & 123 11.1  Results for 16 Ä Ceti observations analyzed with different tests& & & & & & ..130 11.2  Results for 16 Ä Ceti observations analyzed with different tests& & & & & & ..131 11.3  Results for 16 Ä Ceti observations analyzed with different tests& & & & & & ..132 11.4  Results for 16 Ä Ceti observations analyzed with different tests& & & & & & ..134 viii 11.5  Comparison of results with and without wavelength continuity& & & & & & ..137 11.6  Comparison of test results for 628 Ä Ceti observations& & & & & & & & & & 138 11.7  Comparison of test results for HD 217014 observations& & & & & & & & & ..139 ix LIST OF FIGURES Figure Page  SEQ Figure \* ARABIC 1.1 - Image of Jupiter and one of its moons, Io…………………………………………...2 1. SEQ Figure \* ARABIC 2 - Velocity versus time plot for the star 51 Pegasi……………………………………..5 1.3 – Spectrum of Sun produced by an Echelle Spectrometer…………………………….8 1.4 - Spectrum for order 45 at Lick Observatory, 51 Pegasi obsv………………………...9 1.5 – Magnified region of Figure 1.4…………………………………………………….10 1.6 – Sample Point Spread Function (PSF)……………………………………………...15 1.7 – Application of sample PSF to a spectrum………………………………………….16 3.1 – Flow of the Doppler analysis program for a single chunk…………………………23 3.2 – Illustration of the stellar Doppler shift……………………………………………..28 3.3 – Example of a function with local and global minima……………………………...33 3.4 – Illustration showing “neighboring” chunks in an obsv…………………………….37 3.5 – Sample output of the vank program………………………………………………..39 3.6 – Illustration of a chunk-set made from four vd files& & & & & & & & & & & & ..41 3.7  Conceptual illustrations of the diff and sigma arrays& & & & & & & & & & & ...44 4.1 - Overlay of the new and old iodine atlases& & & & & & & & & & & & & & & & .52 4.2  Vank plot from the baseline test of 16 Ä Ceti obsvs& & & & & & & & & & & & ..54 4.3  Vank plot from a test of 16 Ä Ceti obsvs. with the new atlas& & & & & & & & & 55 4.4  RV plot for obsvs. of the star 51 Pegasi, old atlas& & & & & & & & & & & & & 60 4.5  RV plot for obsvs. of the star 51 Pegasi, new atlas& & & & & & & & & & & & ...61 4.6 – Over plotted (r2 fits for the chunks of two vdiod files……………………………..63 5.1 – PSF FWHM contour plot for a vd file from the baseline test……………………...70 5.2 – Plot of the averaged PSF for the chunk with the median (r2 fit……………………70 x 5.3 – Vank plot resulting from application of PSF A……………………………………72 5.4 – Plot of the averaged PSF for the chunk with the median (r2 fit……………………73 5.5 – PSF FWHM contour plot for a vd file creating using PSF A……………………...73 5.6 – Vank plot resulting from application of PSF RF-L………………………………..77 5.7 – Plot of the averaged PSF for the chunk with the median (r2 fit……………………77 5.8 – PSF FWHM contour plot for a vd file creating using PSF RF-L………………….78 6.1 – Vank plot resulting from use of new code from UCB……………………………..84 7.1 – Sample output of the view_gaps program…………………………………………89 7.2 – Sample output of the view_gaps program…………………………………………90 7.3 – Output of the view_gaps program, frozen dispersion invoked…………………….92 7.4 – Vank plot for 100 EMU obsvs., no wavelength continuity…………& & & & & ..94 7.5  Vank plot for 25 Ä Ceti obsvs., no wavelength continuity& & & & & & & & & & 94 7.6  Vank plot for100 EMU obsvs., wavelength continuity& & & & & & & & & & & 96 7.7  Vank plot for 25 Ä Ceti obsvs., wavelength continuity& & & & & & & & & & & .97 7.8  Histogram from the view_gaps program…………………………………………..99 7.9 – Comparison to “actual” wavelength solution for EMU obsv…………………….100 7.10 – Chunk velocities for EMU obsv., no wavelength continuity…………………....102 7.11 – Chunk velocities, EMU obsv., wavelength continuity………………………….102 7.12 – Chunk velocities, Ä Ceti obsv., no wavelength continuity& & & & & & & & & 104 7.13  Chunk velocities, Ä Ceti obsv., wavelength continuity& & & & & & & & & & .104 7.14  1/sigma2 vs. chunk-set, Ä Ceti obsv., no wavelength continuity& & & & & & ...106 7.15  1/sigma2 vs. chunk-set, Ä Ceti obsv., wavelength continuity……………………106 7.16 – 1/sigma2 vs. median of diff entries, no wavelength continuity…………………109 7.17 – 1/sigma2 vs. median of diff entries, wavelength continuity…………………….110 xi 7.18 – Mean vs. median velocity, EMU obsvs., no wavelength continuity…………….112 7.19 – Mean vs. median velocity, Ä Ceti obsvs., no wavelength continuity& & & & & 112 7.20  Mean vs. median velocity, EMU obsvs., wavelength continuity& & & & & & ..113 7.21  Mean vs. median velocity, Ä Ceti obsvs., wavelength continuity& & & & & & .113 7.22  Mean vs. median velocity, 51 Peg. Obsvs., wavelength continuity& & & & & ..114 8.1  Chunk velocity vs. (r2 for obsv. of Ä Ceti, no wavelength continuity& & & & & .118 10.1  Averaged PSF for Ä Ceti obsv. processed with new IPCF file& & & & & & & .125 10.2  PSF FWHM plot for Ä Ceti obsv. processed with new IPCF file& & & & & & .125 11.1  Vank plot from the baseline test of 16 Ä Ceti obsvs& & & & & & & & & & & ..129 11.2  PSF FWHM plot for Ä Ceti obsv. analyzed with PSF RF-L& & & & & & & & ..132 11.3 - Vank plot: new code, newer IPCF, PSF RF-L& & & & & & & & & & & & & ..133 11.4 - Vank plot: new code, newer IPCF, PSF RF-L, wavelength continuity& & & & .135 11.5  PSF FWHM plot for Ä Ceti obsv. with all changes adopted& & & & & & & & .136 11.6  RV plot for 51 Pegasi, analyzed using wavelength continuity& & & & & & & .139 xii Chapter 1 - Background Our Milky Way galaxy is a collection of billions of stars, our Sun being just one of them. For millennium human beings have gazed at the night sky, looking at the small fraction of these stars visible to us, and have pondered a simple question: could some of these stars host other worlds? An answer to that question was finally provided in 1995, when Michel Mayor and Didier Queloz at Geneva Observatory announced that they had discovered a planet orbiting the normal star 51 Pegasi (Mayor 1995). Since their discovery, over 200 planets orbiting stars in our galaxy have been discovered, many of which have been detected though the work of the California Carnegie Planet Search Team (http://www.exoplanets.org). All of the planets discovered to date have been detected using the Doppler analysis technique. This chapter presents some background knowledge in observational astronomy, and on the Doppler analysis technique. A reader who is already familiar with these subjects may skip ahead to Chapter 2. Section 1.1: “Looking” for Planets The telescope was first developed in the early 17th century by Dutch astronomers, and was subsequently improved by Galileo Galilei, who was the first to turn it towards the 1 2 night sky. Since then human beings have used telescopes to study thousands of objects in the night sky. Some of the first objects studied in detail were the planets of our solar system. Modern telescopes allow us to study the planets, and their moons, in far greater detail than during Galileo’s time. Figure 1.1 shows an image of the planet Jupiter, along with one of its moons Io, taken by the Hubble Space Telescope. Galileo was the first person to observe Jupiter’s moons, which through his telescope appeared only as small points of light.  Figure 1.1 - Image of Jupiter and one of its moons, Io. The dark circle to the right of Io is its shadow. Using his telescope, Galileo was the first person to observe that Jupiter has moons. Image courtesy of hubblesite.org. 3 In 1995, astronomers began a new chapter in astronomy by studying planets outside of our solar system. These exoplanets are planets that do not orbit our Sun, but rather, orbit the stars we see in the night sky. But unlike the planets of our solar system, these distant objects cannot be seen through a telescope. The reason for this is two-fold. The first and most fundamental problem is simply one of distance. The further away you are from two nearby objects, the closer they will appear relative to your great distance. Imagine observing the Earth and Sun from very far away. The Sun would look like a small glowing ball and the Earth like a blue point of light nearby it. As you increased your distance, the Sun would shrink in size, and the two objects would appear to be closer together in your field of view. Finally, from a very great distance, the Sun and the Earth would appear to you to occupy the same point. Such is the case with a distant star and its planets. Resolving a planet from a star at our great distance is a difficult process. Though some solutions have been explored, none have yet been successful in finding an exoplanet. The next problem is that planets are not very bright compared to the stars they orbit (host stars). Unlike stars, planets give off no visible light of their own. All of the light visible from a planet is that which is reflected from its parent star, and is much dimmer than the light of the star. For example, the intensity ratio of the light from the Sun to the reflected light from Jupiter is about one-billion to one. Planets do give off light of their own in the 4 infrared spectrum, but this light is also much dimmer than that of the host star (ratios are on the order of one-million to one). Despite these difficult problems, over 200 planets discovered outside of our solar system since 1995. However, none of these planets has been observed directly with a telescope. Instead these planets have been detected by an indirect method, which is described in the following sections. Section 1.2: Stellar Wobble Imagine tying a long piece of string to a tennis ball, and then using the string to swing the ball in a circle above your head. As you swing the ball around and around, your hand will also move around in a circle, though a much smaller one. Your hand wobbles about due to the tension in the string, which pulls on your hand as it keeps the ball moving in a circular path. A phenomenon similar to this occurs on much grander scale, between a star and a planet. In this case, the force of gravity takes the role of the string from our example above, providing an attractive force that keeps a planet in orbit around its star. This force also pulls on the star, causing it to wobble in response, like your hand in the previous example. If we are lucky, from our vantage point on Earth we will be able to detect this wobble as a change in the star’s velocity towards us or away from us over time. 5 Figure 1.2 shows a sample of this velocity variation for the star 51 Pegasi, which is the first star around which a planet was discovered.  Figure 1.2 - Velocity versus Time plot for the star 51 Pegasi. Each dot is a single velocity measurement, obtained from a single observation of the star. Each dot in Figure 1.2 refers to a single velocity measurement, obtained from a single observation of the star. The curve passing through the dots is an extrapolation of how the star’s velocity varies with time. Calculating such a curve to requires the application Newton’s law of gravitation, and Kepler’s laws of planetary motion. Once a best curve is 6 achieved, astronomers are able to infer the characteristic of the orbiting planet. Specifically, astronomers can determine the planet’s mass, the period of its orbit around the star, and characteristics of its orbit such as its eccentricity. According to Figure 1.2, the star 51 Pegasi wobbles with a period of about 4.23 days (4 days, 6 hours) in response to the planet orbiting around it. The data in  REF _Ref138763294 \* MERGEFORMAT Figure 1.1.2 shows 8 such periods, for a planet with a mass close to half that of Jupiter. Before a curve can be determined, many velocity measurements are required, which means that many observations must be taken. Section 1.3: Taking Observations The California Carnegie Planet Search Team takes observations at several telescopes around the world. The observation sites are listed below: Lick Observatory, Mount Hamilton, California Subaru Telescope, Mauna Kea, Hawaii Keck Observatory, Mauna Kea, Hawaii Las Campanas Observatory, Chilean Andes Anglo-Australian Observatory (AAO), New South Wales, Australia We will focus our discussion on the Lick Observatory, though the process of taking stellar observations and performing Doppler analysis is similar for all observatories. 7 The telescopes at Lick Observatory use optical elements, such as lenses and mirrors, to guide and focus incoming light. The light that travels through the telescope ultimately falls onto a device called a CCD, which records it. A CCD is a fine grid of square regions called pixels, each of which is capable of reporting the intensity of light that falls on it. Depending on which telescope at Lick is used, and the brightness of the star being observed, the exposure time can last from a few minutes, to thirty minutes. Since the Earth constantly rotates, the telescope must also move during this time so that it always points at the star it is observing. Because stars are so far away, their images cannot be magnified. No lens can make a star appear any larger than a single point of light. Thus a photograph of a star taken with even the most powerful telescope would not provide astronomers with any useful data. To do useful astronomy, astronomers instead examine the composition of the star’s light. Starlight, like the light from our Sun, is made up of many different colors. Most of this light is not visible to the human eye, so it makes little sense to refer to it by its “color”. But every color of light, whether our eyes can see it or not, has a specific frequency, or wavelength. By separating the light from a star into its different wavelengths, astronomers can learn a great deal about the star. This separation is accomplished by use of a spectrometer. 8 Section 1.4: Studying Starlight A spectrometer is a device that allows the intensity of light to be measured as a function of the light’s wavelength. During an observation, light from a star travels through the optics of the telescope, and ultimately passes through an Echelle Spectrometer. An Echelle Spectrometer separates the light into its component wavelengths (much like a prism creates a rainbow), and then projects the resulting spectrum onto the telescope’s CCD. In particular, the Echelle Spectrometer uses optical elements such as mirrors and prisms to divide the spectrum into evenly spaced, nearly parallel rows called orders. This allows the entire spectrum of starlight to be projected onto the surface of the CCD. The effect looks something like Figure 1.3.  Figure 1.3 – An Echelle Spectrometer separates light into its different wavelengths, and uses optical elements to divide the spectrum into evenly spaced, nearly parallel rows called orders. This figure shows the result when light from our Sun is passed through an Echelle Spectrometer. The black areas are gaps in the spectrum of our Sun. Image courtesy of the National Optical Astronomy Observatory, N.A.Sharp, NOAO/NSO/Kitt Peak FTS/AURA/NSF. 9 Because the spectrum is spread out over the CCD, each pixel on the CCD is only exposed to light from particular wavelength range. At Lick observatory, for example, each pixel on the current CCD is exposed to a wavelength range of about 0.5 nm. Over the course of the observation, the CCD keeps track of the number of photons that strike each of its pixels. These values are recorded into a computer at the end of the observation, thus providing a record of the intensity as a function of wavelength. Figure 1.4 shows a plot of these values for a particular order in an observation of the star 51 Pegasi.  Figure 1.4 - Spectrum for Order 45 of an observation of the star 51 Pegasi, taken at Lick Observatory. 10  EMBED PBrush  Figure 1.5 – A magnified region of Figure 1.4. Section 1.5: Doppler Analysis – Theory As mentioned previously, by gathering multiple velocity measurements of a star over time and analyzing them, astronomers are able to detect new planets. Measuring a star’s velocity from a stellar observation is done using a technique known as Doppler analysis. Doppler analysis is the study of Doppler shift, which is actually a phenomenon familiar to most people. Though a Doppler shift can occur with any type of wave, generally 11 people are most familiar with Doppler shifts in sound waves. For example, imagine you are standing on the sidewalk, and person is driving a car down the street towards you, all the while holding down the horn. The faster the person drives towards you, the higher pitched the horn will sound. This is due to the fact that as the car races down the street, it is moving at a speed that is comparable to the speed of the sound it is producing. As a result, the sound waves produced by the horn get compressed, and so you’ll hear them at a higher frequency. Conversely, if the car is racing away from you, you hear the horn at a lower pitch.  REF _Ref139281264 \h Equation 1 describes this phenomenon mathematically: v / s = ”» / » Equation  SEQ Equation \* ARABIC 1: An equation to calculate Doppler shift. In  REF _Ref139281264 \h Equation 1, v is the velocity of the object producing the sound (e.g. the car), and » is the wavelength of the sound being produced (e.g. the horn). ”» is the change in wavelength due to the object s velocity, which is also known as the Doppler shift. In the car example, this corresponds to the change in pitch of the car s horn. Lastly, s corresponds to the speed of the wave, which for the car example is the speed of sound: 343 m/s or 767 mph. By using this equation, it is possible to predict ”» given v, or conversely to find v given ”». 12 Like sound, light is a wave, and thus Doppler shifts can occur in light as well. With light, the change in frequency from the Doppler shift corresponds to a change in the observed color of the light. So in principle, one could drive towards a red light so fast that the light would appear green. But we never notice Doppler shifts like these in our everyday lives, since light travels so quickly that no car, airplane, or even rocket ship can travel at a speed that is even remotely comparable to the speed of light. Thus the change in wavelength, or color, is not noticeable to our eyes. However, by carefully examining the spectrum of starlight using a spectrometer, it is possible to measure a Doppler shift in light. In particular, one can measure the Doppler shift in starlight caused by the motion of a star. Once a star’s Doppler shift has been measured, its corresponding velocity can be calculated using  REF _Ref139281264 \h Equation 1, with s being set equal to the speed of light. There is a small caveat to the detection of stellar velocities. From our vantage point on the Earth, we can only see light from a distant star that travels along a line-of-sight path. As such, using Doppler analysis, we can only determine the component of a star’s velocity that lies along this path. This is known as the star’s radial velocity. By measuring the Doppler shift found in a star’s light, astronomers are able to determine the radial velocity of star as it wobbles towards, or away from the Earth. This must be taken into account when fitting a curve to the data points (as in  REF _Ref138763294 \* MERGEFORMAT Figure 1.2), and interpreting the result. 13 Section 1.6: Doppler Analysis – Practice If a star is moving towards the Earth, then the wavelengths of light will be decreased in frequency, or blue-shifted. If the star is moving away from the Earth, the wavelengths of light will be increased, or red-shifted. In either case, the wavelengths of light passing through the spectrometer will not fall onto the same pixels as they would if the star were stationary. Stated another way, the orders in Figure 1.3 would be shifted to the left, or right, depending on the amount Doppler shift in the starlight. For a star that is wobbling due to an orbiting planet, the resulting Doppler shifts are so small that they cannot be detected by examining the spectrum with the naked eye. Instead, computers are used to analyze the spectrum and determine the amount of Doppler shift. A computer accomplishes this by comparing the spectrum from an observation to a baseline spectrum taken at a previous date. It examines the entire CCD, comparing the two spectrums, and looking for differences between the two. By examining subtle differences, the computer finally calculates a single Doppler shift for the entire observation. All Doppler shifts calculated in this way produce velocities that are relative to the baseline observation, for which the star is assumed to be stationary. 14 Section 1.7: PSF: The Spreading of Light Devices like telescopes and spectrometers are designed to record images as truthfully as possible. Unfortunately, it is generally not possible to capture an image perfectly. This is because telescope and spectrometers are made up of optical elements such as lenses and diffraction gratings, which affect the light that passes through them. This is best explained with an example. Imagine a tiny dot of light, such as that produced by a laser pointer. Attempts to image this dot through a lens will never result in a perfect image of the dot; the image will always be blurred to some extent. In this case the cause is small imperfections in the lens, which cause the light to spread out of focus. Even a source of light which is truly a point, such as a star, will be captured by a camera as having some measurable size on the film. The effect that an optical instrument has on the light which passes through it can be described mathematically by a Point Spread Function (PSF). Specifically, a PSF describes how an optical system spreads the intensity of a point source of light. Figure 1.6 shows a sample PSF. 15  Figure 1.6 – A Sample Point Spread Function (PSF). The dashed blue line in the center represents the intensity of an actual point of light. The black curve is the PSF, which shows how this intensity is spread out over several neighboring pixels on a CCD. An optical system like an Echelle Spectrometer, which does not utilize lenses, also cannot form perfect images. In this case, the effects of diffraction are responsible. While the diffraction of light is essential to the operation of the spectrometer, it also acts to degrade the information in the light that passes through it, thus degrading the information in the spectra formed. This effect can also be described with a PSF. 16 To provide an idea of how a PSF affects the spectrum of a star, I used a computer program to apply the PSF description shown in Figure 1.6 to the spectrum in  REF _Ref140840720 Figure 1.5. This process is known as convolution. The result is shown in  REF _Ref140840758 Figure 1.7.  Figure 1.7 – Application of the PSF shown in Figure 1.6 to the spectrum shown in  REF _Ref140840720 Figure 1. 1.5. The original spectrum is shown in black, while the spectrum modified by the PSF is shown in red. Note that the spectrum in red is a smoother version of the original spectrum, and thus contains less information. 17 The spectrum in red, which shows the result of convolving the original spectrum with the PSF, is a smoothed out version of the original. This effect is bad for astronomers, since smoothing out the original spectrum removes valuable information from the starlight In total, the PSF that affects an observation is a combination of the PSF from optics in the telescope, and the PSF of the spectrometer. Since materials in the telescope and spectrometer can be affected by environmental factors such temperature and humidity, the PSF description affecting observations can vary from night to night, or even minute to minute. Furthermore, the PSF description is not always the same spatially, and can affect the spectrum differently over different parts of the CCD. If a description of the PSF is known, a spectrum like the one in red can be deconvolved mathematically to regain a spectrum like the original. In general though, the PSF description is not known, and must be experimentally determined. The technique of accurately determining a PSF description, and using it to convolve and deconvolve a spectrum, is of particular importance in Doppler analysis. Chapter 2 – Doppler Precision of 3 m s-1 The previous chapter introduced the concept of detecting a planet by looking for periodic variations in the radial velocity of its host star. The amplitude of these velocity variations depends on the mass of the star, the mass of the orbiting planet, and the orbital characteristics of the planet such as the eccentricity and semi-major axis. Table 2.1 provides examples of the radial velocity amplitude induced on the Sun as a result of planets in our solar system. PlanetRadial Velocity Amplitude of SunJupiter12.4 m/sSaturn2.7 m/sEarth10 cm/s Table 2.1 – Radial velocity amplitudes induced on the Sun from three planets in our solar system. The technique of radial velocity traditionally yielded errors on the order of 1 km/s (Butler et al. 1996), and has slowly improved over time. By 1996, much work had been done to increase precision. According to Paul Butler and collogues, in the 1996 study Attaining Doppler Precision of 3 m s-1 , “current spectroscopic techniques yield Doppler-shift errors of 10 to 50 m s-1. This precision though was still “barely adequate to detect reflex radial-velocities caused by Jupiter-like and lower-mass planets” (Butler 1996). 18 19 In this study, Butler et al. presented a technique to attain Doppler analysis precision of 3 m s-1. While still not enough precision to detect low mass planets like Earth, this technique was a dramatic improvement over existing techniques. The paper describes the technique in some detail, providing mathematical reasoning for the attained velocity precision of 3 m s-1. In addition, the authors describe a computer program that implements the technique, and show the results of using this program to give evidence of their achieved precision. This computer program is still use today by our planet search team to perform Doppler analysis. In this chapter, I will discuss two main techniques used by the authors in their computer program to obtain this precision, as they pertain directly to my work on the Doppler analysis computer program. The first technique described in this study is the use of an Iodine cell. This is a glass cell containing molecular iodine gas. According to Butler et al., “The cell is mounted directly in front of the spectrograph slit in the converging f/36 telescope beam. Starlight passes through the cell just prior to entering the spectrographic slit. The iodine cell thus acts as a transmission filter, imposing thousands of extremely sharp lines on the starlight” (Butler 1996). These lines, which remain the same regardless of a star’s Doppler shift, act as a reference spectrum which aids in precision when performing Doppler analysis. In conjunction with a template observation of the star, this reference serves as a steady backdrop against which the stellar Doppler shift is measured, resulting in increased precision. In addition, the additional spectral information provided by the iodine helps the 20 Doppler analysis program determine the correct wavelength solution during Doppler analysis (see Chapter 3). The second technique important to mention here is that of chunks. Instead of trying to determine the Doppler shift of an entire spectral order at once, the “observation is broken up into several hundred 40-pixel (2 Å) chunks” (Butler 1996). The computer program determines the Doppler shift for each chunk independently, and then statistically averages the results of all chunks to yield a final radial velocity for the observation (see Chapter 3). There are two main reasons for creating chunks. First, the PSF is a function of position on the CCD. While it does vary slowly, this variation must be taken into account during the Doppler analysis process. By breaking up the observation into small chunks, the PSF can be assumed constant over each chunk during the Doppler analysis process. This simplifies dealing with a varying PSF. Second, it sometimes happens that a cosmic ray strikes part of the CCD during an observation. This renders the information in the affected pixels useless. By breaking an observation into chunks, only a few of the several hundred chunks will be affected by this, and as such a radial velocity for the observation can still be determined with little loss in accuracy or precision. Chapter 3 – The Doppler Analysis Computer Program The Doppler analysis programs are written in the IDL programming language, and executed on UNIX-like platforms (typically Sun OS or Apple OS X). An outline of the program is given in this chapter to familiarize the reader. The Doppler analysis program requires many pieces of input information to operate, which will be described in this section as the program is detailed. The main logic of the program is straightforward to describe. The code determines the radial velocity of a star by repeatedly comparing a model spectrum to an observed spectrum. Among the free parameters used in the model spectrum, the Doppler shift is varied until a good fit is achieved with the observed spectrum. Once a good fit is achieved, the radial velocity is calculated from this Doppler shift, and recorded. Section 3.1: Chunks and Observations The most fundamental input needed by the Doppler analysis code is the stellar observation to be analyzed. For our work at SFSU, this comes from the digital output of the CCD at Lick Observatory. The area of the CCD actually utilized is 1851 pixels wide and 1700 pixels high, containing about 90 spectral orders. The first step of converting this data to a useable format is the raw reduction. Here, each spectral order, which has a height of several pixel rows, is traced and mashed into a single value of intensity versus 21 22 pixel. The final result is an array of 90 rows by 1851 columns, where the array entry records the photon count for a particular pixel in an order. This reduced spectrum is passed to the Doppler code where 16 orders (absolute orders 39-53 on the Hamilton spectrometer at Lick) are analyzed to measure radial velocity. These orders span approximately 5090 Å – 5787 Å. They were chosen because most of the stellar radial-velocity information for late-type stars resides between 4000 Å - 6000 Å (Merline 1985). In addition, these orders are processed because iodine produces over one thousand lines in the 5000 to 6000 Å region (Butler 1987). In each order of a Lick observation, 44 contiguous groups of 40 pixels called chunks are established. Each 40 pixel chunk spans roughly 2 Å, and for a typical chunk contains 1-2 spectral lines. The reasons for using chunks are described in Chapter 2. With 16 orders, there are a total of 704 chunks. Doppler analysis is performed for each of these 704 chunks independently, resulting in 704 radial velocity measurements that must be statistically averaged. The individual results from all 704 chunks are stored in a velocity data file, typically referred to as a vd file, or just a vd. With this knowledge we can now present an illustration of how the Doppler analysis program performs its Doppler analysis for each chunk, and can introduce all inputs needed by the program. Please refer to the accompanying diagram in Figure 3.1 during this description. 23  SHAPE \* MERGEFORMAT Figure 3.1 – This diagram illustrates the flow of the Doppler analysis program as it determines the radial velocity for a single chunk. 24 Section 3.2: Wavelength Solution The Doppler analysis code generates a model spectrum for each of the 704 chunks. In order to get started, the code needs the range of wavelengths spanned by each 40 pixel chunk. Each pixel in a chunk corresponds to one resolution element of a stellar spectrum, and has a wavelength associated with it. Exactly which pixel a particular wavelength will fall onto depends on the characteristics of the spectrometer, and on the size, placement, and resolution of the CCD. To characterize this, an initial wavelength solution is needed. The first guess for this solution is drawn from a previous vd file, which is also provided to the Doppler analysis code as input. The wavelength solution for each order has non-negligible quadratic, cubic, and even fourth-order terms. Over the short, 2 Å range of a chunk though, the Doppler analysis code models the wavelength solution as linear. The linear wavelength solution for each chunk is stored with just two numbers: the wavelength zero-point and the dispersion. The zero-point is just the wavelength for the first pixel. The dispersion measures the linear change in the wavelength scale as a function of pixel. Thus to calculate the wavelength at pixel n (where n = 0 to 39), one would use the equation: (n = (0 + (n * d(/dpix) where (n is the wavelength at pixel n, (0 is the wavelength zero-point, and d(/dpix is the wavelength dispersion. 25 Section 3.3: IPCF and vdiods As mentioned, the wavelength solution is provided to the Doppler analysis program by referring it to a previous vd file (the result from a previous run of the program). The Doppler analysis code then extracts the wavelength mapping from this file. At first the method of using an older vd might sound strange. Why not simply record the wavelength mapping once, in some file, and then have the Doppler analysis program use it forever more for the mapping? This would be fine if wavelength solutions never change; but they do. The spectrometer and CCD are influenced by environmental factors such as temperature, and by mechanical settling of the instrument. Therefore, the wavelength solution may slowly drift with time. This is compensated for in two ways. First, the Doppler analysis code is provided with a list of vd files to extract the wavelength solution from, from which it picks the one closest in date to the observation being analyzed. In this way the Doppler analysis code is more likely to start off with a correct wavelength mapping. This list is held in an IPCF (Instrument Profile Configuration File) data structure. Each entry in the IPCF contains information about an observation, including the date of the observation, and the observation’s associated vd file. Second, the Doppler analysis code is allowed to vary the wavelength solution of each chunk as part of its Marquardt fitting routine (more on this later on). Once a best fit to the observed spectrum is achieved, the new wavelength solution for each chunk is recorded 26 in the output vd file. This newly created output vd file is then available for the Doppler analysis code to use for future wavelength solutions. Often the vd files referenced in the IPCF data structure are those created when the Doppler analysis code processed iodine-only observations. These are observations of the iodine cell’s spectral lines, illuminated either by a quartz lamp, or a rapidly rotating B star. Rapidly rotating stars have their spectral features significantly broadened, and as such act like giant lamps which simply illuminate the iodine cell’s spectral features. The advantage of using a B star over a quartz lamp is that the light which illuminates the iodine cell passes through the same telescope optics as the light from any observed star would. Hence, the spectrum of the iodine cell illuminated by a B star will resemble more closely the iodine spectrum seen during other stellar observations. The output vd file generated when processing an iodine only observations is referred to as a vdiod file, and is the preferred type of input for the Doppler analysis program. Section 3.4: The Stellar Template (DSST) With a wavelength solution in hand, the Doppler analysis code can begin to “crank” a model spectrum. Each chunk spans two Angstroms, so an example range is 5787.08 – 5789.02 Angstroms. The code obtains the corresponding spectrum for this same wavelength range from the Deconvolved Stellar Template (DSST). The DSST is a high resolution, high signal-to-noise observation of the star taken at a previous date, at one of 27 the observatories listed in Chapter 1. Typically it will have been taken at Keck or Lick Observatory. This observation serves as a template against which all radial velocities are determined. For practical purposes, the radial velocity of the DSST is set to 0 m/s (i.e., the absolute radial velocity of the star is subtracted), and all other radial velocities are measured with respect to this observation. A template observation differs from an ordinary observation in two ways. First, a template observation is taken without the Iodine cell present, and thus only has spectral lines from the star. Second, the spectrum of the observation has had the instrumental point spread function (PSF) deconvolved, such that the result better represents the spectrum before its light was smeared out through interactions with the optics and CCD. To accomplish this, a full PSF description (see Section 3.5) must be provided to the program dsst.pro. In addition to needing the locations and widths of the component Gaussian curves, the program also requires a height for each Gaussian. The heights are given by passing dsst.pro a list of vdiod files, from which the PSF heights are extracted and averaged for each chunk. If the star has a varying radial velocity, then the spectrum of most any observation of the star will be Doppler shifted from that of the DSST. Hence when forming a model spectrum, the code cannot simply use the pixels of the DSST’s spectrum that corresponds directly to the pixels of the observation (in this example: 5787.08 to 5789.02 Angstroms). Instead, it must take the portion of the DSST’s spectrum which would have been Doppler shifted to this range. This is done using a guess for the Doppler shift, which is another 28 input provided to the program. For example, a radial velocity of 29.8 km s-1 will correspond to a Doppler shift of about: (29,800 m s-1 / c) * 5787.08 Å = 0.580 Å This is illustrated in the Figure 3.2.  SHAPE \* MERGEFORMAT  Figure 3.2 – An illustration of the shift in a stellar spectrum due to Doppler shift. During the modeling process, the Doppler analysis program must select from the stellar template (DSST) the portion the spectrum appropriate to the Doppler shift being modeled. With about 2 Angstroms per chunk (40 pixels), a shift of 0.58 Å corresponds to roughly: (40 pixels / 2 Angstroms) * 0.580 Angstroms = 11.6 pixels 29.8 km s-1 may seem like an odd velocity to use for this example, since we are searching for stars with radial velocities on the order of 10 m s-1. However, a Doppler shift of roughly this size can easily be present due to the Earth’s motion (around the Sun) relative to the star being observed. Hence, this it is this barycentric velocity that must be accounted for when creating a stellar model. In addition, it must be subtracted from the radial velocity calculated in order to yield the radial velocity of the star due its orbiting planet. The barycentric velocity is calculated and recorded at the time of observation, and 29 is passed into the code as an initial guess for the Doppler shift Z. Our measurements for the barycentric velocity have a precision of 0.01 m/s (McCarthy 1995). After obtaining an appropriate section of the DSST’s spectrum, the corresponding piece of spectrum is selected from the Iodine Atlas (5787.08 – 5789.02 Å in our example). The Iodine Atlas is a high-resolution spectrum of the Iodine cell used at the observatory. Typically this atlas is made using a spectrometer at the Kitt Peak National Observatory in Arizona. By multiplying the chunk of the spectrum from the DSST and the Iodine Atlas, a model is generated for the spectrum of the star as observed through the Iodine cell. Section 3.5: The PSF The model is not complete, however, until the combined template and iodine spectrum is convolved with a Point Spread Function (PSF). The PSF is a mathematical function that models the effect of the telescope’s optics on the spectrum. The Doppler analysis program models the PSF by adding together 11 Gaussian functions (1 central Gaussian, and 5 adjacent “little Gaussians” on each side of the central Gaussian). The location and width (sigma) of all 11 Gaussian functions is passed as input to the code for a total of 22 input parameters that describe the PSF and which cannot be varied. These 22 inputs are typically referred to as the “PSF description”. Each Gaussian also has an associated height, which is varied during the Doppler analysis to accurately model the PSF (the height of the central Gaussian is always fixed, however). Because only the height can be 30 varied, it is important that the widths and locations of the Gaussian curves passed as input be selected appropriately so that the model PSF can accurately describe the actual PSF. This is typically determined by trying several different PSF descriptions, to see which works best. For example, a PSF which is thought to be narrow would be modeled by placing the little Gaussians close to the central Gaussian, with small widths. In contrast, a very broad PSF would be modeled by wider little Gaussians, spread further apart. The 11 heights of the Gaussians, in combination with the PSF description, result in 33 parameters that fully specify the shape of the PSF. During its operation, these parameters are passed to a function that constructs a PSF. By convolving the spectrum of the DSST + Iodine Atlas with this PSF, a model spectrum is created that closely matches the observed spectrum. Section 3.6: (r2 Fits In order to determine the goodness of fit between the model and observation, the (r2 (reduced chi-squared) statistic is calculated for each trial. As the model approaches the observed spectrum, the (r2 statistic decreases. In our application, (r2 has the following form:  31 Here observation and model are arrays that hold the observed spectrum and model spectrum (respectively) over the 40 pixels of the chunk. The array weight holds a weight for each pixel, which is based on the photon count for that pixel:  where countn is the photon count for the nth pixel of the chunk, and the epsilon term is the assumed error in the flat field, used to correct electronic variations in the CCD (e.g. from pixel-to-pixel variations in quantum efficiency). The epsilon is a very small value, so it is reasonable to think of the weight as just the inverse of the photon count. Pixels that receive a high photon count receive a low weight, thus reducing their influence on the calculated (r2. The motivation behind this is commonly encountered when constructing models of observations. Even a “perfect” model (one made with correct parameters) will not always match observation due to imperfections in the model inputs (such as the iodine atlas, or stellar template). In general, a smaller (r2 fit value will result when matching the stellar model to a lower signal-to-noise observation (lower photon count). The cause is the noise itself, which helps the imperfect observation and imperfect model resemble each other over the 40 pixels. In contrast, a higher signal-to-noise ratio in the observation (higher photon count), will result in a larger (r2 fit. As a conceptual example of this, consider comparing two photographs of the same subject. Two out-of-focus photos will tend to look more alike than will one blurry photo and one in-focus photo. 32 Finally, the term degf refers to the number of degrees of freedom. Dividing by this term distinguishes the reduced chi-squared fit from a normal chi-squared fit. The number of degrees of freedom is calculated as: degf = (number of pixels) – (number of free parameters) – 2 There are typically 40 pixels per chunk. With 13 free parameters (10 for heights of the PSF Gaussians, 2 for the wavelength solution, and 1 for Doppler shift), this makes degf = 40 – 13 – 2 = 25. A (r2 fit of 0 indicates a perfect match, though in practice this is not possible to achieve due to measurement uncertainties. Larger (r2 values indicating more deviation between the two data sets. The code attempts to minimize the difference between model and observation with the metric of lowest (r2 fit. When the best fit has been found, all 13 model parameters are recorded and the next chunk is analyzed. Section 3.7: The Marquardt Minimization Algorithm To find the lowest value of (r2, and thus the best fit between model and observation, a Marquardt fitting algorithm is employed. To explain the principle behind this algorithm, consider the simple example of a function f(x) shown in the Figure 3.3. 33  Figure 3.3 – An example of a function f(x) with a local minimum (at x2), and a global minimum (at x1). To find the value of x for which f(x) is a minimum, a Marquart algorithm searches for points at which df/dx = 0. To do this, a double-sided gradient search is carried out with x as a free parameter. In this process, stepwise changes are made for the free parameter around the initial best guess. So long as f(x) continues to decrease, the code continues to increment the free parameter in that same direction. If f(x) begins to increase, the code recovers the previous step for the free parameter. While this process can successfully find a location with df/fx = 0, it also allows the algorithm to get caught in a local minimum, such as the one around the value of the free parameter x2. If the stepsize for the free parameter is too small, or the initial best guess too far away, then global minimum (at x1) could be missed. The Marquardt algorithm is not restricted to only one-dimensional functions. It can, in fact, analyze functions of any number of dimensions. The (r2 fit used in the Doppler 34 analysis code is a good example of such a function. It is a function of thirteen different inputs: (r2 = (r2 ((0, d(/dpix, Z, gh1, gh2, …, gh10) Where (0 is the chunk’s wavelength zero-point, d(/dpix is the chunk’s linear wavelength dispersion and gh1 … gh10 are the heights of the little Gaussian functions which describe the PSF. The model spectrum that most closely fits the observed spectrum is found by identifying the thirteen values which minimize the (r2 fit value. To find this minimum, the Marquardt algorithm repeatedly varies all parameters over a thirteen-dimensional function space, continually calculating partial derivatives and exploring the function space. With each new set of parameters, the (r2 fit between the model and observed spectra is calculated. This process continues until a minimum of (r2 is found. Ideally, this minimum will be the global minimum if the model and spectrum are initially very close to each other (due to good initial guesses), It is often the case though that a very good fit between model and observation cannot be found. Sometimes this is the fault of the model, which fails to describe the observation accurately within the function space of its parameters. In the Doppler code, the model can fail because of an incorrect wavelength solution, an inadequate PSF description, or both. 35 Another cause of bad fits is due to the actual observation, and not the model. It is common for a cosmic ray to strike the CCD during an observation, which renders the information in the affected pixels useless. As a result, attempting to match model to observation in the region of the cosmic ray hit will prove fruitless. To accommodate failures in the modeling process, chunks with (r2 fits larger than an arbitrary cut-off are flagged in the output vd file so they can be ignored during future analysis. Typically, this cut-off (r2 value is 4.0. In addition, if the Marquardt algorithm fails to settle on a (r2 minimum within a certain time, the chunk is similarly flagged. Section 3.8: Frozen, Averaged PSF More than one iteration of Doppler analysis is carried out to improve the precision of the model. After the first pass, an output vd file is produced that records information on 704 chunks, including: Pixel and order location of the chunk (Doppler shifted relative to DSST). Pixel and order location of the DSST’s chunk. 10 Gaussian heights (describing the PSF) that gave the best fit. (0 and d(/dpix that gave best fit. Doppler shift that gave best fit. Value of best (r2fit. Value of radial velocity (based on Doppler shift). 36 One of the important parameters in the first pass of the Doppler code (above) is the PSF model. Physically, the PSF changes slowly across the CCD, as a function of the optics. To enforce the physically reasonable assumption of a slowly varying PSF, the PSF model is averaged over a few neighboring chunks. This is achieved with the following algorithm: Run a nominal pass of the Doppler code, passing an IPCF file from which to pick an input vd. Run a second iteration, passing the output vd from step 1 as its input vd (containing improved initial wavelength guesses), and using a frozen, averaged PSF. In the Doppler code, the initial guesses for the wavelength solution are taken from the vd generated in the first call. Thus, the Doppler analysis code begins the second call with a wavelength solution that already gave the best possible (r2 fits. In addition, the Doppler shift for each chunk is taken from the input vd, rather than simply adopting the barycentric velocity, as in the first pass. In the second call, the 10 Gaussian heights that describe the PSF are not varied. Rather, an average of the PSFs from surrounding chunks is adopted and the second Doppler pass is made with this frozen, averaged PSF and with only 3 free parameters (Z, (0 and d(/dpix). It has been experimentally determined that this method results in more homogenous PSFs across the CCD, and lower (r2 fits. Since lower (r2 fit values mean better matches between observation and model, they typically indicate a more accurate value of the Doppler shift, and hence the radial velocity. 37 When running with an averaged PSF, two additional parameters are passed into the second call to the Doppler code: del_ord and del_pix. These parameters specify the size of the area surrounding a chunk over which to calculate the averaged PSF. Typical values for Lick data are del_ord = (3 orders, and del_pix = (60 pixels (2 chunks). Figure 3.4 shows the chunks on the CCD which would correspond to these values (shaded). The center chunk marked with an ‘N’ indicates the chunk for which the averaged PSF is being calculated:  Figure 3.4 – An illustration showing “neighboring” chunks in an observation. The shaded chunks are the neighboring chunks of the chunk marked ‘N’, for the parameters del_ord = 3, del_pix = 60. Neighboring chunks are used to calculate an averaged PSF. Note that when an averaged PSF is calculated, it is always the PSFs from the input vd which are averaged, and never the new, averaged PSFs from the second iteration. 38 Since the PSF can be a function of wavelength, it might seem incorrect to include in the average PSFs from different orders (where the wavelength will be substantially different). However, Dr. Marcy has informed me that the PSF is not strongly dependent on wavelength, and its variation is largely just a function of 2-D position on the CCD. Hence chunks in different orders will have similar PSF descriptions so long as they are close to one another. Section 3.9: Velocity Analysis Code: A vd file contains velocity information about a single observation, but this by itself is seldom useful. Instead, a set of vd files are analyzed together to look for changes in a star’s radial velocity, and to assess uncertainties in each chunk. The velocity analysis code (a.k.a. “vank”) will: Determine a single radial velocity for each vd. Determine the corresponding error, in m/s, of this velocity. Determine the observation-to-observation RMS scatter of all velocities. Plot the results, and save the results in a vst file (velocity structure file). A sample output plot from this analysis is shown in Figure 3.5 for some observations of the star ( Ceti. 39  SHAPE \* MERGEFORMAT  Figure 3.5 – Sample output of the vank program. This program takes as input multiple vd files and calculates a velocity and internal error for each observation. It also calculates the standard deviation of these velocities (rBin). For these 25 observations, the median internal error was found to be ( 4.72 m/s, (rInt), and the observation-to-observation RMS scatter was ( 5.50 m/s (rBin). The observation-to-observation RMS scatter is simply the standard deviation of all the calculated velocities. In the current example, this would be the standard deviation of 25 velocities, which correspond to 25 individual observations. To calculate these velocities and their corresponding error, a weighting scheme is used to assign weights to each chunk. The weighting works as follows: each vd file is loaded into 40 memory as a two dimensional array. The first index of the array identifies a particular vd (or observation), and the second index identifies one of the chunks in that vd. For example, the nth chunk of the oth observation/vd would be accessed as: chunk = vdarray[o,n] Once the vd array has been created, the procedure vdclean identifies “bad” chunks. These are chunks which: Have DSST photon limited error above the 99th percentile. Belong to a chunk-set having a median (r2 fit above the 99th percentile. Section 3.10: Chunk-Sets and Weights: The second requirement calls for a new definition. A chunk-set is a collection of chunks, one from each vd, that are all located at the same position on the CCD (same order and starting pixel). A good way to visualize this is to imagine vds stacked on top of one another (like a pancake stack), and imagine a narrow beam of light (pencil beam) piercing through the stack perpendicularly. All of the chunks intersected by the beam would constitute a single chunk-set. Figure 3.6 shows an example of a chunk-set of 4 chunks created from 4 separate vd files. 41  Figure 3.6 – An illustration of a chunk-set made from four vd files. The concept of a chunk-set is also used when weighting chunks, which is the next step carried out in the vdclean procedure. Weighting chunks is important, as it allows poorly modeled chunks to be rejected when determining an average velocity for an observation. The idea behind weighting chunks is simple: assign higher weights to chunks that give consistent results from observations to observation. If chunks in a particular chunk-set always give consistent results, then we want to weight these higher than we would chunks in a chunk-set that give erratic results. In this section I’ll explain how this is accomplished. 42 If a planet were orbiting the star we were analyzing, we would not expect that chunks in a chunk-set would give consistent results, since the velocity of the star would change between observations. It is assumed in the code, though, that the star’s true radial velocity can be estimated by taking the median velocity of all chunks in an observation. By then subtracting out this median velocity from each observation, the actual stellar motion of the star will be subtracted out, and thus chunks in a chunk-set can be fairly compared to each other. As a first step towards this end, for each observation the median velocity of all chunks in the observation is calculated and stored. Then for each chunk in that observation, this median velocity is (temporarily) subtracted out. The resulting difference is recorded in an array called diff. The IDL code for this is: FOR ob = 0, num_obsvs - 1 do begin medvel = median( vdarr[ob,*].rvel ) diff[ob,*] = vdarr[ob,*].rvel – medvel ENDFOR Here the .rvel expression just means that we’re accessing the part of the chunk which stores the radial velocity calculated by the Doppler analysis program. Also, in IDL the * character is a wildcard that lists all array entries, rather than just a single entry. Hence the 43 argument to the function median is a list radial velocities, one from each chunk in the observation ob. Once the diff array has been calculated, the sigma array is calculated. This array has one entry for each chunk-set, and is populated with the following IDL code: FOR n = 0, num_chunks - 1 do begin igood = where( vdarr[*,n].rvel lt 99. ) sigma[n] = stdev( diff[igood,n] ) ENDFOR This piece of code loops through all chunk-sets (indexed by the variable n), and for each calculates the standard deviation of diff[igood, n], where igood is a list of chunks in the set that were not flagged as “bad”. Typically this list contains all chunks in the set, and so we can rewrite this statement more clearly as: sigma[n] = stdev( diff[*,n] ) Recall that diff records the differences between each chunk’s velocity, and the median velocity of the observation that the chunk belongs to. The array sigma, then, records the standard deviation of these differences for all chunks in each set. Both arrays can be visualized in Figure 3.7. 44  Figure 3.7 – Conceptual illustrations of the diff and sigma arrays, used by the vank code to determine chunk weights. For each chunk in each observation, the diff array stores the difference between the velocity for that chunk, and the median velocity of all chunks in its respective observation. The sigma array holds the standard deviation of these values for each chunk-set. The first table in Figure 3.7 above represents the diff array. Conceptually, each column is a chunk-set, and each row contains all chunks in an observation. Each horizontal blue line then represents a single entry in the diff array. The longer the line is, the larger the deviation between the chunk’s velocity and the median velocity of its observation (all deviations are represented in Figure 3.7 as absolute values). As you scan down a single 45 column then, you are seeing by how much the velocity deviates from the median along a chunk-set. The second table represents the sigma array, which stores the standard deviation of these differences for the chunk-set. If most chunks in a chunk-set have velocities that are consistently different from the median velocity for their respective observations, then this standard deviation will be small (look at chunk-sets 3 and 5 as examples). If, on the other hand, these differences across the chunk-set are erratic, then the standard deviation will be larger, as in chunk-set 2. With sigma calculated, the vdclean procedure finally calculates the aforementioned weights using the following IDL code: fudgefac = 2. FOR ob = 0, num_obsvs - 1 do begin const = median( abs( diff[ob,*] )/ sigma ) sigmaob = const * fudgefac * sigma for n=0, num_chunks - 1 do begin vdarr(ob,n).weight = 1./sigmaob(n)^2 endfor ENDFOR The heart of this code is on this line: vdarr(ob,n).weight = 1./sigmaob(n)^2 46 The array sigmaob is simply the sigma array (depicted in the Figure 3.7 above) multiplied by some constant factors (to be discussed later on). Hence, the weight assigned to a chunk n in an observation is proportional to 1/(sigma[n])^2. This means that chunks belonging to a chunk-set n with a small value of sigma[n] will be assigned a large weight, whereas chunks belonging to a chunk-set m with a large value of sigma[m] will be assigned a lower weight. Due to the fact that the weights are assigned as an inverse-square of sigma[n], very low values of sigma[n]result in very high weights. This type of weighting achieves exactly what was described earlier. Chunks that give consistent results from observation to observation (have similar diff entries along their chunk-set) are assigned high weights. Conversely, chunks that give erratic results from observation to observation (i.e. have erratic diff entries along their chunk-set) are assigned low weights. Chunks such as these would likely give erratic behavior due to bad modeling parameters (e.g., a suspect wavelength solution), and/or a difficult spectrum to model (e.g., one without many spectral lines, or perhaps one with too many lines). Before moving on, the constant factors mentioned earlier must be explained: const = median( abs( diff[ob,*] ) / sigma ) sigmaob = const * fudgefac * sigma 47 The value fudgefac was introduced at some point in the code’s history to increase the sigma values for each chunk in order to provide a better match to the observation-to-observation RMS errors. This is essentially an admission that at some level the assignment of observational uncertainty must be tuned up to be physically reasonable (e.g., match our expected uncertainty of 3 m/s). With current improvements in the Doppler code, this factor has recently been assigned to a value of 1.0 for the Keck velocities. However, during my thesis work it has always remained at the value 2.0. The work described here will help to assess if a different scaling factor is needed at Lick. Each entry in the array const measures where a chunk lies in the distribution of diff values for the chunk-set it belongs to. For example, such a ratio for a chunk in observation ob, and chunk-set n, can be calculated as: ratio = [ diff[ob, n] – mean( diff[*,n] ) ] / sigma[n] A ratio of less than 1 means that a particular chunk’s difference value is within one standard deviation of other chunks in the set, whereas a value greater than 1 indicates that the chunk’s difference value is outside one standard deviation. By taking the median ratio for all chunks in an observation, and multiplying this by the previously calculated sigma value, the code changes the weight of all chunks in an observation by a fixed factor. Observations with a small median ratio contain many chunks which agree well with the other chunks in their respective chunk-sets, and thus have their sigma reduced (and hence 48 weight increased by a factor of 1/const^2). Conversely, observations with a large median ratio contain many chunks which do not agree well with the other chunks in their respective chunk-sets, and thus have their weight reduced. Comments in the code describe this process as comparing “actual to typical discrepancy”. Note that this ratio will only compare apples to apples if the mean value of diff for a chunk-set is first subtracted from each chunk’s diff value (as shown in the example equation). This is because the standard deviation (sigma here) is always measured from the mean value of a distribution. Interestingly, this subtraction is not done in the code: const = median( abs( diff[ob,*] ) / sigma ) However, I’ve often seen that it’s typical for this value to be close to 0, and thus the calculation still works. In other words, I’ve often seen that for any given chunk-set, on average the difference of a chunk’s velocity from the median velocity of its respective observation is close to 0. This is not surprising to me. After all, it is implied in the code’s weighting scheme that we expect the true radial velocity of the star to be close to the median velocity of all chunks in an observation. Nonetheless, I suggest that Dr. Fischer or Dr. Marcy investigate this matter as well, to determine if this calculation should be changed in the code to match the example calculation. 49 Section 3.11: Mean Velocity, and Errors After the weighting of each chunk has been determined, a weighted mean velocity is calculated for each observation, accomplished with the following IDL code: velo = vdarr[ob,gd].rvel wt = vdarr[ob,gd].weight fit = vdarr[ob,gd].fit vind = where(fit lt 99) ; Ignore any “bad” chunks. velo=velo(vind) outcf[ob].mnvel = total(velo*wt)/total(wt) ; Weighted mean vel. The first several lines of code simply gather the velocity, weight, and (r2 fit for all “good” chunks in the observation. In this case “good” means 99% of all chunks, with the remaining 1% having the largest absolute difference between their radial velocity and the median velocity. The last line of the code simply calculates a weighted mean velocity for these chunks, storing the result in the outcf array. Finally, the internal errors are calculated. There are actually two types of errors to be calculated: Internal error for each observation. Observation to observation RMS scatter. 50 As mentioned earlier, the latter error (rBin on the plot) is just the standard deviation of the weighted mean velocities: rBIN = stdev( outcf[*].mnvel ) The internal error for each observation is calculated as: outcf[ob].errvel = 1./sqrt( total(wt) ) ;weighted error in mean Where the arrays wt and outcf are the same as those referenced previously. If we expand this expression, we obtain:  Where C is the const ratio, and f is the scaling factor for the internal errors. The final internal error reported on the plot (rInt) is just the median of the internal errors from each observation: rINT = median( outcf[*].errvel ) Chapter 4 – New Iodine Atlas In the Summer of 2005, I began investigating specific aspects of the Doppler analysis used at Lick Observatory. Dr. Fischer reported that the end orders in the iodine range showed gradually worsening model fits, particularly after the introduction of the newest CCD detector used on the Hamilton spectrometer in mid 2002. Because planet detection improves with higher precision, the general goal of this project was to: Reduce the internal errors of observations. Reduce the RMS scatter of velocities for consecutive observations (on timescales too short for the star’s radial velocity to have changed due to an orbiting planet). In order to achieve this goal, software visualization tools were developed and comparison tests were made for a few specific changes (testing one at a time) to the Doppler analysis. In particular: Comparison of velocity precision was made after introducing a higher resolution atlas, recently obtained from a scan of the Lick iodine cell with the Fourier Transform Spectrometer at Kitt Peak National Observatory. Comparisons were made using different PSF descriptions for the high resistivity CCD installed in 2002. Comparisons were made before and after enforcing wavelength continuity across the artificial chunk boundaries. 51 52 In this chapter I will discuss the first of these tests. The subsequent chapters of this thesis will address the remaining tests listed, as well as some other techniques that were tested. We tested the effect of changing the input iodine atlas to the Doppler code. A higher resolution scan of the iodine cell was obtained with the Kitt Peak Fourier Transform Spectrometer (FTS). Figure 4.1 shows an overlay of the higher resolution, new atlas (ftslick05), and the old atlas (ftseso50), over a range of about 1.5 Angstroms. The new atlas is plotted in red.  Figure 4.1 - An overlay of the higher resolution, new iodine atlas (ftslick05), and the old iodine atlas (ftseso50), over a range of about 1.5 Angstroms. The new atlas is plotted in red. 53 From the plot it can be seen that the new atlas has much higher resolution, and thus more spectral information, than the old atlas. The expectation was that integrating the new atlas would be an easy task to carry out, and that an atlas with more information should improve our current results. Section 4.1: A Baseline Test To compare the atlases a set of 16 observations of Ä Ceti (HD 10700) were analyzed with both atlases. Ä Ceti was chosen because it does not have any planets orbiting with a short period that are detectable within our current velocity precision limits (3 m/s). As such, radial velocity measurements taken on this star over several nights should generally yield similar results (within the limits of our precision). Both tests used the same PSF description (Gaussian locations and widths). The results of our test with the older atlas are shown first in Figure 4.2. These results served a baseline to compare all of or future results to. 54  Figure 4.2 – The vank plot from the baseline test of 16 Ä Ceti observations, analyzed using the old iodine atlas, with no modifications to the Doppler analysis program at SFSU. As indicated in Figure 4.2, the comparison parameters for this baseline test on HD 10700 were: Median Internal Error (rInt): 3.39 m/s Standard Deviation of Velocities (rBin): 3.35 m/s 55 Section 4.2: First Test with the New Atlas Having established this baseline, the new, higher resolution iodine atlas was tested. For the analysis of stars with the new FTS Iodine atlas, new vdiod files were generated for the test with the new atlas (referenced in a new IPCF file), along with a new stellar template (DSST) created using these vdiod files. The results are shown in Figure 4.3.  Figure 4.3 – The vank plot from a test of 16 Ä Ceti observations, analyzed using the new iodine atlas, with no other modifications to the Doppler analysis code at SFSU. 56 Table 4.1 shows the results as compared to our baseline test. Old Atlas (ftseso50)New Atlas (ftslick05)rInt3.39 m/s4.67 m/srBin3.35 m/s4.81 m/s Table 4.1  Comparison of rInt and rBin from analyzing 16 Ä Ceti observations using the old and new iodine atlases at SFSU. As can be seen in these results, the introduction of the new iodine atlas did not improve the velocity results for this test of 16 Ä Ceti observations. Both the internal error and observation-to-observation velocity scatter increased by more than 1 m/s. This result was counter-intuitive since the new, higher resolution iodine atlas has better wavelength information than the old atlas. One aspect of this project that is important to mention is that different parameters can be related. Modifying the FTS iodine atlas can requires a new input PSF description, so it was hard to assess this one result without making additional adjustments. To explain why a new PSF description might have been needed, consider the following. Let: IC = Pure spectrum of the iodine cell (no PSF). PK = Function that applies the PSF of the Kitt Peak Fourier Spectrometer. PL = Function that applies the PSF of Lick Observatory. 57 Given these definitions, then: PK(IC) = Our iodine atlas. PL(IC) = The observed iodine spectrum at Lick Observatory. When analyzing an iodine-only observation (to create a vdiod file), the Doppler analysis program determines a PSF P' for which (ideally):   PL(IC) = P'( PK(IC) ) In words, the PSF P' maps the iodine atlas to the observed spectrum of the iodine cell at Lick Observatory. Thus the PSF P' determined by the code is not equal to PL, the actual PSF of Lick Observatory, and in fact depends on the PSF of the Kitt Peak Fourier Transform Spectrometer. In the limit that PK is very narrow, such that is has no effect on the spectrum of the iodine cell, then P' will approach PL. The new iodine atlas was made using a much higher resolution scan than the older atlas, implying that the PSF applied at Kitt Peak was much narrower than before. Based on the previous argument, the Doppler code should then find a new P’. If this new P’ is sufficiently different from the old P’, a new input PSF description might be necessary for the Doppler analysis program to model it well. 58 Section 4.3: Fine Tuning the New Atlas To follow up on this result, a test program was written to overlay the spectrum of the older atlas on top of the new one. This demonstrated that the two spectra were always slightly offset from one another, and had slightly different dispersions. While a constant offset would not be a problem for our analysis, the observed difference in dispersion would be. The difference in dispersion was further verified by comparing the FTS scan of the Keck iodine cell against the old Lick atlas, which matched perfectly. An adjustment to the dispersion for the new atlas was made in the FTS-reading program. This program contains two atlas-specific values: the wave number for pixel 0 of the atlas, and the linear dispersion factor. By adjusting these two values slightly, the new iodine atlas matched up well with both the older atlas, and the atlas from Keck Observatory. The new and old values are shown in the Table 4.2. ftslick05Wave Number of Center of Pixel 0 (cm-1)Dispersion (Å / pixel)Original Values20983.3012700.01334083214Updated Values20983.2012700.01334078000 Table 4.2 – Original and updated wave number and dispersion values for the ftslick05 iodine atlas. The PSF description was modified to account for the intrinsically narrower iodine lines, and the new FTS atlas was used as input for another trial of the 16 Ä Ceti observations. The velocity precision did not improve much, and was not as good as the original analysis with the old atlas. The results are compared in Table 4.3. 59 Old AtlasNew AtlasNew Atlas (with fix)rInt3.39 m/s4.67 m/s4.46 m/srBin3.35 m/s4.81 m/s4.41 m/s Table 4.3  Comparison of rInt and rBin from analyzing 16 Ä Ceti observations using the old iodine atlas, and the new iodine atlas both with and without our fixes. Unwilling to give up on the new FTS atlas, we adopted the zero-point and dispersion change for the new atlas and continued testing other aspects of the Doppler code. I did not return to my analysis of the iodine atlas until 9 months later, as explained below. Section 4.4: Further Tests with the New Atlas Over the next 9 months I tried applying various techniques to improve our precision and to reduce both rInt and rBin to values lower than obtained with the old atlas. Using techniques discussed in Chapter 5, I was able to reduce one or the other of these values below their baseline value, but never both at once. In April of 2006, I attempted to port some of my experimental changes to an unmodified, older version of the Doppler analysis code at UC Berkeley. As a first step, the only change made to this code was the introduction of the new FTS iodine atlas. I used this code to analyze observations of the star 51 Pegasi (HD 217014), both with the unmodified code (old atlas), and with the slightly modified code (new atlas). The radial velocity plot from the test with the unmodified code is shown in the Figure 4.4. It is important to note that because observations of this star spanned multiple 60 dewars, observations were analyzed one dewar at a time. Figure 4.4 below shows only velocities calculated from observations taken with dewar 8:  Figure 4.4 – Radial velocity plot for observations of the star 51 Pegasi, analyzed with the unmodified code at UC Berkeley (old iodine atlas). Only observations taken on dewar 8 are shown. The results of my test with the new atlas are shown in Figure 4.5. As before, only the observations taken with dewar 8 are plotted. 61  Figure 4.5 – Radial velocity plot for observations of the star 51 Pegasi, analyzed with the code at UC Berkeley modified to use the new iodine atlas. Only observations taken on dewar 8 are shown. With the new atlas, the RMS scatter of observations from their Keplerian fit was larger than with the old atlas (8.83 m/s versus 6.60 m/s). Table 4.4 summarizes the median internal error and RMS scatter for the two tests. HD 217014 – Dewar 8Unmodified Code at UCBSame Code, Using New AtlasMedian Internal Error4.47 m/s5.56 m/sRMS Scatter6.60 m/s8.83 m/s Table 4.4 – Comparison of HD 217014 analysis done with the unmodified code at SFSU (old iodine atlas), versus this same code using the new atlas. Results shown are for observations taken on dewar 8. 62 For observations taken with dewar 8, the median internal error increased by more than 1 m/s. Results from other dewars also had larger RMS fits, as well as larger median internal errors, as shown in the Tables 4.5 and 4.6. HD 217014 – Dewar 13Unmodified CodeSame Code, Using New AtlasMedian Internal Error4.92 m/s5.71 m/sRMS Scatter5.47 m/s9.02 m/s Table 4.5 – Comparison of HD 217014 analysis done with the unmodified code at SFSU (old iodine atlas), versus this same code using the new atlas. Results shown are for observations taken on dewar 13. HD 217014 – Dewar 6Unmodified CodeSame Code, Using New AtlasMedian Internal Error5.51 m/s6.44 m/sRMS Scatter4.80 m/s5.66 m/s Table 4.6 – Comparison of HD 217014 analysis done with the unmodified code at SFSU (old iodine atlas), versus this same code using the new atlas. Results shown are for observations taken on dewar 6. Continuing to apply my other code modifications to the code at UC Berkeley, in addition to using the new atlas, gave similar results. As a final comparison between the two iodine atlases, I created a vdiod file using each. As discussed in Chapter 3, a vdiod file is the result of running the Doppler analysis code on an observation of a BStar or quartz lamp, such that the only spectral lines in the observation are those of the illuminated iodine cell. Doppler analysis on such an 63 observation fits a spectral model to these iodine lines using the iodine atlas. The resulting (r2 fits in the vdiod file indicate the goodness of fit between the atlas and the iodine cell spectra for each chunk. This provides a measure of how well the atlas can be used to model the iodine cell. Both vdiod files were generated using the same input parameters (PSF description, BStar observation, and IPCF entry). In Figure 4.6 I over plotted the (r2 fits of all 704 chunks from each vdiod file.  Figure 4.6 – Over plotted (r2 fits for the chunks of two vdiod files. One vdiod was created using the old iodine atlas (*), and the other with the new iodine atlas (+). Both vdiod files were created using the Doppler analysis code at SFSU. 64 The red + symbols indicate (r2 fits from the vdiod made with the new atlas, and the black * symbols indicate (r2 fits from the vdiod made with the old atlas. The plot shows that in chunks with larger chunk numbers (higher orders), the (r2 fits for the new atlas are typically worse. An identical test using a quartz lamp observation (instead of a BStar) gave very similar results. In hindsight, this vdiod analysis should have been performed when the new iodine atlas was first introduced, with efforts made then to see if PSF modifications could have improved the (r2 fits. In future work, I suggest applying this technique when other atlases are introduced. Due to time limitations, I decided to revert back to the old atlas for the remainder of my tests at SFSU. After doing so, I found that the rBin values in all tests I had been running improved. In particular, I was able to achieve both rInt and rBin values below the baseline values achieved in the analysis of the 16 Ä Ceti observations. Additionally, a better radial velocity fit was obtained for the observations of HD 217014 than I had analyzed with the unmodified code at UC Berkeley. These results will be presented in more detail in my Chapter 11. Note that unless otherwise stated, all tests discussed in this thesis before Chapter 11 were done while still using the new atlas. My final conclusion regarding the new FTS atlas is that it cannot be incorporated into our Doppler analysis code without much further investigation. It is even possible that the 65 atlas can never be successfully integrated. I was unable determine the source of the problem, but saw consistently that with the new atlas I could not obtain both rInt and rBin values as good as obtained with the old atlas. Chapter 5 – Polynomial Fitting, PSF Modifications, and IPCF Entries Section 5.1: Polynomial Fitting As mentioned in chapter 3, the Doppler analysis code extracts the initial wavelength solution from a vd file that is provided as input. This is accomplished via a call to the function sm_wav. For each order, this function extracts the wavelength zero points of all chunks and fits a weighted polynomial curve to them. The weighting favors chunks that obtained a low (r2 fit in the first pass of the Doppler analysis. The result is a polynomial equation that gives wavelength as a function of pixel. This smoothed, continuous wavelength solution is used when extracting ranges of spectra from the DSST and iodine atlas. In addition, the sm_wav function uses this wavelength solution to generate and return a vd structure. In each chunk of this vd, the two wavelength parameters reflect the smoothed wavelength solution at that location. This vd structure later provides the initial guesses for the Marquardt minimization algorithm. It is important to mention that the wavelength solution across an order is not linear. Approximating it as linear within a 2 Angstrom wide chunk is acceptable; however one cannot model it as linear across the entire order. Typically the solution has non-ignorable quadratic and cubic terms, and as such the default maximum order in the sm_wav function is pord=3. 66 67 We tried changing this value to pord=4, to see if increasing the resolution of our wavelength solution might improve our results with the new atlas. As a first test, we changed pord=4 and re-analyzed the 16 Ä Ceti observations. In addition, we also regenerated the DSST used with pord=4. Note that before doing this test we had begun our initial experiments on varying the PSF description (see next section), and as such the rInt and rBin values for the default pord=3 do not match those shown in Figure 4.3. The results of this test were: Value of pord Parameter3 (original default)4rInt4.45 m/s4.38 m/srBin5.49 m/s4.64 m/s Table 5.1  Results of changing pord from 3 to 4 when analyzing 16 Ä Ceti observations. The decrease in rBin made us confident at the time that increasing pord was a good idea. In addition, a plot of the PSF curve for a chunk with the median (r2 fit looked more symmetrical and Gaussian shaped with pord=4. Given these results we adopted this change in the code. As was to be seen, much better improvements could be obtained by further modification of the PSF description. Despite its promise of improvement, this modification of pord was unintentionally abandoned when we later switched to using an updated version of the code obtained from UC Berkeley. This fact was only noticed at a much later date. 68 Section 5.2: PSF Modifications – First Pass During Doppler analysis, portions of the iodine atlas are convolved with a PSF function to make them resemble the observed spectrum of the iodine cell. The newly introduced atlas, ftslick05, has much narrower spectral lines than the old atlas. At the time, we believed this meant the Doppler analysis code should use a narrower PSF function during the creation of spectral models to ensure that they closely matched observation. In retrospect this decision seems wrong, based on the discussion of the PSF in Section 4.2. With PK narrower, but PL fixed, P’ should have actually needed to be broader. Regardless, as will be seen later on, it was ultimately a narrower PSF description that gave the best results. I am unclear now on why a narrower PSF description gave the best results. We attempted to characterize an appropriate PSF description for the new atlas by processing the same set of observations with multiple PSF descriptions. To compare the results, we used the measures of internal error (rInt) and observation-to-observation velocity scatter (rBin) on the same set of 16 Ä Ceti observations used before. An accurate PSF description will produce good spectral models during the Doppler analysis, and so should reduce both of these measures towards their lowest values expected from intrinsic error (around 3 m/s). 69 In addition to these measures, I developed a tool to visualize the final PSF curve determined by the Doppler analysis code. As discussed in chapter 1, a PSF curve is typically Gaussian in shape and thus has a well defined full-width at half-max (FWHM). In addition, the PSF is a smoothly varying function of position, and so the curve and its FWHM should vary slowly across the CCD. My tool, called view_psf, visualizes the FWHM as a function of position on the CCD. The tool takes a vd file as input and extracts from it the PSF curve determined for each chunk. It then determines the FWHM of this PSF, and plots it on the vertical axis. Finally, the code uses built-in IDL graphics functionality to create a smooth topology over the plotted points. Figure 5.1 below shows a sample output of the view_psf program. The input to view_psf for this plot was the vd for observation rf87.32 from the baseline test of the 16 Ä Ceti observations analyzed with the old iodine atlas. Accompanying Figure 5.1 is a plot of the PSF function for the chunk in this vd with the median (r2 fit (Figure 5.2). 70  Figure 5.1  PSF FWHM contour plot for a vd file (observation rf87.32) of the baseline test of 16 Ä Ceti observations, analyzed with the old iodine atlas.  Figure 5.2  Plot of the averaged PSF for the chunk with the median (r2 fit. The vd file is for observation rf87.32, and is from the baseline test of 16 Ä Ceti observations analyzed with the old iodine atlas. 71 The FWHM does vary smoothly across much of the CCD, but shoots up near the end of each order, resembling a potato chip. Since the PSF should not change so dramatically over such a short distance, it was hypothesized that this feature was caused by bad spectral modeling near the end of each order, likely due to a bad wavelength solution. I ran tests with 22 different PSF descriptions, each with different values for the positions and widths for the Gaussians curves. After a first round of tests, I chose PSF description “A” (PSF A) as the best candidate based on the values it gave for internal error (rInt) and observation-to-observation velocity scatter, as well as its FWHM and PSF curve plots. The PSF descriptions for PSF A, and the old atlas, are: PSF A: pix = [0.00,-3.25,-2.60,-1.95,-1.30, -0.65, 0.65, 1.30, 1.95, 2.60, 3.25] sig = [1.00, 0.40, 0.40, 0.40, 0.40,  0.40, 0.40, 0.40, 0.40, 0.40, 0.40] Old Atlas PSF: pix = 1.3*[0.00,-3.25,-2.60,-1.95,-1.30, -0.65, 0.65, 1.30, 1.95, 2.60, 3.25] sig = [1.0, 0.40, 0.40, 0.40, 0.40, 0.40, 0.40, 0.40, 0.40, 0.40, 0.40] The pix array determines the locations of the Gaussian curves, and the sig array determines their widths. Compared to the PSF description used with the old atlas, PSF A has the Gaussian curves less widely spaced. 72 The results obtained with PSF A are shown in Figures 5.3 – 5.5.  Figure 5.3 – The vank plot from a test of 16 Ä Ceti observations, resulting from the application of PSF A. 73  Figure 5.4  Plot of the averaged PSF for the chunk with the median (r2 fit. The vd file is from a test of 16 Ä Ceti observations analyzed with PSF A.  Figure 5.5  PSF FWHM contour plot for a vd file from the test of 16 Ä Ceti observations, analyzed with PSF A. 74 I chose this PSF description to proceed testing with for 3 reasons. First, of all the PSF descriptions tried it gave one of the lowest pair of values for rInt and rBin, and a reasonable looking vank plot. The second two reasons that I chose PSF A were actually incorrect reasons at the time. I had provided the view_psf program with the wrong PSF description, resulting in very unusual looking plots. Of all these incorrect plots though, the plots for PSF A were the most reasonable looking, and so I picked this PSF. Despite this mistake, fortuitously I still made a good choice. The plots in Figures 5.4 and 5.5 were both generated correctly, and compared to the baseline results in Figures 5.1 and 5.2, both seem like reasonable results. The contour plot in Figure 5.5 also appears a bit flatter than its counterpart in Figure 5.1. It is important to mention that we do not know exactly what this FWHM plot is supposed to look like. However, given that the PSF is a slowly changing function of position, we do not expect it to undergo large changes over a short distance. Hence we hypothesized that a more correct PSF description will lead to a flatter FWHM plot. Section 5.3: IPCF Entries Prior testing on vdiod file creation had shown that the initial wavelength solution given to the Doppler analysis code (via the IPCF file) made a noticeable difference in the resulting PSF curves and (r2 fits. In addition, we knew that the PSF description and wavelength solution changed over time. Thus when proving an initial guess for the wavelength 75 solution, it made the most sense to use a vd file created from an observation temporally close to the observation being analyzed. The IPCF entries being used in our tests at this time referenced vd files generated long before the 16 Ä Ceti observations, taken on an older CCD. Given the sensitivity to the input wavelength solution seen, Dr. Fischer and I believed that our results could only benefit by referencing more recent vd files in the IPCF structure. I decided to update our IPCF entries to reference vd files created from observations that were temporally closer to the observations being analyzed. Rather than using a vd file resulting from just any observation, I created a vdiod file from a BStar observation taken on the same observatory run as the 16 Ä Ceti observations. My new IPCF file included only one entry, which referenced this vdiod file. The hypothesis was that this would provide a better initial wavelength solution, and thus improve results. We re-processed our previous test, done with PSF A, with this added improvement. The results of our test gave rInt = 4.31 m/s, and rBin = 3.78 m/s. These results were not much different than those obtained with the introduction of PSF A. However, both Dr. Fischer and I believed that using a more recent IPCF entry could only help matters, and so we adopted this change. This choice led to a future decision to always apply this logic when creating IPCF files, and will be discussed in Chapter 10. 76 Section 5.4: PSF Modifications – Second Pass With the changes to pord and our IPCF entry adopted, I decided to further refine our PSF description by trying 12 variations on PSF A. Of the results, I chose PSF description “RF L” (PSF RF-L) as the best. The description for PSF RF-L is: pix=1.1*[0.00,-3.25,-2.60,-1.95,-1.30, -0.65, 0.65, 1.30, 1.95, 2.60, 3.25] sig=-0.2+[0.9, 0.40, 0.40, 0.40, 0.40,  0.40, 0.40, 0.40, 0.40, 0.40, 0.40] Compared to PSF A, PSF RF-L spaces the Gaussians out a little more, but reduces their widths. The results obtained with PSF RF-L are shown in the Figures 5.6  5.8. The PSF plots are from the observation rf87.85. 77  Figure 5.6  The vank plot from a test of 16 Ä Ceti observations, resulting from the application of PSF RL.  Figure 5.7  Plot of the averaged PSF for the chunk with the median (r2 fit. The vd file is from a test of 16 Ä Ceti observations analyzed with PSF RF-L. 78  Figure 5.8  PSF FWHM contour plot for a vd file from the test of 16 Ä Ceti observations, analyzed with PSF RF-L. As with PSF A, this PSF description was chosen for three reasons. First, in conjunction with the other changes applied, it produced an rBin value (3.28 m/s) that was smaller than the rBIN value obtained with the old atlas (3.39 m/s). This was a major breakthrough. Also, the velocity versus time plot gave data points and error bars plausible for a star with radial velocity amplitude near 0 m/s, which we expect for Ä Ceti. Second, the PSF curve for the chunk with the median (r2 fit is nicely shaped and symmetric. 79 The last reason was the FWHM contour plot. This plot is far less flat than those shown in Figures 5.2 and 5.5. However, we hypothesized that this potato chip shape, though not physically correct, could be explained given other problems present in our analysis. To explain why, consider how the Doppler analysis code treats the PSF during modeling. During the Marquardt fitting, the Doppler analysis code can modify the applied PSF to account for nightly changes in the actual Lick PSF. Often though, in an effort to find the lowest possible (r2 fit, the PSF is instead contorted in an effort to account for noise, or even a poor wavelength solution. When studying the FWHM contour plot from the old atlas baseline test, we hypothesized that the rise in FWHM near the end of each order was caused by contortion of the PSF to account for high (r2 fits there. We further hypothesized that the cause of these high (r2 fits was a bad wavelength solution, which due to bad linear dispersion values caused the (r2 fits to grow larger over the length of each order. It was commonly seen that data from Lick Observatory, processed with the old atlas, would typically give higher (r2 fits near the ends of each order. A poor wavelength solution would explain these high (r2 fits and a contorted PSF, which could possibly give rise to the potato chip shaped FWHM plot seen. A bad PSF description, on the other hand, would have result in high (r2 fits for all chunks, regardless or position on the CCD, which was not seen. 80 At this point we felt confident that our PSF description was adequate and that the larger problem of the wavelength solution needed to be addressed before any more PSF refinements were made. While our PSF changes reduced rBIN dramatically compared to the original PSF, the reduction in rInt for all PSF descriptions tried was always less dramatic. The reduction of 0.43 m/s (4.67 m/s ’! 4.24 m/s) obtained with PSF RF-L was among the largest reductions seen. Chapter 6 – New Code and New Observations from UC Berkeley Section 6.1: EMU Observations Before proceeding with our tests on the wavelength solution, we decided to expand our data set. I made use of an IDL program written by John Johnson, a graduate student of Dr. Geoff Marcy at UC Berkeley. This program generates “end to end mockup” (EMU) observations based on the NSO solar atlas. In other words, it makes fake observations, using the spectrum of our sun and the iodine atlas as input. The caller of the program can generate these observations with any Doppler shift desired. The program adds noise, according to a given SNR parameter, and convolves the observation with a PSF curve that is also provided as input. Alternately, a simple Gaussian PSF curve can be requested. Typically, one has the program generate observations with Doppler shifts equal in value to the star’s random barycentric velocities. With a star’s barycentric velocity as input, the Doppler analysis code should then conclude that the star’s radial velocity is close to 0 m/s (as close as inherent errors will allow). I used this program to create 100 EMU observations, each with a radial velocity equal to its barycentric velocity, and convolved with a simple Gaussian PSF. All tests done with wavelength continuity were performed on this data set. Since the radial velocity of the 81 82 EMU star was known ahead of time to be 0 m/s, this data set gave me confidence in knowing whether the results of our modifications were sending us in the right direction. Section 6.2: Updated Code During my meeting with John Johnson at UC Berkeley to discuss his EMU program, I was informed that he was using a different version of the Doppler analysis code than I was. His version was more recent, and had certain improvements, probably added by Paul Butler. Because I was working with John Johnson’s EMU observations, which had been tested using this code, I decided to adopt this code for the remainder of my testing at SFSU. This newer code has 3 major differences from our code at SFSU: The program has an option to apply a smoothed, frozen wavelength dispersion. If this option is selected, the code will smooth the wavelength scale with a function called sm_disp (similar to sm_wav). It then freezes this dispersion so that it cannot be varied during the Marquardt fitting. Typically this option is selected on the second pass of the Doppler analysis program, along with the frozen, averaged PSF. The code still calls sm_wav to extract the wavelength solution from the input vd, but does not give the resulting smoothed wavelength solution to the Marquardt minimization algorithm as an initial guess. Instead, this guess comes either directly from the input vd file, or from the result of sm_disp. 83 During Doppler analysis, the code attempts to recover from “Bombs” (failure to fit model to observation) by re-analyzing a chunk using the modeling parameters of the previous chunk. At the time, I simply adopted this code and tested all my modifications against the results it gave initially (a new baseline). In addition, I stopped using the 16 Ä Ceti observations, and initially used only the 100 EMU observations for testing. As such, I never compared how this new code performed relative to our older code at SFSU. Recently, I checked this by re-running our baseline test of 16 Ä Ceti observations using the UCB code. The UCB code used in this test was exactly the same as the code I received from John Johnson, with no modifications. All parameters such as PSF description, iodine atlas (ftseso50), DSST, and IPCF file, were kept the same. The only difference in the test with the new code, besides the program itself, is that frozen dispersion was used in the second pass of the Doppler analysis. Figure 6.1 shows the results. 84  Figure 6.1 – The vank plot from a test of 16 Ä Ceti observations, analyzed with the newer, unmodified Doppler analysis code obtained from UC Berkeley (frozen dispersion invoked). Table 6.1 compares the results of this test to the original baseline test, done with our original code. Baseline Ä Ceti TestSFSU CodeUCB Code (with frozen dispersion)rInt3.39 m/s2.96 m/srBin3.35 m/s4.08 m/s Table 6.1  Results of analyzing 16 Ä Ceti observations with our original unmodified code at SFSU, versus the newer, unmodified Doppler analysis code obtain from UC Berkeley (frozen dispersion invoked). 85 The test shows that for this data set, the use of the new code with frozen dispersion invoked reduces internal error by 0.43 m/s. The rBin measurement, however, increased substantially by 0.73 m/s. Chapter 7 – Wavelength Continuity Section 7.1: Lack of Wavelength Continuity During Doppler analysis with both atlases, it was common to see (r2 fits that rose in value towards the ends of each order. It was hypothesized that this was due to a poor wavelength solution, which deviated from the true wavelength solution along each order due to an inaccurate linear dispersion values in the chunks. It was also believed that this was the cause for the potato chip shape of our FWHM contour plots. We hypothesized that correcting this wavelength solution would improve these (r2 fits, and as a result reduce our internal errors (rInt) and observation-to-observation velocity scatter (rBin). Because the spectrum projected onto the CCD is the output of an optical system, the wavelength solution is a continuously changing function across the CCD. However, one might not think this to be the case when examining a vd file created by our Doppler analysis. Table 7.1 below examines two neighboring chunks from an actual vd file. Chunk #4Chunk #5Start pixel: 198 Order: 38 (0: 5795.0933976 Angstroms d(/dpix: 0.0494791 Angstroms / pixelStart pixel: 238 Order: 38 (0: 5797.0695572 Angstroms D(/dpix: 0.0493706 Angstroms / pixel Table 7.1 – Information held in chunk #4 and chunk #5 of a vd file. The information shown specifies the location of the chunk on the CCD, and the wavelength solution for the chunk. 86 87 Upon close inspection there appears to be a problem with the wavelength solution. Using the values of (0 and d(/dpix from chunk 4 to calculate what (0 should be for chunk 5 gives: (0(5) = (0 (4) + (40 pixels * d( (4)/dpix) = 5795.0933976 A + (40 pixels * 0.0494791 A / pixel) = 5797.0725616 A This number is different from the actual value of (0 (5) by 0.0030044 Angstroms. Though this may seem like a small value, consider the velocity that a Doppler shift of this amount would correspond to at this wavelength: v / c = 0.003044 A / 5797 A v = 155.4 m/s Physically, the wavelength solution is a continuous function of position across the CCD. This means that the (0 predicted with the data in chunk 4 should match up perfectly with the actual (0 recorded in chunk 5. Because the Doppler analysis code is allowed to vary (0 and d(/dpix as part of its fitting routine though, there is no enforcement of this continuity in the wavelength solution. Factors such as noise, low photon count, and poor PSF descriptions can thus result in incorrect wavelength solutions that are chosen because they provide the lowest (r2 fit. As a result, discontinuities appear in the wavelength solution. 88 In one form of this discontinuity, the ending wavelength of a chunk is less than the starting wavelength of the next chunk. In the other form, the ending wavelength of a chunk overlaps the starting wavelength of the next chunk (as in the example). For both cases, the difference between the starting wavelength of a chunk and the ending wavelength of the previous chunk is referred to as a wavelength “gap”. Section 7.2: Visualizing the Gaps I wrote an IDL program to visualize these gaps across each order. The program takes a vd file as input, and from it determines the gaps between adjacent chunks. It then plots the gaps across each order. It additionally provides a histogram of gap sizes. A sample output is shown in Figure 7.1 and Figure 7.2. In this case, the input file was a vd file from my baseline test of 16 Ä Ceti observations using the old iodine atlas. Though not clearly visible, the y-axis reads  wvln gaps , is measured in Angstroms, and spans both positive and negative values. The x-axis ranges from 0 to 703, indicating which chunk pair the gap is being plotted for. 89  Figure 7.1 – Sample output of the view_gaps program. This output plots the wavelength gaps between neighboring chunks for each order. The gaps shown are of a vd file from my baseline test of 16 Ä Ceti observations using the old iodine atlas. 90  Figure 7.2  Sample output of the view_gaps program. This output is a histogram of wavelength gaps between neighboring chunks for all orders. The gaps shown are of a vd file from my baseline test of 16 Ä Ceti observations using the old iodine atlas. The program also provides statistical information about the gaps. Note that this information is calculated using the absolute value of all gaps. The output for this vd is shown in Table 7.2, with all values given in Angstroms. Largest gap size: 0.016898096Smallest gap size: 0.000000000Median gap size:0.001343250Standard deviation: 0.001761757 Table 7.2 – Wavelength gap information for a vd file from the baseline test of 16 Ä Ceti observations using the old iodine atlas. 91 Section 7.3: Freezing the Wavelength Solution The ability of the newer code to smooth and freeze dispersion was similar in concept to our goal of making a continuous wavelength solution. However, a simple application of the view_gaps program shows that freezing only the dispersion does not result in a continuous wavelength solution. Figure 7.3 shows the output of the view_gaps program when run on a vd file that had frozen dispersion invoked during the second Doppler analysis pass. This test was performed using the old iodine atlas. 92  Figure 7.3 – Output of the view_gaps program for a vd file generated with the new code from UC Berkeley. Frozen dispersion was invoked on the second pass of the Doppler analysis. 93 Table 7.3 shows the statistical analysis on the absolute gap sizes, calculated by the program, with all values given in Angstroms. Largest gap size: 0.022881981Smallest gap size: 0.000000000Median gap size:0.003153950Standard deviation: 0.003755130 Table 7.3 - Wavelength gap information for a vd file generated with the new code from UC Berkeley. Frozen dispersion was invoked on the second pass of the Doppler analysis. These results show that even with a smoothed, frozen dispersion, the ability of the Marquardt minimization algorithm to change the wavelength zero-point of each chunk can still result in a discontinuous wavelength solution. In fact, the mean absolute gap size is approximately 2.4 times larger than the value obtained in our baseline test of 16 Ä Ceti observations. Section 7.4: New Baseline Tests Before applying wavelength continuity, I analyzed 100 EMU observations with the new code to establish a baseline. Later on, I also performed this test using the 16 Ä Ceti observations, as well as an additional 9 Ä Ceti observations, for a total of 25. For each data set, frozen dispersion was invoked on the second pass of the Doppler analysis. The velocity plots from each of these tests are shown in Figures 7.4 and 7.5. 94  Figure 7.4  The vank plot from a test of 100 EMU observations, processed with the new code from UC Berkeley.  Figure 7.5  The vank plot from a test of 25 Ä Ceti observations, processed with the new code from UC Berkeley. 95 Section 7.5: Methods to Enforce Wavelength Continuity Dr. Fischer and I decided on and tested a method to enforce wavelength continuity. First, we would freeze not only the dispersion, but also the wavelength zero-points in each chunk. Second, we would only freeze these parameters after a call to the sm_wav function, which would guarantee a continuous wavelength solution. Lastly, instead of the usual two passes of the Doppler analysis program, we would make three passes: First pass: Run the Doppler analysis as usual. Second pass: Smooth the wavelength scale, and then freeze the wavelength dispersion and zero-points for all chunks. Run the Doppler analysis with the remaining free parameters (for PSF and Doppler shift). Unfreeze the wavelength paramet