Temperature dependence of the eShel spectrograph for precise radial velocity measurements

(1)

University of Amsterdam

Anton Pannekoek Institute, FNWI

Bachelor Thesis

Temperature dependence of the eShel

spectrograph for precise radial velocity

measurements

Author:

Justin van Lierop 10262180 Supervisors: Prof. Dr. H.F. Henrichs MSc. S.M. Straal Examiners: Prof. Dr. H.F. Henrichs Prof. Dr. L. Kaper 31-03-2014 - 12-07-2015

(2)

Abstract

At the Anton Pannekoek Observatory, the two most recent radial velocity measurements of the exoplanet hosting star υ And using the fibre-fed Shelyak eShel spectrograph and 50 cm RCOS telescope were unsuccessful, suffering from unexpectedly large errors. This research investigates whether this was caused by temperature changes affecting the spec-trograph. New radial velocity measurements were done of the exoplanet hosting star τ Boo, during which spectrograph temperatures were measured. Additionally, wavelength calibrations were done at different temperatures. By processing a single τ Boo spectrum using these different calibrations, the direct wavelength shift induced by the temperature changes could be determined by processing the spectra as radial velocity measurements. The resulting wavelength shift expressed in velocity is -0.29±0.43 km/s at a ∆T of 0.5◦C and -1.46±0.53 km/s at a ∆T of 1.2◦C. The radial velocity measurements of τ Boo were not accurate enough to measure the expected temperature effect. These large errors, and probably the large errors of the earlier measurements at APO appear to be caused by a less than optimal focusing of the collimator of the spectrograph where the light exits the fiber. This has likely affected all subsequent measurements of the eShel spectrograph since 2 November 2012, and should urgently be corrected.

(3)

English Summary

At the Anton Pannekoek Observatory (APO) (University of Amsterdam) four attempts have been made to determine the orbit of exoplanets using the fiber-fed Shelyak eShel spectrograph. The first attempt to be successful was of τ Boo by J.Camphuijsen et al.

(2012), while the subsequent two of υ And by Boots & Meinster (2012) and Barendse & van Lierop (2013) were not. For measurements to succeed, accuracy and stability of the spectrograph are crucial. Different measurements of the same object may be equally accurate, yet still yield different results. This can be caused by a lack of stability, on which temperature is the dominant influence. This is the main investigation of this research. Which wavelengths correspond to the given pixels of the CCD camera is de-termined with the wavelength calibration. This is an exposure of a Th-Ar lamp of which the spectral lines are well-known. Because the mechanical structure of the spectrograph expands and contracts under influence of temperature, the lightpath inside the spectro-graph changes as well. If such a change occurs between the wavelength calibration and science exposure, the wavelengths of the spectral lines would become uncertain. Both the shift in wavelength of the Th-Ar calibration and τ Boo observations as a function of temperature have been investigated. The eShel spectrograph was attached to the 50cm RCOS of the Anton Pannekoek Observatory. A thermometer placed on the outside of the box shielding the eShel from stray light was used to measure all temperatures. Four exposures of τ Boo with suitable S/N were obtained using this setup, with aver-age spectrograph temperatures of 24.25, 24.0 and 23.75◦C. These were processed using a Mathematica 10.1 cross-correlation program written for Doppler spectroscopy along with nine other spectra from similar observations done over the course of the month after. Unfortunately, the spectra appeared to be slightly out of focus, something that should have, but was not noticed for over two years. This resulted in errors too large to determine a significant temperature effect. This was most likely caused by the eShel fiber not being properly focused on the collimator. This appeared to have happened to the observations of Boots & Meinster (2012) and Barendse & van Lierop (2013) as well, as the spectra of the latter were similarly out of focus, which would explain the large errors they suffered. In addition, Th-Ar exposures at different temperatures were

(4)

iv

used as wavelength calibrations for a single τ Boo spectrum, resulting in three spectra calibrated at different temperatures, 21.0, 21.5 and 22.2◦C. The relative velocities of these spectra were calculated using the same Mathematica 10.1 program and resulted in a shift of -0.29±0.43 km/s at a ∆T of 0.5◦C and -1.46±0.53 km/s at a ∆T of 1.2◦C. This is similar to the results fromGarde(2015), who determined the pixelshift of Th-Ar spectra as a function of temperature. Translated to velocities, this resulted in 1.6 and 2.3 km/s for a∆T of 0.5 and 1◦C respectively.

(5)

Dutch Summary

De planeten van ons zonnestelsel zijn niet de enige planeten in het heelal. Sinds de jaren ’90 worden er ook planeten gevonden die in een baan rond andere sterren bewe-gen en deze worden exoplaneten bewe-genoemd. Deze zijn echter niet zomaar zien met een telescoop, want de ster overstraalt het licht van de exoplaneten. Hoewel een ster veel zwaarder is dan een planeet, oefent de planeet ook een beetje zwaartekracht uit op de ster. Hierdoor beweegt de ster ook, in een baan om het massamiddelpunt. Dit betekent dat als het baanvlak niet haaks op onze kijkrichting staat, de ster soms naar de aarde toe beweegt, en soms van de aarde weg. Deze beweging zorgt voor een Doppler-effect in het spectrum van de ster. Net als geluidsgolven van een auto worden samengedrukt of uitgerekt als deze langsrijdt en daardoor hoger of lager klinken, zorgt het Doppler-effect dat licht roder of blauwer wordt als de lichtbron weg van of naar de waarnemer beweegt, omdat de golflengte langer of korter wordt. Omdat een ster een spectrum uitzendt met karakteristieke absorptielijnen, is dit Doppler-effect te zien met een spectrograaf als donkere (absorptie) lijnen die over het spectrum heen bewegen. Omdat zo’n ster re-latief langzaam beweegt ('100 m/s), is de verschuiving van deze absorptielijnen heel erg klein, en is er een hele nauwkeurige spectrograaf nodig om dit waar te nemen. De eShel spectrograaf van het Anton Pannekoek Observatorium zou hier nauwkeurig genoeg voor moeten zijn, maar de laatste twee pogingen om de beweging van de ster υ Andromeda waar te nemen zijn niet gelukt. De metingen waren te onnauwkeurig, en het vermoe-den is opgekomen dat dit zou kunnen komen door fluctuaties in de temperatuur van de spectrograaf. Door Doppler-spectroscopie toe te passen op de ster τ Bo¨otis, en door temperatuursveranderingen van de golflengtecalibratie te bestuderen is het de bedoeling om de temperatuursafhankelijkheid van de spectrograaf te testen. Welke golflengte bij welke pixels horen wordt bepaald met een golflengtecalibratie. Het spectrum van een Th-Ar lamp heeft karakteristieke spectraallijnen met bepaalde golflengten, waardoor met geschikte software de betrekking tussen golflengte en pixels bepaald kan worden. Als de spectrograaf warmer wordt zet deze uit, en wordt het lichtpad in de spectro-graaf anders. Hierdoor komt licht met een bepaalde golflengte op andere pixels van de CCD terecht, waardoor niet meer te zien is of spectraallijnen zijn verschoven door

(6)

Contents vi

het Doppler-effect of door deze temperatuursverandering. Door golflengtecalibraties te doen bij verschillende temperaturen, en deze vervolgens te gebruiken om een τ Boo spec-trum te calibreren worden verschillende spectra met spectraallijnen die onderling alleen verschoven zijn door de temperatuursverschillen verkregen. Dit resulteerde in een ver-schuiving van -0.29±0.43 km/s voor een ∆T van 0.5◦C en -1.46±0.53 km/s voor een ∆T van 1.2◦C. Dit is dus een grote verschuiving, maar door de golflengtecalibratie voor en na elke waarneming te doen wordt dit effect sterk verminderd. De temperatuureffecten van de snelheidsmetingen van τ Boo vielen echter binnen de fout en waren dus niet goed waarneembaar. Dit komt omdat de spectra niet scherp genoeg waren, doordat de eShel niet goed stond afgesteld. Het lijkt er zeer sterk op dat dit ook het geval was bij de vorige twee onderzoeken met deze spectrograaf.

Figure 1: Dit is een voorbeeld van een eShel spectrum. Het totale spectrum is opgedeeld in verschillende banden. De kleine verticale donkere strepen in de witte

(7)

Introduction

While mankind has gazed up to the night sky for thousands of years, for a long time it was not known whether earth and the other planets in our solar system were the only planets in existence. The first planet orbiting a star other than our sun was discovered byWolszcan & Frail(1992), and the first planet orbiting a solar-type star was discovered byMayor & Queloz(1995). These are known as extrasolar planets, or exoplanets. Since then, many more exoplanets have been found as technological advances have been made. As of now, it is believed almost every star has at least one planet (Cassan et al. 2012).

1.1 Exoplanet Detection

Due to the faintness of the planet compared to its host star and the distances involved, for a long time it was impossible to directly observe an exoplanet from earth. This method of detection has been successful only quite recently, the first time being presented by

Chauvin et al. (2004). There are several options to directly observe exoplanets despite their low luminosity compared to their host star. With a coronagraph for example, it is possible to block most of the starlight. The luminosity of the star is also relatively less in the infrared region compared to exoplanet’s luminosity. Large planets orbiting far from their star are easier to observe directly. Levine & Soummer(2009).

The number of exoplanets observed directly may have increased since then, but many exoplanets are still being discovered by either the transit or radial velocity method, as shown in Figure 1.1. The transit method entails the monitoring of the brightness of a star. Should an exoplanet orbiting the star move in between the earth and the star it would block part of the star’s light from reaching earth, noticeable as a slight decrease in the brightness of the star for a relatively short time. This is known as a transit (Figure

(10)

Chapter 1. Introduction 2

Figure 1.1: The number of exoplanets found per detection method since the first discover, based on information of the Open Exoplanet Database. Database (2014) Those found with the transit method is shown in green, radial velocity method in blue and direct imaging in red. Yellow are those found by timing variations and orange for

microlensing. The high peak in 2014 is caused by the Kepler discoveries.

1.2). Because the planet orbits the star periodically, these transits appear with the same periodicity, which would confirm the existence of an exoplanet.

Figure 1.2: A schematic depiction of an exoplanet (black) in orbit around a star, as well as the characteristic light curve of an exoplanet transit. As the planet(black) moves in front of the star, the star’s brightness decreases. Picture from Hendersson

(2013)

When a planet orbits a star, both are actually orbiting the center of mass of the system due to their mutual gravitational attraction. This means that as seen from the earth, the star periodically moves towards and away from the earth if its orbital plane is not perpendicular to the line of sight. This movement induces the well-known Doppler

(11)

effect in the spectrum of the star, making the spectral lines appear blue- and redshifted. This shift in wavelength can be measured with a spectrograph, and from this the radial velocity of the star can be derived. If the radial velocity is found to vary periodically, it is a strong indication that an (or more) exoplanet orbits the star. Using Kepler’s laws, a lower limit of the mass ratio of the planet and star can be calculated. If the mass of the star is known, a lower limit of the planet’s mass can be determinedMayor & Queloz

(1995). The method described above is the radial velocity method, and it requires very accurate measurements and equipment as typical velocity amplitudes are smaller than 100m/s.

Figure 1.3: An exoplanet (brown) in orbit around a star (yellow). The black dot is the center of mass of the system, and the arrows indicate the motion of the star and planet around it. As the star moves away from the observer, its spectral lines become redshifted, as shown at the left spectrum. When the star moves away from the observer, the spectral lines become blueshifted, as shown at the right spectrum. The shift of the absorption lines is exaggerated, as often the amounts shift less than a pixel

on the detector.

As the limiting factor for measuring exoplanet orbit is the accuracy and stability of the spectrograph, it is very important to have clear insight which factors determine these properties. For instance, Buil (2012) claims with similar instrumentation a highest accuracy of 50m/s over several months, allowing the determination of 4 exoplanets with a high radial velocity amplitude. Attempts to replicate these results at the Anton Pannekoek Observatory (APO) have had mixed success. J.Camphuijsen et al. (2012) has succeeded in confirming the orbit of τ Bo¨otis (τ Boo) (Figure1.4), with typical errors of 100 m/s. Similar measurements of υ Andromeda (υ And) byBoots & Meinster(2012) and Barendse & van Lierop (2013) failed however (Figure 1.5), for unknown reasons, possibly related to the long-term changes of the instrument. This is the subject of this investigation.

(12)

Chapter 1. Introduction 4 0.0 0.2 0.4 0.6 0.8 1.0 -_1.0 -0.5 0.0 0.5 1.0 Phase Velocity Hkm s L

APO 2011

Figure 1.4: Successful measurement of τ Boo byJ.Camphuijsen et al.(2012). This is a phasediagram with radial velocities shown in km/s, showing errors as one standard

deviation, typically 100 m/s.

Figure 1.5: Unsuccessful measurement of υ And byBarendse & van Lierop (2013). This is a phasediagram with radial velocities shown in km/s, showing errors as one

(13)

This study investigates the temperature dependence of the Shelyak eShel spectrograph at the APO. Four observations of τ Boo, which hosts an exoplanet, from 2014 were used. A decrease in temperature from 24.5 to 23.5◦C was noticed during these observations. By calculating the radial velocities from these spectra, it can be determined whether a change in temperature influences the expected results.

The outline of this thesis is as follows: The second chapter gives a description of the radial velocity method in theory. In chapter 3, the specifications and use of the instrumentation used in this research is given, as well as a description of the acquisition of spectra. Chapter 4 then describes how the acquired spectra are processed and how a radial velocity is determined. It also elaborates on the crosscorrelation of spectra. The results of this research are given in chapter 5, after which these are discussed in chapter 6. This chapter also outlines the problems encountered during the observations. Chapter 7 concludes the research.

(14)

Chapter 2

Theoretical Background

The mutual gravitational attraction of a star and a planet causes both to orbit the center of gravity of the system. During this periodical movement, the star moves both towards and away from the earth if the orbital plane is not perpendicular to the line of sight. This induces a Doppler shift of a spectral line at laboratory wavelength λ0, which relates

to the radial velocity.

λ − λ0

λ0

= v

c (2.1)

Here the λ is the observed wavelength, c the speed of light and v the radial velocity. If the radial velocity as a function of time varies with a well defined period, this corresponds to the period of the orbit of the exoplanet and parent star around the center of gravity of the system. Accurate determination of the radial velocity is therefore important, as the size of the errors determines the limit of detectability of any periodicity. With the period P determined, the third law of Kepler can be used to find a, the semimajor axis of the orbit of the exoplanet.

a3 P2 =

G

4π2(M + Mp) (2.2)

With M the mass of the star and Mp the mass of the planet and G the gravitational

constant. The radial velocity as a function of time also provides a semi-amplitude K, or half the change in radial velocity from the maximum to the minimum over the course of an orbit. Together with the period this gives the so-called mass function, which gives an estimate for Mp: P K3 2πG = Mpsin3i (1 + q)2 (2.3) 6

(15)

Chapter 2. Theoretical Background 7

Here q is the mass-ratio of the star and planet and i is the inclination angle of the system. Because the inclination is often unknown, only a lower mass limit can be determined using this method.

(16)

Chapter 3

Instrumentation

All observations were done at the Anton Pannekoek Observatory (APO) of the University of Amsterdam, located at Science Park. The telescope used for observations is a 50 cm f /8.2 RCOS Ritchey-Chr´etien telescope on a 10micron GM4000 mount, with dome automation (ACE) which allows slaved dome operation. The telescope has a port selector (Van Slyke Sidewinder) which allows three light paths, usually assigned to an imager, an eyepiece holder and the guiding unit of the eShel fibre-fed spectrograph. Due to necessary repairs of the Sidewinder, for some observations, a flipmirror was used instead. A focal reducer f /6.3 is placed between the sidewinder and the guiding unit to match the focal ratio of the telescope with the acceptance cone of the fiber.

Figure 3.1: Photograph of the telescope, showing the guiding unit, port selector and eyepiece.

(17)

Chapter 3. Instrumentation 9

Figure 3.2: Schematic depiction of the components used in the observatory. The red line represents the 50µm fiber transporting light from the telescope to the spectrograph, the green line represents the 200µm fiber transporting light from the calibration unit

to the guiding unit.

3.1 eShel Spectrograph

Spectra were taken using the Shelyak eShel fibre-fed spectrograph, which has a spectral resolution R=_∆λλ =10000 over the wavelength range 4450-6750 ˚A. The spectral resolu-tion indicates the smallest difference in wavelength ∆λ distinguishable at wavelength λ (Rayleigh criterium). This spectrograph is suited for radial velocity measurements due to its broad range of wavelengths. This provides more spectral lines that can be used for the radial velocity measurements, resulting in a higher accuracy as compared to a higher resolution instrument like the LHIRES III with R=17000, which is higher, but with a wavelength range 10 times smaller. Light collected by the telescope is directed into the eShel Guiding unit, which is attached to the telescope at the port selector or flipmirror. Here the light falls upon a mirror that reflects the light into the guiding camera. A 50µm wide hole in the mirror allows light to pass through into a 50µm optical fiber that guides the light to the spectrograph.

The spectrograph is placed in a temperature controlled room (by the building), with the temperature roughly between 20 and 25◦C. The spectrograph is placed in a box to

(18)

Figure 3.3: Shelyak eShel spectrograph and calibration unit in the APO control room on the left, schematic cross-section of eShel on the right. (Cochard & Thizy 2014)

shield it from any stray light in the room. During measurements, the room is kept dark. The calibration unit is next to the spectrograph. It contains a thorium-argon lamp for a wavelength calibration and a tungsten lamp with led light for flatfield exposures. During these calibrations, the light of the calibration unit is transported to the telescope through a 200µm optical fiber which enters the eShel guiding unit and follows the same path as any light collected by the telescope. The light entering the spectrograph first passes through a diffraction grating. The dimensions of this grating are such that most of the light is reflected in the higher orders which strongly overlap. The light then passes through a so called cross-disperser, to separate the orders perpendicular to the dispersion direction. In this case, a prism is used. In order for constructive interference to occur for light passing the grating, the path difference between light from adjacent grooves must be an integer multiple, m, of the wavelength λ. The initial dispersion occurs according to the grating equation:

mλ = d(sin θi+ sin θr) (3.1)

d is the distance between the grooves and θi and θr are the angles of the incident and

outgoing light respectively, as measured from a line perpendicular to the grating. Instead of a single spectrum, the result is several parallel spectra which are then projected onto a QSI 516M2 CCD camera. The CCD has 9µm square pixels, roughly corresponding to a velocity of 6km/s per pixel. The CCD is cooled to -10◦C to minimize the noise. A lower temperature would be better, but could not be reached due to the CCD limitations. The wavelength ranges of the individual orders overlap to allow flux calibration at the edges. This design of spectrograph is known as echelle (French for staircase) design, due to the

(19)

Figure 3.4: Schematic depiction of a front and side view of the eShel guiding unit. It shows the path of the light inside. Red arrows depict the path of light captured with

the telescope, green arrows the light emitted by the calibration lamps.

raw spectra looking like a ladder or staircase, each order representing a rung. The name eShel is a contraction of this and Shelyak.

3.2 Acquisition

During the nights of observation, τ Boo was visible high in the sky for the entire night. During early observations, τ Boo crossed the southern meridian, and reached the me-chanical limit of the german mount in its daily motion. Therefore, a pier flip was nec-essary during some observations. Dark, bias and flat-field calibrations were done during the day or morning after observations were done. Wavelength-calibrations were done right before and after each individual observation. The software used for autoguiding and dome-control is Maxim DL v5.24. This was the first time at APO autoguiding with the eShel spectrograph was done. The guiding unit of the spectrograph hosts a Watec 120N+ video camera. The guiding unit divides the light between the guiding camera and the entrance fiber using a mirror with a hole. To counteract static friction (stiction) in the motion of the telescope, the Y anti-stiction is set to 1. The X-aggressiveness is set

(20)

Figure 3.5: The orders of an echelle spectrum shown in its translated colours (Cochard

& Thizy 2014).

to 4, Y-aggressiveness to 3. The option for pier-flip correction was used after a pier-flip had been made, but sometimes a recalibration was necessary instead. Exposures were done using the software Audela v2.0.0. This software controls the spectrograph, guiding and calibration unit. Audela allows direct feedback of the exposure quality. After an exposure is completed, Audela immediately renders an image of the obtained spectrum. While this image can be visually inspected for any defects in the spectrum, some tools are available to inspect the quality of the spectrum. A useful tool to determine the pixel saturation of the spectrum is the section tool, which enables the drawing of a line onto the image. This renders a graph of a cross cut through the image along this line, measured in counts. We aimed to stay below a maximum saturation per pixel of 60000 counts to avoid saturation. To acquire a clear spectrum with a sufficient S/N we aimed for a maximum 40000 counts per pixel.

(21)

Table 3.1: Journal of observations

Date Time in UT Exp. time S/N circumstances

20-04-2014 Telescope out of focus.

02-05-2014 Problems with Video Camera

and USB ports.

16-05-2014 23:17 1200 59 17-05-2014 23:46 2400 65 17-05-2014 01:53 2400 49 17-05-2014 02:37 2400 61 17-05-2014 03:18 2400 63 20-05-2014 00:17 2400 49 20-05-2014 01:49 2400 43 20-05-2014 02:46 2400 49 20-05-2014 03:43 2400 60

30-05-2014 00:28 2400 Calibration unit had to be

30-05-2014 01:40 3000 replaced because Th-Ar lamp

30-05-2014 02:49 3600 was working improperly.

01-06-2014 00:01 3600 15 01-06-2014 02:00 3600 22 02-06-2014 00:02 4000 43 02-06-2014 01:56 4000 29 12-06-2014 23:52 4000 27 12-06-2014 01:20 4600 67 12-06-2014 02:45 4600 13

(22)

Chapter 4

Analysis

4.1 Calibrations

The software we used for processing eShel spectra is Audela v2.0.0, the same as used for the acquisition. To transform a raw image into a digitized spectrum, the software needs several calibration images. These are called dark, flatfield, bias and wavelength calibration images. A dark calibration is an exposure with a duration equal to that of the observation, but without any light entering the CCD. The resulting image contains the same average level of thermal and electronic noise (readout noise of the CCD) as the science exposure and is subtracted from the science image. A simple subtraction will multiply the intrinsic noise by a factor√2. Using the average of more dark frames reduces the final noise incrementally. However, a single dark frame was used in our case. A flatfield calibration is an exposure of a led and tungsten lamp giving a more or less uniform spectrum over the whole wavelength range. After normalization to unity, any difference in sensitivity between the pixels of the CCD can be detected and compensated for in this way. Like in the case with the dark frames, dividing by a flatfield frame adds to the noise. Therefore, 5 flatfield frames were averaged to obtain a final flatfield image. The bias calibration is a zero-second exposure to calibrate the zero point of the final fluxes and is only used when dark exposures are not of the same duration as the science exposures. The wavelength calibration is an exposure of a Thorium-Argon (Th-Ar) lamp, of which the wavelengths of the emission lines are well known. By carefully analyzing spectral lines in this image, this calibration determines the relation between pixels and wavelength. Its accuracy is obviously vital for obtaining accurate radial velocity measurements, as this would shift the spectral lines an unknown value in an unknown direction. Accuracy is independent of stability. Spectral calibrations taken at different epochs where the physical condition of the instrumental setup is different

(23)

Chapter 4. Analysis 15

may give equally accurate but different results. It is therefore very important that the circumstances under which the wavelength calibrations take place are therefore as identical as possible to those of the observations. It is not the calibration, it is the temperature of the eShel spectrograph itself that causes mechanical changes leading to changes in the calibration. In practice, wavelength frames are taken before and after each observation to ensure the physical conditions are as similar as possible. Both are then processed by Audela, which takes an average of these frames. Audela then gives as output 17 fits files, each one containing the part of the spectrum of the corresponding order.

4.2 Radial Velocities

To further process the data, a program was written in Wolfram Mathematica 10 to apply heliocentric corrections and perform a cross-correlation. To determine the radial velocity with respect to a standard, the relative shift in wavelength between the two spectra is considered. This shift must have an accuracy that corresponds to a velocity smaller than 100 m/s, which roughly equates to 1/60th of a pixel. This shift can be determined using the cross-correlation method. The cross-correlation of functions f (t) and g(t) is calculated with (Weisstein 2014b):

(f ? g)(τ ) = Z

f∗(t)g(t + τ )dt (4.1)

Here, f ∗ (t) is the complex conjugate of f (t) and τ represents the shift. In this con-text this process overlaps the spectra and then shifts them across each other on the wavelength axis which is sampled in pixels. At each successive shift, both functions are multiplied. The resulting function is then normalized to unity at maximum zeroshift. This gives the cross-correlation function, which has a value between zero and one. The cross-correlation function peaks at the shift where both spectra are the most similar, as two overlapping positive (or negative) peaks result in a high value for the cross-correlation function. Hence the maximum of the cross-cross-correlation function gives the relative shift of the spectra. The peak of the cross-correlation function does not have to be at a shift equal to a complete pixel. Therefore, this method provides the necessary sub-pixel accuracy. This works especially well for spectra with many absorption lines, as this results in a narrower peak of the cross-correlation function, of which the top can be determined more accurately. To calculate the cross-correlation, a method can be used that does not require solving the integral of eq. 4.1. The cross-correlation of f (t) and

(24)

g(t) is equivalent to the convolution of f∗(−t) and g(t).

f (t) ? g(t) = f∗(−t) ∗ g(t) (4.2)

The convolution theorem shows that the convolution of two functions is equal to the product of the Fourier transforms (F ) of both (Weisstein 2014a).

F (f (t) ∗ g(t)) = F (f (t))F (g(t)) (4.3) For a cross-correlation, one can multiply the Fourier transform of the first spectrum with the complex conjugate of the Fourier transform of the second spectrum, and take the inverse Fourier transform afterwards. Using fast Fourier algorithms, this method is suited for a numerical approach to cross-correlation. Additionally, a Gaussian (or any other suitable function) can be fitted to the cross-correlation function to determine the location (in pixels) at its maximum. The accuracy of this fit is used to determine the errors of the found shift.

Figure 4.1: This figure shows two functions in the upper graph, and the cross-correlation function in the lower. The values on the axis are arbitrary. The green and blue functions are identical, with the exception that the green lines has been shifted 9 units to the right. The height of the cross-correlation function corresponds to mul-tiplication of both functions for each value the blue spectrum is shifted to the right. When the blue function is shifted 9 units, both functions completely overlap. The cross-correlation function peaks at this value. In this example the two functions are

identical, in which case it is called the autocorrelation function.

For each spectrum, all 17 orders were cross-correlated with those of the comparison spectrum. The shift (xi) and error (σi) of each of these orders (N) was averaged using

(25)

the weighted average (xwav) and standard deviation (σwav) from Taylor(1982):

xwav= Σwixi Σwi with wi= 1 σ2 i (4.4) σwav = σi √ N (4.5)

4.3 Temperature Dependence

Mechanical stability of the instrument, which limits its final stability, is determined by temperature and air pressure conditions, of which the former is the most important. It is also determined by the direction in which the light enters the the fiber, because the scrambling for cylindrical fibres is not perfect. We consider here only the temperature dependence. The theoretical limit of the accuracy of the eShel spectrograph with a temperature stabilization of 0.05◦C is 50 m/s (Bouchy et al. 2001).

To determine the mechanical changes of the eShel spectrograph as a function of tem-perature, three reference calibrations were taken at different times, at different temper-atures, determined with a thermometer placed at the outside of the eShel box. This box shields the spectrograph from stray light that may possibly disturb any measurements. Temperature measurements were done before and after calibrations were made. These calibrations consist of bias, flatfield and wavelength calibration. Dark calibrations were taken separately without temperature measurements after observations were done. The calibrations were done at average temperatures of 21.0, 21.5 and 22.2◦C. A τ Boo spec-trum from 01-06-2013 from the APO archive was processed three times using each of these calibrations once. This resulted in three different processed spectra, of which the first was cross-correlated with the subsequent to determine the shift in radial velocity as a result of the change in temperatures. To determine the thermal stability of the spectrograph, the four observations done on 17-05-2014 were used. Before and after these observations, the temperature of the spectrograph was measured with the ther-mometer placed on the outside of the eShel box. The temperature slightly cooled over night, with the first observation done at an average temperature of 24.5◦C. The second and third observations were both done at a temperature of 24.0◦C, and the fourth at 23.5◦C, for a total decrease of 1◦C. These spectra were then processed as usual with the radial velocity method, to see whether the observations would yield different results as a consequence of the change in temperature.

(26)

Figure 4.2: This is a raw image of the spectrum used. Dark vertical lines in the white bands indicate absorption lines. The dark lines in the upper-right corner are a clear

(27)

Chapter 5

Results

5.1 Radial Velocity Measurements

Figure 5.1shows the radial velocity measurements of τ Boo (17-05-2014) with errors of 1 standard deviation, and red dots showing the expected values. Also shown are the relative temperatures in◦C. 5.40 5.45 5.50 5.55 5.60 5.65 5.70 -0.4 -0.2 0.0 0.2 0.4 HJD-2456790 V e lo ci ty (km /s )

Boo 17-18 May 2014

T=0.25 T=0.25 T=0.5 T=0

Figure 5.1: Radial velocity measurements of τ Boo (17-05-2014), with red dots show-ing the expected values. Also shown are the relative temperatures in ◦C. The shown

errors are 1σ, with a largest error of 250 m/s.

(28)

Chapter 5. Results 20

Figure5.2shows the radial velocities of all measurements done from 17-05 to 12-06 2014, with errors of 1σ, the largest being 1.3 km/s. The expected radial velocity curve is also shown as a red line.

0.0 0.2 0.4 0.6 0.8 1.0 -2 -1 0 1 2 3 4 Phase Velocity (km /s ) τ Boo 2014

Figure 5.2: Radial velocity measurements of τ Boo (17-05 to 12-06 2014). The red line is the expected radial velocity curve. The shown errors are 1σ, with a largest error

of 1.3 km/s.

5.2 Mechanical Temperature Dependence

Figure 5.3 shows the velocities that correspond to the shift in pixels as a result of the mechanical change induced by the difference in temperatures. The first is an autocor-relation of the first spectrum (T=21.0◦C), which obviously has a shift of 0. The second spectrum (T=21.5◦C) has a shift of -0.29±0.43. The third spectrum (T=22.2◦C) has a shift of -1.46±0.53. Important to note is that these are the averages of the temperatures before and after each exposure. An average change of 0.63±0.2◦ was noticed for each measurement.

(29)

Chapter 5. Results 21 !" #!$ -3.0 -2.5 -2.0 -1.5 -1.0 -0.5 0.0 0.5 %&'()* V e lo ci ty (km /s ) Boo

Figure 5.3: Mechanical change as a function of temperature. The shift in pixels expressed as (radial) velocity in km/s. The shown errors are 1σ, with a largest error of

(30)

Chapter 6

Discussion

6.1 Temperature Dependence

Figure 5.1 shows that all of the measured values lie very close to the expected values. The errors however, are larger than 250 m/s and probably have an instrumental cause. The expected temperature effect lies within the error range, and cannot be determined with this data. The failure of Boots & Meinster (2012) and Barendse & van Lierop

(2013) to successfully determine the radial velocity curve of υ And might be caused by the different method they used to determine the wavelength shift and the correspond-ing errors. Instead of cross-correlatcorrespond-ing suborders of each spectrum and determincorrespond-ing the errors by fitting a Gaussian, they cross-correlated complete orders and determined the errors by averaging and calculating the standard deviation from the shifts of each or-der. However, figure5.2shows that determining the radial velocity curve of τ Boo was unsuccessful as well. This is most likely caused by the spectra being out of focus for both this research and that of Barendse & van Lierop (2013) as shown in figure 6.1. This could also explain the large errors of Boots & Meinster (2012), but none of their spectra could be obtained for comparison. This also explains well the (increasingly) long exposure times needed to obtain a spectrum with desired S/N. There could be several causes for the long exposure times, such as that the optical fiber used to transport light from the telescope to the spectrograph is damaged or faulty or that while the stellar image was in focus on the camera, it was not in focus on the fiber, as shown in figure

6.2. The most probable cause is that the fiber entering the eShel is not properly focused on the collimator. This would reduce the resolution, increase the necessary exposure times, and would explain the lack of focus of the spectra. Additionally, to better deter-mine the thermal stability of the spectrograph, a star without a radial velocity should be observed. The difference in wavelength shift between two measurements caused by

(31)

Chapter 6. Discussion 23

thermal instability is obscured by any extra shift in wavelength caused by motion of the star.

Figure 6.1: Three similar parts of spectra from van Woerkom (2011), Barendse & van Lierop(2013) and this research. Only the first is properly in focus, as seen by the

blurriness of the (dark) spectral lines in the white bands.

Figure 6.2: The light path inside the guiding unit. The placement of the components are not to scale and exaggerated. The ideal situation is option number 1. The light entering the guiding unit (light yellow) is focused precisely on the fiber, with only a minimal loss of flux due to a small part reflecting off the mirror to the guiding camera. The second situation, where the distance between the mirror and the eShel fiber is increased, would appear exactly as the first situation on the guiding camera, but a lot of the light for the spectrograph would be lost because of the dispersion, significantly

diminishing the flux.

6.2 Mechanical Temperature Stability

Figure 5.3 shows that the wavelength shift expressed in velocity is -0.29±0.43 km/s at a ∆T of 0.5◦C and -1.46±0.53 km/s at a ∆T of 1.2◦C. This shows that temperature differences of at least 1◦C have a considerable effect on the structure of the spectrograph.

(32)

It is therefore important that the circumstances under which the wavelength calibration and the science exposures occur must be kept as similar as possible. A similar research was conducted by Garde (2015), who studied the direct shift in pixels of the spectral lines of a Th-Ar spectrum. He used an ATIK 460EX camera with the CCD-chip cooled to a temperature of -5◦C, and two PT 100 thermometers to measure both the ambient temperature and the temperature of the metal exterior of the spectrograph. Both were connected to a PID which registered the temperatures each minute of both electronic thermometers with an USB-connection. For one and a half hours, Th-Ar exposures of 5 seconds each are taken continuously, with images taken every three minutes, resulting in a total of 30 images. Using the Mira Pro software, he measured the centrepoint shift of hydrogen lines in orders 34 (H-α), 46 (H-β) and 52 (H-γ) of the spectrum. Figure

6.3 shows the plot of X and Y pixelshift and temperature as a function of time, of the 46th order. The same plots of other orders look similar. The spectrograph temperature rises steeply the fist 30 minutes, and rises about 3◦C in total. This is probably due to the cooler of the CCD producing heat during the initial cooling of the CCD. A similar increase in temperature was noticed during the wavelength calibrations done at APO. Figure6.4is a table showing the maximum pixelshift at a temperature difference of 1.04 and 0.56◦C for both the X and Y direction. For the wavelength calibration, only the X direction is important. This shows that the average shift is 0.26 pixels at a ∆T of 0.5◦C and 0.38 pixels at a ∆T of 1◦C, which corresponds to 1.6 and 2.3 km/s respectively. These are similar results, confirming that the wavelength calibration is very sensitive to temperature.

Figure 6.3: Graph showing the temperature (left axis) of the environment (blue) and the spectrograph (yellow) and the shift in pixels (right axis) in both the X (green) and

(33)

Figure 6.4: The pixelshift in the X and Y direction of all three orders at two different spectrograph temperatures. Adapted from (Garde 2015)

(34)

Chapter 7

Conclusion

Temperature has a large effect on the wavelength calibration of the fiber-fed eShel spec-trograph, as the mechanical expansion and contraction of the spectrograph changes the lightpath inside. This is very important for radial velocity measurements of exoplanet orbits, where an accuracy smaller than 100 m/s is desired. The wavelength shift induced by changing temperatures, expressed in velocities is -0.29±0.43 km/s at a ∆T of 0.5◦C and -1.46±0.53 km/s at a ∆T of 1.2◦C. Similar results were achieved by Garde (2015), who found shifts of ∆T of 1.6 and 2.3 km/s for a ∆T of 0.5 and 1◦C respectively, by considering the physical shift of spectral lines on the CCD as a function of temperature. This effect can be calibrated for by taking wavelength calibrations before and after each science exposure. However, different exposures under similar conditions may yield dif-ferent results despite a high accuracy. This stability effect was investigated as well by doing radial velocity measurements of τ Boo, but any temperature effect fell within the error range. This was caused by a lack of accuracy due to the eShel spectra being out of focus, probably because the fiber is not properly focused on the eShel collimator. However, this might have been the case for previous radial velocity measurements at APO by Boots & Meinster (2012) andBarendse & van Lierop (2013) and explains the large inaccuracy of their results as well. The temperature stability of the eShel spectro-graph could be investigated further with more accurate observations during which the temperature of the spectrograph is measured carefully, preferably using an electronic thermometer. Additionally, a radial velocity standard star would be better for this type of research as the difference in wavelength shift between two measurements caused by thermal instability is obscured by any extra shift in wavelength caused by motion of the star. Temperature fluctuations with an average of 0.63±0.2◦ were noticed during wavelength calibrations. This is probably because the CCD-cooler releases heat while cooling the CCD, and this heat is allowed to build up inside the box, instead of being vented into the room. Figure7.1 shows a suggested improvement to the eShel setup to

(35)

Chapter 7. Conclusion 27

reduce this problem. By lowering the top of the box and allowing the cooler to stick out of a hole, the spectrograph is protected from these temperature fluctuations.

Figure 7.1: If CCD and cooler stick out of the box containing the eShel, any heat produced by the cooling process is vented into the room instead of the box, reducing

(36)

Chapter 7. Conclusion 28

I would like to thank Samayra Straal, Huib Henrichs, Abigail Stevens and Jasmijn van Vulpen for their guidance, help and encouragement throughout this endeavor, which wouldn’t have been successful without.

(37)

Appendix A

Reduction Program

For further reduction, a program was written in Mathematica 10.1 that uses a crosscor-relation function to correlate the shifted spectra. It imports the fits-files of the spectra as well as a file containing the heliocentric correction and expected phase and radial velocity corresponding to the spectra. In the program a number of variables must be manually set for the cross-correlation to function. These are the number of spectra to be correlated and the number of orders in each spectra, the number of suborders each order must be divided into, and the number of pixels the spectra must be shifted (-6 to 6 is used). The heliocentric correction can also be calculated with the program. Each oder of each spectrum is subsequently divided into 4 suborders of even length. The first spectrum if then used as reference spectrum, and each suborder of this spectrum is then cross-correlated with the corresponding suborder of all other spectra. An 8th order polynomial with even orders is then fitted to each of the resulting cross-correlation function, of which the maximum is determined, as well as the accuracy of this maximum. These values are then converted to velocities in km/s. Next, the (weigthed) average shift and standard deviation (σ) are calculated. All suborders with a σ less than a set limit (0.05) are discarded. The data is further selected by removing all points of which the absolute value with the mean subtracted is larger than 2σ. Of the remaining points, the mean and standard deviation are calculated for each cross-correlation and subsequently plotted.

(38)

Appendix B

Reduction Program Code

1)1)1) Clear the memory, initialize a package for the use of error bars, and select a directoryClear the memory, initialize a package for the use of error bars, and select a directoryClear the memory, initialize a package for the use of error bars, and select a directory Quit

QuitQuit

SetSystemOptions[“DataOptions” → “ReturnQuantities” → False]; (*forHJD : newinversion10*) SetSystemOptions[“DataOptions” → “ReturnQuantities” → False]; (*forHJD : newinversion10*)SetSystemOptions[“DataOptions” → “ReturnQuantities” → False]; (*forHJD : newinversion10*)

Needs[“ErrorBarPlots`”] Needs[“ErrorBarPlots`”]Needs[“ErrorBarPlots`”] clight = 299792.458; clight = 299792.458;clight = 299792.458; SetDirectory[“Directory”] SetDirectory[“Directory”]SetDirectory[“Directory”]

2)2)2) Defining the variables, tables and functions:Defining the variables, tables and functions:Defining the variables, tables and functions:

noo Number of orders

nos Number of spectra

split The number of parts the orders are going to be split up in

np Thee number of pixels the spectrum is to be shifted left and right during correlation

fns A table containing the filenames in the directory. We select one using “fns[[1,Spectrum##,Order##]]”

plots The table that will be filled with the correlation graphs shift The table that will be filled with the velocity shifts graphs

error The table that will be filled with the errors for the velocity shifts graphs HJDmidExpo The table that will be filled with the mid-exposure Heliocentric Ju-lian Date for each spectrum

vhelio The list of Heliocentric Velocity Corrections, imported from helio.txt in the selected directory

(39)

Appendix B. Reduction Program Code 31

raList The list of Right Ascension coordinates of τ -Boo decList The list of Declination coordinates of τ -Boo

phase The list containing the phase of τ -Boo during the recording of each spec-trum

key A function used for extracting variables from the imported datafiles JulianDate A function for calculating the Julian Date

ModifiedJulianDate A function for calculating the Modified Julian Date

HeliocentricCorrectionTime A function for calculating the heliocentric correction of the Julian Date

HeliocentricJulianDate A function for calculating the Heliocentric Julian Date noo = 17;noo = 17;noo = 17;

nos = 13;nos = 13;nos = 13; split = 4;split = 4;split = 4; np = 6;np = 6;np = 6;

fns = FileNames[“*.fit”];fns = FileNames[“*.fit”];fns = FileNames[“*.fit”]; fns//TableFormfns//TableFormfns//TableForm

plots = Table[0, {i, noo ∗ split}, {k, nos}];plots = Table[0, {i, noo ∗ split}, {k, nos}];plots = Table[0, {i, noo ∗ split}, {k, nos}]; shift = Table[0, {i, noo ∗ split}, {k, nos}];shift = Table[0, {i, noo ∗ split}, {k, nos}];shift = Table[0, {i, noo ∗ split}, {k, nos}]; error = Table[0, {i, noo ∗ split}, {k, nos}];error = Table[0, {i, noo ∗ split}, {k, nos}];error = Table[0, {i, noo ∗ split}, {k, nos}]; HJDmidExpo = Table[0, {k, nos}];HJDmidExpo = Table[0, {k, nos}];HJDmidExpo = Table[0, {k, nos}];

Print[“vhelio = ”, vhelio = Drop[Import[“vhelio phase.txt”, “Table”], 1][[All, 3]], “ km/s”]Print[“vhelio = ”, vhelio = Drop[Import[“vhelio phase.txt”, “Table”], 1][[All, 3]], “ km/s”]Print[“vhelio = ”, vhelio = Drop[Import[“vhelio phase.txt”, “Table”], 1][[All, 3]], “ km/s”] Print[“phase = ”, phase = Drop[Import[“vhelio phase.txt”, “Table”], 1][[All, 4]]]Print[“phase = ”, phase = Drop[Import[“vhelio phase.txt”, “Table”], 1][[All, 4]]]Print[“phase = ”, phase = Drop[Import[“vhelio phase.txt”, “Table”], 1][[All, 4]]]

Print[“vrad = ”, vrad = Drop[Import[“vhelio phase.txt”, “Table”], 1][[All, 5]]/1000., “ km/s”]Print[“vrad = ”, vrad = Drop[Import[“vhelio phase.txt”, “Table”], 1][[All, 5]]/1000., “ km/s”]Print[“vrad = ”, vrad = Drop[Import[“vhelio phase.txt”, “Table”], 1][[All, 5]]/1000., “ km/s”]

raList = {13, 47, 15.7};raList = {13, 47, 15.7};raList = {13, 47, 15.7}; decList = {17, 27, 24.8};decList = {17, 27, 24.8};decList = {17, 27, 24.8};

key[dat , keyword ]:=(keypos = Position[dat, keyword][[1]]; ++keypos[[−1]];key[dat , keyword ]:=(keypos = Position[dat, keyword][[1]]; ++keypos[[−1]];key[dat , keyword ]:=(keypos = Position[dat, keyword][[1]]; ++keypos[[−1]]; If[Length[Extract[dat, keypos]] == 0,If[Length[Extract[dat, keypos]] == 0,If[Length[Extract[dat, keypos]] == 0,

Extract[dat, keypos], Extract[dat, keypos][[1]]]);Extract[dat, keypos], Extract[dat, keypos][[1]]]);Extract[dat, keypos], Extract[dat, keypos][[1]]]);

JulianDate[{y, m, d, h, mint, sec}]:=367y − IntegerPart[7(y + IntegerPart[(m + 9)/12])/4]−JulianDate[{y, m, d, h, mint, sec}]:=367y − IntegerPart[7(y + IntegerPart[(m + 9)/12])/4]−JulianDate[{y, m, d, h, mint, sec}]:=367y − IntegerPart[7(y + IntegerPart[(m + 9)/12])/4]−

(40)

IntegerPart[275m/9] + d + 1721028.5 + (h + mint/60 + sec /3600)/24; IntegerPart[275m/9] + d + 1721028.5 + (h + mint/60 + sec /3600)/24;IntegerPart[275m/9] + d + 1721028.5 + (h + mint/60 + sec /3600)/24;

ModifiedJulianDate[{y , m , d , h , min , sec }]:=JulianDate[{y, m, d, h, min, sec}] − 2400000.5; ModifiedJulianDate[{y , m , d , h , min , sec }]:=JulianDate[{y, m, d, h, min, sec}] − 2400000.5;ModifiedJulianDate[{y , m , d , h , min , sec }]:=JulianDate[{y, m, d, h, min, sec}] − 2400000.5; HeliocentricCorrectionTime[date List, ra List, dec List]:=

HeliocentricCorrectionTime[date List, ra List, dec List]:=HeliocentricCorrectionTime[date List, ra List, dec List]:=

Module[{ALPHA, DELTA, PICON, TT, P, RL, G, AJ, GLM, GLP, TGLP, TGLM, X, Y, A, DEL}, Module[{ALPHA, DELTA, PICON, TT, P, RL, G, AJ, GLM, GLP, TGLP, TGLM, X, Y, A, DEL},Module[{ALPHA, DELTA, PICON, TT, P, RL, G, AJ, GLM, GLP, TGLP, TGLM, X, Y, A, DEL}, ALPHA = ra[[1]] + ra[[2]]/60 + ra[[3]]/3600;

ALPHA = ra[[1]] + ra[[2]]/60 + ra[[3]]/3600;ALPHA = ra[[1]] + ra[[2]]/60 + ra[[3]]/3600; DELTA = dec[[1]] + dec[[2]]/60 + dec[[3]]/3600; DELTA = dec[[1]] + dec[[2]]/60 + dec[[3]]/3600;DELTA = dec[[1]] + dec[[2]]/60 + dec[[3]]/3600; PICON = 2Pi/360;

PICON = 2Pi/360;PICON = 2Pi/360;

TT = (JulianDate[date] − 2400000 − 15020.)/36525.; TT = (JulianDate[date] − 2400000 − 15020.)/36525.;TT = (JulianDate[date] − 2400000 − 15020.)/36525.; P = (1.396041 + 0.000308 ∗ (TT + 0.5)) ∗ (TT − 0.499998); P = (1.396041 + 0.000308 ∗ (TT + 0.5)) ∗ (TT − 0.499998);P = (1.396041 + 0.000308 ∗ (TT + 0.5)) ∗ (TT − 0.499998); RL = 279.696678 + 36000.76892 ∗ TT + 0.000303 ∗ TT∧2 − P ; RL = 279.696678 + 36000.76892 ∗ TT + 0.000303 ∗ TTRL = 279.696678 + 36000.76892 ∗ TT + 0.000303 ∗ TT∧∧2 − P ;2 − P ; G = 358.475833 + 35999.04975 ∗ TT − 0.00015 ∗ TT∧2; G = 358.475833 + 35999.04975 ∗ TT − 0.00015 ∗ TTG = 358.475833 + 35999.04975 ∗ TT − 0.00015 ∗ TT∧∧2;2; AJ = 225.444651 + 2880.0 ∗ TT + 154.906654 ∗ TT; AJ = 225.444651 + 2880.0 ∗ TT + 154.906654 ∗ TT;AJ = 225.444651 + 2880.0 ∗ TT + 154.906654 ∗ TT; If[G > 360, G-=360]; If[G > 360, G-=360];If[G > 360, G-=360]; If[RL > 360, RL-=360]; If[RL > 360, RL-=360];If[RL > 360, RL-=360]; RL = RL ∗ PICON; RL = RL ∗ PICON;RL = RL ∗ PICON; G = G ∗ PICON; G = G ∗ PICON;G = G ∗ PICON; AJ = AJ ∗ PICON; AJ = AJ ∗ PICON;AJ = AJ ∗ PICON; GLM = G − RL; GLM = G − RL;GLM = G − RL; GLP = G + RL; GLP = G + RL;GLP = G + RL; TGLP = 2.0 ∗ G + RL; TGLP = 2.0 ∗ G + RL;TGLP = 2.0 ∗ G + RL; TGLM = 2.0 ∗ G − RL; TGLM = 2.0 ∗ G − RL;TGLM = 2.0 ∗ G − RL;

X = .99986 ∗ Cos[RL] − 0.025127 ∗ Cos[GLM] + 0.008374 ∗ Cos[GLP] + .000105 ∗ Cos[TGLP]+ X = .99986 ∗ Cos[RL] − 0.025127 ∗ Cos[GLM] + 0.008374 ∗ Cos[GLP] + .000105 ∗ Cos[TGLP]+X = .99986 ∗ Cos[RL] − 0.025127 ∗ Cos[GLM] + 0.008374 ∗ Cos[GLP] + .000105 ∗ Cos[TGLP]+

.000063 ∗ TT ∗ Cos[GLM] + .000035 ∗ Cos[TGLM] − .000026 ∗ Sin[GLM − AJ] − .000021 ∗ TT ∗ Cos[GLP]; .000063 ∗ TT ∗ Cos[GLM] + .000035 ∗ Cos[TGLM] − .000026 ∗ Sin[GLM − AJ] − .000021 ∗ TT ∗ Cos[GLP];.000063 ∗ TT ∗ Cos[GLM] + .000035 ∗ Cos[TGLM] − .000026 ∗ Sin[GLM − AJ] − .000021 ∗ TT ∗ Cos[GLP]; Y = .917308 ∗ Sin[RL] + .023053 ∗ Sin[GLM] + .007683 ∗ Sin[GLP] + .000097 ∗ Sin[TGLP]−

Y = .917308 ∗ Sin[RL] + .023053 ∗ Sin[GLM] + .007683 ∗ Sin[GLP] + .000097 ∗ Sin[TGLP]−Y = .917308 ∗ Sin[RL] + .023053 ∗ Sin[GLM] + .007683 ∗ Sin[GLP] + .000097 ∗ Sin[TGLP]−

.000057 ∗ TT ∗ Sin[GLM] − .000032 ∗ Sin[TGLM] − .000024 ∗ Cos[GLM − AJ] − .000019 ∗ TT ∗ Sin[GLP]; .000057 ∗ TT ∗ Sin[GLM] − .000032 ∗ Sin[TGLM] − .000024 ∗ Cos[GLM − AJ] − .000019 ∗ TT ∗ Sin[GLP];.000057 ∗ TT ∗ Sin[GLM] − .000032 ∗ Sin[TGLM] − .000024 ∗ Cos[GLM − AJ] − .000019 ∗ TT ∗ Sin[GLP]; A = ALPHA ∗ 15.0 ∗ PICON;

A = ALPHA ∗ 15.0 ∗ PICON;A = ALPHA ∗ 15.0 ∗ PICON; DEL = DELTA ∗ PICON; DEL = DELTA ∗ PICON;DEL = DELTA ∗ PICON;

−0.0057755 ∗ ((Cos[DEL] ∗ Cos[A]) ∗ X + (.4337751 ∗ Sin[DEL] + Cos[DEL] ∗ Sin[A]) ∗ Y )]; −0.0057755 ∗ ((Cos[DEL] ∗ Cos[A]) ∗ X + (.4337751 ∗ Sin[DEL] + Cos[DEL] ∗ Sin[A]) ∗ Y )];−0.0057755 ∗ ((Cos[DEL] ∗ Cos[A]) ∗ X + (.4337751 ∗ Sin[DEL] + Cos[DEL] ∗ Sin[A]) ∗ Y )];

(41)

HeliocentricJulianDate[date List, ra List, dec List]:=HeliocentricJulianDate[date List, ra List, dec List]:=HeliocentricJulianDate[date List, ra List, dec List]:=

JulianDate[date] + HeliocentricCorrectionTime[date, ra, dec];JulianDate[date] + HeliocentricCorrectionTime[date, ra, dec];JulianDate[date] + HeliocentricCorrectionTime[date, ra, dec]; 3)

3)3) For each spectrumFor each spectrumFor each spectrum k, kk, kk, k =1,19, order=1,19, order=1,19, order i, ii, ii, i = 1,17= 1,17= 1,17 and fraction of each orderand fraction of each orderand fraction of each order j, jj, jj, j = 0,= 0,= 0, splitsplitsplit calculate the correlation with respect to the spectrum in fitdata2 as a function of the shift, and find the maximum of the correlation.calculate the correlation with respect to the spectrum in fitdata2 as a function of the shift, and find the maximum of the correlation.calculate the correlation with respect to the spectrum in fitdata2 as a function of the shift, and find the maximum of the correlation. Monitor[Monitor[Monitor[

Do[Do[Do[ Do[Do[Do[

fitdata1 = Import[fns[[k]], {{“Metadata”, “Data”}}];fitdata1 = Import[fns[[k]], {{“Metadata”, “Data”}}];fitdata1 = Import[fns[[k]], {{“Metadata”, “Data”}}]; fitdata2 = Import[fns[[1]], {{“Metadata”, “Data”}}];fitdata2 = Import[fns[[1]], {{“Metadata”, “Data”}}];fitdata2 = Import[fns[[1]], {{“Metadata”, “Data”}}]; object1 = key[fitdata1, “OBJNAME”];object1 = key[fitdata1, “OBJNAME”];object1 = key[fitdata1, “OBJNAME”];

orderID1 = key[fitdata1[[1, 36 + i]], “EXTNAME”];orderID1 = key[fitdata1[[1, 36 + i]], “EXTNAME”];orderID1 = key[fitdata1[[1, 36 + i]], “EXTNAME”]; seriesID1 = key[fitdata1, “SERIESID”];seriesID1 = key[fitdata1, “SERIESID”];seriesID1 = key[fitdata1, “SERIESID”];

expoTime1 = key[fitdata1, “EXPOSURE”];expoTime1 = key[fitdata1, “EXPOSURE”];expoTime1 = key[fitdata1, “EXPOSURE”]; expoStart1 = key[fitdata1, “DATE-OBS”];expoStart1 = key[fitdata1, “DATE-OBS”];expoStart1 = key[fitdata1, “DATE-OBS”]; expoEnd1 = key[fitdata1, “DATE-END”];expoEnd1 = key[fitdata1, “DATE-END”];expoEnd1 = key[fitdata1, “DATE-END”]; wav11 = key[fitdata1[[1, 36 + i]], “CRVAL1”];wav11 = key[fitdata1[[1, 36 + i]], “CRVAL1”];wav11 = key[fitdata1[[1, 36 + i]], “CRVAL1”]; disp1 = key[fitdata1[[1, 36 + i]], “CDELT1”];disp1 = key[fitdata1[[1, 36 + i]], “CDELT1”];disp1 = key[fitdata1[[1, 36 + i]], “CDELT1”]; npixels1 = key[fitdata1[[1, 36 + i]], “NAXIS1”];npixels1 = key[fitdata1[[1, 36 + i]], “NAXIS1”];npixels1 = key[fitdata1[[1, 36 + i]], “NAXIS1”];

object2 = key[fitdata2, “OBJNAME”];object2 = key[fitdata2, “OBJNAME”];object2 = key[fitdata2, “OBJNAME”];

orderID2 = key[fitdata2[[1, 36 + i]], “EXTNAME”];orderID2 = key[fitdata2[[1, 36 + i]], “EXTNAME”];orderID2 = key[fitdata2[[1, 36 + i]], “EXTNAME”]; seriesID2 = key[fitdata2, “SERIESID”];seriesID2 = key[fitdata2, “SERIESID”];seriesID2 = key[fitdata2, “SERIESID”];

expoTime2 = key[fitdata2, “EXPOSURE”];expoTime2 = key[fitdata2, “EXPOSURE”];expoTime2 = key[fitdata2, “EXPOSURE”]; expoStart2 = key[fitdata2, “DATE-OBS”];expoStart2 = key[fitdata2, “DATE-OBS”];expoStart2 = key[fitdata2, “DATE-OBS”]; expoEnd2 = key[fitdata2, “DATE-END”];expoEnd2 = key[fitdata2, “DATE-END”];expoEnd2 = key[fitdata2, “DATE-END”]; wav12 = key[fitdata2[[1, 36 + i]], “CRVAL1”];wav12 = key[fitdata2[[1, 36 + i]], “CRVAL1”];wav12 = key[fitdata2[[1, 36 + i]], “CRVAL1”]; disp2 = key[fitdata2[[1, 36 + i]], “CDELT1”];disp2 = key[fitdata2[[1, 36 + i]], “CDELT1”];disp2 = key[fitdata2[[1, 36 + i]], “CDELT1”]; npixels2 = key[fitdata2[[1, 36 + i]], “NAXIS1”];npixels2 = key[fitdata2[[1, 36 + i]], “NAXIS1”];npixels2 = key[fitdata2[[1, 36 + i]], “NAXIS1”]; maxwav = Max[wav11, wav12];maxwav = Max[wav11, wav12];maxwav = Max[wav11, wav12];

fdata1 = fitdata1[[2, 36 + i]];fdata1 = fitdata1[[2, 36 + i]];fdata1 = fitdata1[[2, 36 + i]]; fdata2 = fitdata2[[2, 36 + i]];fdata2 = fitdata2[[2, 36 + i]];fdata2 = fitdata2[[2, 36 + i]];

(42)

dropwavdata1 = Drop[fdata1, 10 ∗ (maxwav − wav11)]; dropwavdata1 = Drop[fdata1, 10 ∗ (maxwav − wav11)];dropwavdata1 = Drop[fdata1, 10 ∗ (maxwav − wav11)]; dropwavdata2 = Drop[fdata2, 10 ∗ (maxwav − wav12)]; dropwavdata2 = Drop[fdata2, 10 ∗ (maxwav − wav12)];dropwavdata2 = Drop[fdata2, 10 ∗ (maxwav − wav12)]; dropzeros1 = Length[Position[dropwavdata1, 0.]]; dropzeros1 = Length[Position[dropwavdata1, 0.]];dropzeros1 = Length[Position[dropwavdata1, 0.]]; dropzeros2 = Length[Position[dropwavdata2, 0.]]; dropzeros2 = Length[Position[dropwavdata2, 0.]];dropzeros2 = Length[Position[dropwavdata2, 0.]]; ddata1 = Drop[dropwavdata1, −dropzeros1]; ddata1 = Drop[dropwavdata1, −dropzeros1];ddata1 = Drop[dropwavdata1, −dropzeros1]; ddata2 = Drop[dropwavdata2, −dropzeros2]; ddata2 = Drop[dropwavdata2, −dropzeros2];ddata2 = Drop[dropwavdata2, −dropzeros2]; mini = Min[Length[ddata1], Length[ddata2]]; mini = Min[Length[ddata1], Length[ddata2]];mini = Min[Length[ddata1], Length[ddata2]];

data1t = Take[ddata1, mini]; data1t = Take[ddata1, mini];data1t = Take[ddata1, mini]; data2t = Take[ddata2, mini]; data2t = Take[ddata2, mini];data2t = Take[ddata2, mini];

For[j = 0, j < split, j++, For[j = 0, j < split, j++,For[j = 0, j < split, j++,

data1 = Take[data1t, {(split − j − 1) ∗ IntegerPart[Length[data1t]/split] + 1, data1 = Take[data1t, {(split − j − 1) ∗ IntegerPart[Length[data1t]/split] + 1,data1 = Take[data1t, {(split − j − 1) ∗ IntegerPart[Length[data1t]/split] + 1, (split − j) ∗ IntegerPart[Length[data1t]/split]}];

(split − j) ∗ IntegerPart[Length[data1t]/split]}];(split − j) ∗ IntegerPart[Length[data1t]/split]}];

data2 = Take[data2t, {(split − j − 1) ∗ IntegerPart[Length[data2t]/split] + 1, data2 = Take[data2t, {(split − j − 1) ∗ IntegerPart[Length[data2t]/split] + 1,data2 = Take[data2t, {(split − j − 1) ∗ IntegerPart[Length[data2t]/split] + 1, (split − j) ∗ IntegerPart[Length[data2t]/split]}];

(split − j) ∗ IntegerPart[Length[data2t]/split]}];(split − j) ∗ IntegerPart[Length[data2t]/split]}]; data1s[o ]:=RotateRight[data1, o];

data1s[o ]:=RotateRight[data1, o];data1s[o ]:=RotateRight[data1, o]; data2s[o ]:=RotateRight[data2, o]; data2s[o ]:=RotateRight[data2, o];data2s[o ]:=RotateRight[data2, o];

correl = Table[{o, Correlation[data1, data2s[o]]}, {o, −np, np}]; correl = Table[{o, Correlation[data1, data2s[o]]}, {o, −np, np}];correl = Table[{o, Correlation[data1, data2s[o]]}, {o, −np, np}]; nlm = NonlinearModelFit[correl,

nlm = NonlinearModelFit[correl,nlm = NonlinearModelFit[correl,

a + b(x − q)∧2 + c(x − q)∧4 + d(x − q)∧6 + e(x − q)∧8, {a, b, c, d, q, e}, x]; a + b(x − q)a + b(x − q)∧∧2 + c(x − q)2 + c(x − q)∧∧4 + d(x − q)4 + d(x − q)∧∧6 + e(x − q)6 + e(x − q)∧∧8, {a, b, c, d, q, e}, x];8, {a, b, c, d, q, e}, x]; plots[[j + 1 + (i − 1) ∗ split]][[k]] =

plots[[j + 1 + (i − 1) ∗ split]][[k]] =plots[[j + 1 + (i − 1) ∗ split]][[k]] =

Plot[Normal[nlm], {x, −np, np}, ImageSize → Scaled[.2], PlotStyle → Blue, Plot[Normal[nlm], {x, −np, np}, ImageSize → Scaled[.2], PlotStyle → Blue,Plot[Normal[nlm], {x, −np, np}, ImageSize → Scaled[.2], PlotStyle → Blue,

PlotLabel → “spectrum = ” <> ToString[k] <> “, order = ” <> ToString[i] <> “, frac = ” <> ToString[j], PlotLabel → “spectrum = ” <> ToString[k] <> “, order = ” <> ToString[i] <> “, frac = ” <> ToString[j],PlotLabel → “spectrum = ” <> ToString[k] <> “, order = ” <> ToString[i] <> “, frac = ” <> ToString[j], Epilog → {Point[correl]}];

Epilog → {Point[correl]}];Epilog → {Point[correl]}];

shift[[j + 1 + (i − 1) ∗ split]][[k]] = vhelio[[k]] − vhelio[[1]]+ shift[[j + 1 + (i − 1) ∗ split]][[k]] = vhelio[[k]] − vhelio[[1]]+shift[[j + 1 + (i − 1) ∗ split]][[k]] = vhelio[[k]] − vhelio[[1]]+

((clightdisp1)/(maxwav + (j + 1/2)IntegerPart[Length[data1t]/split] ∗ disp1))nlm[[1, 2, 5, 2]]; ((clightdisp1)/(maxwav + (j + 1/2)IntegerPart[Length[data1t]/split] ∗ disp1))nlm[[1, 2, 5, 2]];((clightdisp1)/(maxwav + (j + 1/2)IntegerPart[Length[data1t]/split] ∗ disp1))nlm[[1, 2, 5, 2]]; (*shift1[[i]][[k]] = vhelio[[k]] + ((300000disp1)/(wav11 + disp1npixels1/2))

(*shift1[[i]][[k]] = vhelio[[k]] + ((300000disp1)/(wav11 + disp1npixels1/2))(*shift1[[i]][[k]] = vhelio[[k]] + ((300000disp1)/(wav11 + disp1npixels1/2)) Solve[D[Normal[nlm], x] == 0, x][[3, 1, 2]]; *)

Solve[D[Normal[nlm], x] == 0, x][[3, 1, 2]]; *)Solve[D[Normal[nlm], x] == 0, x][[3, 1, 2]]; *) error[[j + 1 + (i − 1) ∗ split]][[k]] = ((clightdisp1)/ error[[j + 1 + (i − 1) ∗ split]][[k]] = ((clightdisp1)/error[[j + 1 + (i − 1) ∗ split]][[k]] = ((clightdisp1)/

(maxwav + (j + 1/2)IntegerPart[Length[data1t]/split] ∗ disp1))nlm[“ParameterErrors”][[5]]; (maxwav + (j + 1/2)IntegerPart[Length[data1t]/split] ∗ disp1))nlm[“ParameterErrors”][[5]];(maxwav + (j + 1/2)IntegerPart[Length[data1t]/split] ∗ disp1))nlm[“ParameterErrors”][[5]];

(43)

];];];, {i, 1, noo}]; (* end j for loop *), {i, 1, noo}]; (* end j for loop *), {i, 1, noo}]; (* end j for loop *)

HJDmidExpo[[k]] = (HeliocentricJulianDate[DateList[expoStart1], raList, decList]+HJDmidExpo[[k]] = (HeliocentricJulianDate[DateList[expoStart1], raList, decList]+HJDmidExpo[[k]] = (HeliocentricJulianDate[DateList[expoStart1], raList, decList]+ HeliocentricJulianDate[DateList[expoEnd1], raList, decList])/2HeliocentricJulianDate[DateList[expoEnd1], raList, decList])/2HeliocentricJulianDate[DateList[expoEnd1], raList, decList])/2

, {k, 1, nos}], {k, 1, nos}], {k, 1, nos}]

, ProgressIndicator[17 ∗ (k − 1) + (i − 1), {1, 17 ∗ (nos − 1) + (noo − 1) − 1}]], ProgressIndicator[17 ∗ (k − 1) + (i − 1), {1, 17 ∗ (nos − 1) + (noo − 1) − 1}]], ProgressIndicator[17 ∗ (k − 1) + (i − 1), {1, 17 ∗ (nos − 1) + (noo − 1) − 1}]] Print[“Done”];Print[“Done”];Print[“Done”];

4) Plot the calculated velocity as a function of the number of the suborder, for each spectrum 4)4) Plot the calculated velocity as a function of the number of the suborder, for each spectrumPlot the calculated velocity as a function of the number of the suborder, for each spectrum k.k.k. norm[a , b , c , w , x ]:=a + bExph−1

2 c−x

w

2i norm[a , b , c , w , x ]:=a + bExp

h

−1₂ c−x_w 2i norm[a , b , c , w , x ]:=a + bExph−1

2 c−x

w

2i

nlm1 = NonlinearModelFit[correl, norm[a, b, c, w, x], {a, b, c, w}, x]nlm1 = NonlinearModelFit[correl, norm[a, b, c, w, x], {a, b, c, w}, x]nlm1 = NonlinearModelFit[correl, norm[a, b, c, w, x], {a, b, c, w}, x] nlm1[“ParameterErrors”]nlm1[“ParameterErrors”]nlm1[“ParameterErrors”]

nlm1[“BestFitParameters”]nlm1[“BestFitParameters”]nlm1[“BestFitParameters”]

Plot[{Normal[nlm1] − Normal[nlm]},Plot[{Normal[nlm1] − Normal[nlm]},Plot[{Normal[nlm1] − Normal[nlm]},

{x, −np, np}, ImageSize → Scaled[.2], PlotStyle → Blue,{x, −np, np}, ImageSize → Scaled[.2], PlotStyle → Blue,{x, −np, np}, ImageSize → Scaled[.2], PlotStyle → Blue,

PlotLabel → “spectrum = ” <> ToString[k] <> “, order = ” <> ToString[i] <>PlotLabel → “spectrum = ” <> ToString[k] <> “, order = ” <> ToString[i] <>PlotLabel → “spectrum = ” <> ToString[k] <> “, order = ” <> ToString[i] <> “, frac = ” <> ToString[j], Epilog → {Point[correl]}]“, frac = ” <> ToString[j], Epilog → {Point[correl]}]“, frac = ” <> ToString[j], Epilog → {Point[correl]}]

Plot[Normal[nlm], {x, −np, np}, ImageSize → Scaled[.2], PlotStyle → Blue,Plot[Normal[nlm], {x, −np, np}, ImageSize → Scaled[.2], PlotStyle → Blue,Plot[Normal[nlm], {x, −np, np}, ImageSize → Scaled[.2], PlotStyle → Blue, PlotLabel → “spectrum = ” <> ToString[k] <> “, order = ” <>PlotLabel → “spectrum = ” <> ToString[k] <> “, order = ” <>PlotLabel → “spectrum = ” <> ToString[k] <> “, order = ” <>

ToString[i] <> “, frac = ” <> ToString[j], Epilog → {Point[correl]}]ToString[i] <> “, frac = ” <> ToString[j], Epilog → {Point[correl]}]ToString[i] <> “, frac = ” <> ToString[j], Epilog → {Point[correl]}]

((nlm[[1, 2, 5, 2]] − nlm1[[1, 2, 3, 2]])/nlm[[1, 2, 5, 2]])6000((nlm[[1, 2, 5, 2]] − nlm1[[1, 2, 3, 2]])/nlm[[1, 2, 5, 2]])6000((nlm[[1, 2, 5, 2]] − nlm1[[1, 2, 3, 2]])/nlm[[1, 2, 5, 2]])6000 nlm[“ParameterErrors”]nlm[“ParameterErrors”]nlm[“ParameterErrors”]

nlm[“BestFitParameters”]nlm[“BestFitParameters”]nlm[“BestFitParameters”]

Length[shift]/17.Length[shift]/17.Length[shift]/17.

orderfit = Table[ListPlot[shift[[All, i]], PlotLabel → i, Frame → True,orderfit = Table[ListPlot[shift[[All, i]], PlotLabel → i, Frame → True,orderfit = Table[ListPlot[shift[[All, i]], PlotLabel → i, Frame → True,

FrameLabel → {“Suborder”, “velocity(km/s)”}, PlotRange → {−15, 15}], {i, 1, nos}]FrameLabel → {“Suborder”, “velocity(km/s)”}, PlotRange → {−15, 15}], {i, 1, nos}]FrameLabel → {“Suborder”, “velocity(km/s)”}, PlotRange → {−15, 15}], {i, 1, nos}] orderfiterror = Table[ListPlot[error[[All, i]], PlotLabel → i, Frame → True,orderfiterror = Table[ListPlot[error[[All, i]], PlotLabel → i, Frame → True,orderfiterror = Table[ListPlot[error[[All, i]], PlotLabel → i, Frame → True,