First Radial Velocity Results From the MINiature Exoplanet Radial Velocity Array (MINERVA)

The SBIG ST-i cameras originally used in the FAUs had two major problems. First, their small field of view (36 × 27) coupled with a surprisingly quick degradation of the telescope pointing meant that we could not blindly point to a target and be confident it would fall on the detector. We had to redo the pointing model weekly to ensure the pointing was sufficient for robust acquisition—a time-consuming task that must be done manually at night, when the telescopes would otherwise be carrying out science observations. The source of the pointing degradation is not clear, and the telescope manufacturers have not seen such degradation for other users, suggesting a problem with the robotic control software that we have not been able to fully investigate.

Second, the SBIGs had a high failure rate. During the initial month-long spectrograph commissioning, three out the four cameras in use experienced critical failures.

The manufacturers were aware that this problem affected a small batch of cameras and repaired them. After those repairs, the cameras performed better, but over the three years that followed, several more failures occurred. Given that replacing a failed camera requires a site visit, at a cost significantly greater than that of the cameras themselves, we decided to replace the SBIG cameras with ZWO ASI 174 cameras.

In the six months of daily use on four telescopes since their installation, we have not had a single camera failure. These cameras have a CMOS detector with 1936 × 1216, 5.86 μみゅーm pixels, and is similarly priced. This provides us an 85 × 54 field of view that allows us to robustly acquire our targets despite the pointing degradation of the telescopes. One downside to these cameras is that they are incompatible with our Black Box USB extenders, and so we had to move our control computers into the domes in order to use them.

The science fiber was originally purchased from CeramOptec and has a 50 μみゅーm octagonal core and 330 μみゅーm cladding. As is typical for CeramOptec, they actually have two slightly different claddings for octagonal fibers—one deposited onto the octagonal core to make it circular, and another with the core drilled out that they plug in with the circularized octagonal core. As was the case with our fiber, the boundary between the two claddings can be problematic if the indices are not well matched, since it can guide starlight through the cladding. Light transmitted through the cladding reduces the resolution of the spectrograph and the instrumental profile can vary dramatically as a function of how much light couples with the cladding—both of which are catastrophic for precision RV measurements. In addition, the cladding was too large to pack together at the focal plane with the required core-to-core spacing, and the standard Ferrule Connector (FC) connectors at each end, done by CeramOptec, suffered from severe Focal Ratio Degradation (FRD) and thus led to a major loss in throughput.

As a short-term solution, we re-terminated the original fibers to improve the FRD, and coupled the fibers into a short section of fiber with a 50 μみゅーm circular core and a 125 μみゅーm cladding to remove light from the cladding of the octagonal fiber. One side is butt-coupled to each of the four telescope fibers, and the other four ends combine into a V-groove at the spectrograph end, spaced 220 μみゅーm center to center.

As a long-term solution, we ordered a new fiber with a custom preform (a macro-sized piece of glass from which the fiber is drawn), which had a 50 μみゅーm octagonal core and 110 μみゅーm cladding from Polymicro, with the intention of packing seven of them cladding to cladding to allow for future expansion should we decide the significant (∼10%) crosstalk from such tightly packed fibers is manageable. Polymicro deposits the entire cladding onto the octagonal core to avoid the secondary cladding issue. However, our cladding was much thicker than typical, and during the lengthy deposition, the core melted and mixed into the cladding, creating the same effect as a secondary cladding. Light was still transmitted through the cladding of the fiber. A second attempt was no better, at which point they would not attempt the expensive process again. We could not afford a thinner cladding because we split the expensive custom preform with two other groups that required a thicker cladding.

So, our short term solution has become our final fiber solution. While the butt-couple is lossy, it makes the installation easier, it provides an easy point to add a double scrambler if we decide it is necessary in the future, and a change in the fiber geometry improves the near-field scrambling (Halverson et al. 2015). Meanwhile, MINERVA Australis (Addison et al. 2019) has created a fiber similar to our desired long-term solution, with enough for a spare if we decide it is worth the effort to replace in the future.

The KiwiSpec spectrograph was installed in 2015 December, and is a commercial adaptation of the spectrograph designed and described in Barnes et al. (2012) and Gibson et al. (2012), with a new camera designed by Prime Optics.³⁴ With a few exceptions highlighted below, it is as we described in Swift et al. (2015). Instead of the simultaneous etalon or thorium argon wavelength reference described in Swift et al. (2015), we used the thorium argon lamp during the installation as a rough wavelength solution and now rely solely on the iodine to provide the exact wavelength solution. The simultaneous wavelength reference in addition to iodine is unnecessary, and removing them allowed us to reduce the scattered light and increase the spacing between fibers, reducing the crosstalk to ≲0.1%.

We determined our total system throughput to be ∼5% using Doppler Tomography observations of KELT-24b without the iodine cell (Rodriguez et al. 2019). We computed the expected flux from the V = 8.33 host star between 6175 and 6185 Å and compared it to the actual flux in the extracted 1D spectrum at the same bandpass (at the peak of the blaze). This throughput estimate includes all losses, including the atmosphere, coatings, beam splitter, fiber coupling, Echelle, and charge-coupled device (CCD).

Light from each of the four telescopes is focused directly onto our 50 μみゅーm octagonal fibers at f/6.6 in our FAUえーゆー (see Section 3.1). Three meters before the spectrograph entrance, each of the fibers are butt-coupled to circular 50 μみゅーm fibers, which are then arranged into a V-groove at the entrance to the spectrograph, separated by 220 μみゅーm center to center. The light exits the four fibers and is collimated. A pupil mask truncates the beam to ensure the beam is precisely f/6, allowing for some focal ratio degradation within the fibers. The collimated light travels to the iodine stage, where the iodine gas cell can be moved into or out of the beam. The fibers are then re-imaged onto the entrance slit, and the light follows the path to the detector shown in Figure 1 of Barnes et al. (2012).

We empirically determined the resolving power per resolution element of the spectrograph by forward modeling high signal to noise spectra taken of the daytime sky and numerically solving for the FWHM of the fitted IP for each chunk. As expected, there is a slight wavelength dependence in our resolving power. The best-fit line of the resolving power per resolution element as a function of wavelength is R = 84,000 + (λらむだ − 5500 Å) × 10/Å, in good agreement with our theoretical expectation. We also plot the best-fit dispersion per resolution element as a function of wavelength for each chunk in Figure 1 for a representative night on a representative telescope.

Zoom In Zoom Out Reset image size

We have always had an exposure meter inside the spectrograph that picks off the reflection of the nearly collimated beam from the vacuum window. We have since added a V-band filter to approximate the bandpass of the spectrograph. The major downside to this design is that the exposure meter reports an average flux from all four telescopes, so we cannot use it to apply a per-telescope barycentric correction. Instead, we use the guide images from the FAUえーゆー (at ∼5 s cadence), with an aperture the size of the fiber drawn around the measured fiber position to determine the relative flux during the exposure for each telescope. We have confirmed that, when only one telescope is used, we can use the FAUえーゆー guide images reproduce the relative exposure meter flux to within the uncertainties, in order to compute a per-telescope barycentric correction.

The FAUえーゆー design described in Section 3.1 and Swift et al. (2015) flexes depending on its rotation angle and the telescope's altitude. This flexing causes the apparent position of the fiber on the acquisition camera to move by ∼10 μみゅーm over the sky—or 20% of the fiber diameter, which would be a significant source of light-loss if left uncorrected. We knew this would be a problem and the FAUえーゆー was designed to be able to locate the fiber tip on the acquisition unit by backlighting the fiber, but we had not yet fully fleshed out a solution at the time we wrote the Swift et al. (2015) paper. We considered using the exposure meter to refine the star's position, but that would dramatically increase our acquisition time since it would have to be done serially with each telescope and it can be difficult to make such a procedure robust during variable weather conditions. Ultimately, we added a disk of LEDs that swings in front of the V-groove to illuminate the fibers from inside the spectrograph. We do this before each exposure to refine the reference pixel to move the star to, and after to evaluate how much throughput might have been lost due to drift during the exposure. By evaluating a large number of these backlight images, we may be able to map the flexure and eliminate this step and/or compensate for drift during an exposure in the future.

Because the fibers do not provide much signal to noise in the wings of their profile, it is difficult to determine the pixel-to-pixel variations in the spectrograph detector with flat fields illuminated through the fiber. We added a light, mounted on the iodine stage, that illuminates a slit where the fibers are re-imaged. We use this flat field to correct for the pixel to pixel variations in the detector, as described in Section 6. The flat field lamp simply shines onto the entrance slit, with no attempt to match the f/6 science beam. We see no significant scattered light contamination with this approach, but we have yet to perform a detailed investigation.

While our iodine cell was designed to have counter-rotated wedges to eliminate fringing with minimal beam deflection, we believe the parallelism of the iodine cell faces was not within specification. As a result, when the iodine cell is in place, the position of the fibers shifts by almost the entire diameter of the fiber in the dispersion direction. Originally, the slit was only slightly oversized relative to the fiber size. We replaced it with a much wider slit to accommodate both the undeflected and deflected beams. While this significantly degrades the resolution of our slit flat fields, the flat only varies slowly as a function of color and its change across the degraded resolution is negligible. The resolution of our science images is set by the fiber size, not the slit width and therefore widening the slit has no impact on our science images.

Figures 2(a) and (b) show the cross section of the bluest and reddest orders, respectively, of the slit flat (without the iodine cell) overlayed on the same cross section of a daytime sky spectrum (with the iodine cell) to show that the slit flats give us adequate signal to calibrate the pixel to pixel variations under the science fibers, despite the deflection of the science fibers due to the iodine cell.

Zoom In Zoom Out Reset image size

In Swift et al. (2015), we described robotic control software based on the Robo-AO control software written in C (Riddle et al. 2012). In the following years, we determined that the growing popularity of Python, the many easily importable libraries, and vendor-provided APIs made it an attractive alternative to write and maintain the code while simultaneously allowing more complex features and capabilities. Our entire operational code base, written in Python, is hosted on Github.³⁵

A computer called "main" runs a 64 bit Ubuntu operating system and is responsible for most of the high-level operations. On startup, it begins three continuous functions. First, it operates an NTP server to which all other MINERVA clocks sync. It syncs itself to one of several stratum 2 time servers in the Western United States. Second, it runs a watchdog routine to monitor the temperature of the spectrograph in many locations and alert us via email if any are out of their operating range. Third, it runs a "domeControl" daemon that monitors the weather from several local weather stations: one at the MEarth building about 300 m away (Irwin et al. 2009), one at the HAT building about 230 m away (Bakos et al. 2002), one at the FLWO 1.2 m robotic telescope (home of KeplerCam) about 60 m away (Szentgyorgyi et al. 2005), and one we installed at the MINERVA building. It automatically evaluates if it is safe to open the domes based on cloud coverage, rain, humidity, wind speed, and Sun altitude, allowing overrides to open during cloudy weather or during the day for engineering. If it has been below freezing and wet (i.e., a potential for ice or snow), it sends us an email notifying us that manual approval is required to open the domes. Snow or ice on the roof can fall on the telescopes or overload the motors that open the shutters, which can prevent further robotic or remote control. Snow at the base of the enclosure can prevent it from fully opening. It must either melt or be cleared by the local site staff before we can safely open.

All automated safety checks must pass for 30 consecutive minutes and it must be requested to open before it will actually open. Once the enclosures are open, the criteria for closing are somewhat looser. These two requirements prevent rapid cycling of the enclosure during marginal conditions.

The domeControl daemon runs through its safety checks, sends a "heartbeat," and updates a status file every 15 s. The heartbeat is a firmware-level safety feature that protects us against a variety of potential failures. If a minute has elapsed and the enclosure has not received a heartbeat, it will automatically close, independent of any other activities. Should a failure of some kind prevent the enclosure from closing, emergency text messages are sent to several people to investigate immediately.

Finally, the "main" machine orchestrates the observations, which start each night at 4 pm local time via cron job, which we will describe in detail later in the next subsection.

Each of the four telescopes is controlled by its own computer running a 64 bit Windows 7 Professional operating system. Windows is required to run MaximDL for camera control and the PlaneWave Interface (PWI) software for telescope control. MaximDL controls our Andor and Apogee imagers and filter wheels for photometry (see Swift et al. 2015), as well as our ZWO imagers (see Section 3.1). We wrote our own server that runs locally on each Windows machine and can relay commands from our main control computer on the network to MaximDL. All images are saved to their own control computer on a drive that is cross mounted on the main computer. This allows us to run more complex image analysis like Source Extractor (Bertin & Arnouts 1996) and astrometry.net (Lang et al. 2010) to perform automated world coordinate solutions for acquisition and guiding or automatic exposure time adjustment during sky flats.

PWI hosts its own server that can be controlled by any computer on the network through simple XML commands. PlaneWave provided several example functions in Python, which we integrated into our software. Each windows computer has "scheduled tasks" (the Windows equivalent to a cron job) that reboot the computers daily and start the servers.

The spectrograph is controlled by two additional Windows 7 computers provided with the KiwiSpec spectrograph from KiwiStar Optics (a business unit of Callaghan Innovations). One computer is dedicated to the thermal control servo that maintains the spectrograph temperature and runs independently of all others.

The other computer operates much like the telescope control computers, and runs our server to relay commands from the main computer to the hardware connected directly to the spectrograph control computer. While KiwiStar Optics provided software to control the spectrograph manually, there was no API to interact with it robotically. We wrote our own spectrograph control software to enable robotic operations. This computer is responsible for operating the iodine stage, iodine heater, the spectrograph detector, the backlight, the flat field lamp, the exposure meter, and the vacuum pump and valves. The server also doubles as a watchdog that emails us if the vacuum pressure goes out of its operating range.

The observations begin at 4 pm local time. Our software computes the time it takes for a standard suite of biases, darks and flats for the spectrograph necessary for calibrating our RV observations. If photometric observations are desired, we upload a schedule file that contains the observations and corresponding calibrations. At 4 pm, the software computes how long the requested photometric calibrations will take, then begins the calibrations so they will finish 10 minutes before sunset.

Under normal spectroscopic observations, a dispatch scheduler reads active targets tabulated in a Google spreadsheet, and computes a score for each target based on the current time, the target's visibility, when it was last observed, how many times it has been observed that night, and how many times we would like it to be observed each night. In addition, it computes a weight for a single B-star observation that grows throughout the night until it is observed to ensure we obtain one B-star observation per night for calibrations. Further details about the MINERVA's scheduler can be found in S. A. Johnson et al. (2019, in preparation).

When photometric observations are requested for any subset of telescopes, we schedule our RV observations around the allocated times for the allocated telescopes. When only a subset of telescopes are scheduled for photometric observations, the others continue obtaining RV observations. Each telescope within the subset is capable of observing a distinct photometric target while the other telescopes obtain RVs on a single target.

During a typical spectroscopic observation, with all telescopes that are in RV mode, we slew to the target, turn on the backlight inside the spectrograph to illuminate the fibers, and expose the FAUえーゆー camera. This provides us with the precise pixel location of the fiber on the acquisition camera. Next, we do a fine acquisition to put the star onto the fiber. Because our targets are so bright, it is a safe assumption to move the brightest star in the field to the position of the fiber. Then we perform an autofocus and begin guiding to keep the star onto the fiber. While we use an Alt/Az telescope, we do not correct for field rotation, opting to keep the target star on the fiber and letting all other stars rotate about it to minimize the change in the pupil illumination during the exposure.

During a typical photometric observation, we can either cycle through a list of filters throughout some observing window (e.g., a predicted transit window), observe continuously in one filter, or take some number of exposures in each of several filters. While we have an off-axis guider for the imager, MaximDL does not allow us to control three cameras simultaneously, nor does it provide an API to switch between them robotically. However, the tracking performance of the CDK700s is superb and the direct drive motors have no periodic error, allowing us to take 5 minutes exposures unguided without any measurable trailing. Therefore, instead of the off-axis guider, we use the previous science image to correct for any long-term drift in tracking. This also has the advantage of not being subject to flexure or differential field rotation between the off-axis guider and the science camera, allowing us to easily maintain sub-pixel guiding accuracy throughout an hours-long transit—a capability that is critical to obtaining precise differential photometry.

We observe either RV targets or photometric targets throughout the night as desired, all the while monitoring the status of the dome and pausing if it closes. At the end of the night, we perform another set of calibrations, and close the dome. The data are backed up to our local RAID6 NAS (which can suffer two simultaneous drive failures without data loss), and our spectroscopic reduction pipeline is initiated at 10 am local time.

We obtain ${F}_{{I}_{2}}$ from a high resolution Fourier Transform Spectrometer (FTS) scan of the gaseous iodine cell as it is illuminated by a high signal-to-noise-ratio (S/N) continuum light source. We have two FTS scans of the MINERVA iodine cell. The first one was obtained at the Pacific Northwest National Laboratory a few years ago together with the CHIRON iodine cell (Tokovinin et al. 2013). The second FTS scan was done by Dr. Gillian Nave's group at NIST (e.g., see Nave 2011; Crause et al. 2018). Both scans were taken with the iodine cell at its operating temperature specification of 55 °C. Unfortunately, the two scans disagree in terms of line depths and line depth ratios, and we are further investigating the origin of this discrepancy. For concreteness, the results shown here use the second FTS scan, though both produce similar results.

The FTS scan is sufficient for determining our fiducial wavelength solution λらむだ(x) because the resolving power of the FTS (R ≈ 300,000) is about a factor of 4 greater than that of our KiwiSpec spectrograph (R ≈ 80,000). Because the molecular iodine lines span from 500 to 630 nm, the wavelength solution is determined solely within this range.

In our Doppler code, we choose between two functional forms for our model (IP(x)) of the spectrograph's instrumental profile. One form is a time-invariant spline function that is introduced in Section 7.5 and the other is a time-varying summation of satellite Gaussian profiles stacked on one central Gaussian profile that is described in Section 7.7. When using the former, we characterize our instrumental profile by using observations of a continuum light source with a high S/N while it illuminates the iodine gas cell in our KiwiSpec spectrograph. For reasons discussed later, we use the scattered sunlight of the daytime sky as our light source. For the Gaussian-like IP(x), however, we characterize the instrumental profile simultaneously with each stellar spectrum during our forward modeling. As a precaution, we also observe a B-type star each night (with iodine cell in place) to allow a more precise characterization of our instrumental profile as it changes over long periods of time.

A reference stellar spectrum is needed to determine the magnitude of the Doppler-shift seen when observing the science target. We use reference stellar spectra previously constructed by the CPS group using Keck/HIRES. They find the references by observing the science target without contamination from the iodine gas absorption lines. In this case, the IP(x) can be deconvolved with this observed (iodine-free) spectrum to get a reference stellar spectrum. In other words, F_⋆(λらむだ_ref(x)) is solved for via F_obs,ref(x) = F_⋆(λらむだ_ref(x))∗IP(x). This IP(x) and λらむだ_ref(x) here however are found using observations taken immediately before and after the iodine-free observation. The iodine-free observation of the science target is bracketed by iodine-calibrated observations of a nearby B-type star. The bracketed observations are particularly helpful if the instrumental profile is known to fluctuate on very short timescales, which is true in the case of Keck/HIRES. An IP(x) and λらむだ_ref(x) is evaluated for each iodine-calibrated B-type star observation and subsequently averaged. The resultant IP(x) is then deconvolved with the spectrum from the iodine-free science-target observation F_obs,ref(x) to get F_⋆(λらむだ_ref(x)), whose wavelength solution is assumed to be the averaged λらむだ_ref(x). The aforementioned CPS group refers to the reference stellar spectrum as the Deconvolved Stellar Spectral Template (DSST).

Using these DSSTs may limit our ability to accurately model F_obs(x) because the two spectrographs may suffer from different systematics. Furthermore, the observatory locations (Hawaii and Arizona in the U.S.) have significantly different water columns, dramatically changing the telluric features in our spectra. These differences may be a source of systematic error in our forward modeling procedure via the DSSTs. We will investigate the extent of these errors in the future. Meanwhile, we find that the DSSTs are sufficient for the first RV results of our RV survey. We have yet to derive our own reference stellar spectra because their development requires a substantial amount of observing time for each of our targets. The DSSTs are derived from observations with a higher S/N, and Keck's large aperture allows it to obtain such high S/N observations in a shorter time compared to MINERVA, which minimizes complications due to barycentric motion. In addition, unlike our fiber size, Keck's slit width is adjustable, allowing higher resolution templates which is helpful in developing the template.

The forward modeling procedure finds the best fit to F_obs(x) in Equation (3). Our Doppler code uses a least-squares algorithm to evaluate the best fit parameters. One of the parameters is the Doppler shift of the stellar spectrum F_⋆(λらむだ'(x)). The Doppler-shifted wavelength solution can be decomposed as λらむだ' = λらむだ_ref · (1 + z), where z is the Doppler shift from the radial velocity of the star and the motion of the telescope with respect to the star. This relative motion of the telescope is our barycentric velocity and it is dominated primarily by the Earth's rotation (∼0.5 km s⁻¹) and orbital motion around the Sun (∼30 km s⁻¹). The methods introduced in Wright & Eastman (2014) are used in our Doppler pipeline to correct for the telescope's barycentric velocity and subsequently determine the radial velocity of the star—which may or may not contain information about a planetary companion.

Our optical system can be divided into three general components: the telescope, the optical fibers, and the spectrograph. As the stellar rays trace this path, the optics distort the star's image. To determine the radial velocities from our observed spectra, we must know the manner of distortion that occurred en route to the spectrograph's CCD. The IP(x) is the shape a delta function would have when distorted by the entire optical system. The convention for determining an IP(x) is to assume it follows some function comprised of Gaussian structures that extend in the dispersion direction. This section describes how we have taken a unique approach.

In developing our new IP(x), we do not use a B-type star nor do we use calibration lamps. Instead, we take spectra of the sky during the daytime. Unlike our nightly stellar spectra, these daytime sky spectra yield a S/N/pixel > 150. With our stable spectrograph, the IP(x) characterized from daytime sky spectra should not change by the time we take stellar spectra at night. While the daytime sky is uniformly illuminated and may mask IP variations due to imperfect scrambling, such spectra give us a starting point to evaluate IP variations due to changes beyond the fiber.

The daytime sky spectra have proven to be a reliable source of data for determining the timescales at which our instrumental profile is stable. To find this timescale, we create an IP(x) that is time-invariant, which we refer to as the "fixed IP." By deducing the instrumental profile in these spectra, we can determine when and why it evolves.

We use the profile in the cross-dispersion direction to model the instrumental profile's shape in the dispersion direction. While the circular fiber makes this approximately correct, the distortions caused by the spectrograph's optical design certainly invalidate this assumption in detail. However, we can still evaluate the stability of the IP regardless, and we suspected that as long as the IP was stable and systematically wrong the same way each time, it would not impact the RV precision. Indeed, the results of long-term stability described later in this section justify this assumption.

The two-dimensional echellogram from MINERVA has four traces per order as shown in Figure 5. There is a trace for each telescope, and we divide each of the 18 orders into 15 "chunks." Each chunk consists of 128 columns of the trace. The total number of chunks in a frame is the number of chunks per order times the number of telescopes times the number of orders in the frame—1080. We define a distinct IP(x) for each of these chunks because the length of each chunk acts as a characteristic length scale for which the intrinsic instrumental profile changes. For this reason, we apply Equation (3) only over 128 pixels in the dispersion direction. Thus in practice, our forward modeling procedure is repeated for each chunk.

Zoom In Zoom Out Reset image size

To characterize the fixed IP, we first split a chunk (of two-dimensional daytime sky spectra) into 128 columns, or "crosscuts." As shown at the top of Figure 5, the traces are not perfectly horizontal. We therefore find the centroid of each crosscut and align the crosscuts' centroids so that the chunk is essentially as horizontal as the bottom of Figure 5. Ultimately, we want to normalize these crosscuts such that they collectively constrain the shape of our true instrumental profile. To align and normalize them, we start by assuming each crosscut can be modeled as a Gaussian,

$$\begin{eqnarray}&&G(x)=A\exp \left[-{\left(\displaystyle \frac{x-{x}_{0}}{\sigma }\right)}^{2}\right]+b,\end{eqnarray} \tag{ 4 }$$

where A is the amplitude, x is the pixel position in the cross-dispersion direction, x₀ is the centroid, σしぐま is the width, and b is the background of the raw spectra. We use a least-squares optimizer that follows the Levenberg–Marquardt algorithm to find the best fit parameters. We then subtract b from the crosscut data and integrate over this background-subtracted crosscut. After dividing the background-subtracted crosscut by this integral we can obtain the normalized background-subtracted crosscut. In other words, D_norm(x) = (D_raw(x) − b)/N, where D_norm is the normalized background-subtracted crosscut, D_raw is the original crosscut, and N is the normalization factor calculated by the aforementioned integration.

Once D_norm(x) is found for each of the 128 crosscuts, we fit a spline function of the third degree to all of them simultaneously. This cubic spline has breakpoints that are each separated by 6/10 of a pixel from each other. The spline acts as our preliminary fixed IP: IP_f(x) = spline(D_norm(x)). To find the optimal fixed IP, we perform an iterative process of modeling the D_raw(x) and subsequently evaluating a new D_norm(x) and IP_f(x).

Instead of a Gaussian, we use the previous spline fit model during the iterative process:

$$\begin{eqnarray}&&M(x)=N\times {\mathrm{IP}}_{f}(x-{\rm{\Delta }}x)+b,\end{eqnarray} \tag{ 5 }$$

where M(x) is the model for one crosscut and Δでるたx is a translational shift parameter. Now, the least-squares optimizer has only three parameters to evaluate: Δでるたx, N, and b. For the iterative process, we repeat the following procedure: find D_norm(x) of each crosscut, define the IP_f(x) for the chunk, optimize the 3 parameters of M(x) for each crosscut, calculate the reduced χかい² of all the crosscuts' data and M(x) models collectively, and lastly evaluate the difference between the previous iteration's reduced χかい² and the current iteration's reduced χかい². The most important distinction between iterations is the differing fixed IPs; when the χかい² gets lower, we conclude that the current iteration's IP_f(x) is better at modeling the spectrograph's instrumental profile than the previous iteration's IP_f(x). As the IP_f(x) gets better with each iteration, the difference in reduced χかい² values lessens. Once this difference is less than 10⁻⁴, any changes made to the IP_f(x) in the subsequent iterations are insignificant. The final iteration's IP_f(x) then becomes our nominal fixed IP.

Figure 6(a) is an example of the final fixed IP. The blue data points represent the D_norm(x) for all crosscuts of the final iteration. The orange line is the final IP_f(x). In this case, the data come from one chunk in one daytime sky exposure taken on 2018 June 19. The bottom plot illustrates the residuals, D_norm(x) − IP_f(x). The residuals are greatest near the center, where the shot noise is greatest, but they show no systematic structure. This suggests a good fit to the data.

Zoom In Zoom Out Reset image size

To test the longevity of our instrumental profile's stability for as long as possible, we construct the fixed IP with data taken at the time when we began this stability test and we used this fixed IP on spectra taken days, months, and a full year prior to the commencement of this test. We commenced this test after the end of our first full-season observing campaign (see Section 2). We then tested the fixed IP on spectra taken as far back in time as we saw fit for this test. We fit the same fixed IP to daytime sky spectra taken on 2017 October 2—about nine months away from the construction of the fixed IP. The result is presented in Figure 6(b) and it has the same general pattern of noise in its residuals as Figure 6(a). This implies that the instrumental profile has not changed within that time period. Note that for Figure 6(b) the 2018 fixed IP was used to model the 2017 spectra (via Equation (5)) and subsequently normalize its crosscuts.

To extend the timeline of this test further, we tried to use daytime sky spectra taken at the very beginning of that first observing campaign (2017 September 14). Unfortunately, our daytime sky spectra taken within those first two weeks, between 2017 September 14 and October 2, were of poor quality and had a relatively low S/N until we resolved the issue. Therefore, we tried to use daytime sky spectra taken before the 2017 monsoon season and thus before our first full-season observing campaign. Fortunately for this test, we took many daytime sky spectra back in 2017 March and April. We therefore extend the timeline of this test to 2017 April.

Figure 6(c) shows how our 2018 fixed IP is used to model data taken on 2017 April 3. The residuals here show strong systematic structure when compared to Figures 6(a) and (b). From these three examples, it is clear that our spectrograph's instrumental profile was stable from 2017 October 2 to 2018 June 19 but not from 2017 April 3 to October 2.

The spectrograph's instrumental profile evolved significantly within the window of six months between 2017 April 3 and October 2. As explained in Section 4 and shown in Figure 3(a), during this time, the pressure rose dramatically for an extended period of time during the 2017 monsoon season after a brief power outage. This event permanently altered the spectrograph in such a way that the environment could not naturally return to its original instrumental profile when the pressure returned to its original operating specification. This means that the instrumental profile might have been stable for longer than nine months if the power outage and subsequent pressure instability did not occur. When characterizing the instrumental profile with a fixed IP, a new fixed IP must be used whenever an event such as this occurs. If this is not done, a situation like that of Figure 6(c) is likely to occur.

We explored the possibility of modeling a fixed component of the IP with a time-varying component to account for changes in the IP due to perturbations like that which is seen in Figure 6(c) and explained in Figure 3. We can formulate such an IP as

$$\begin{eqnarray}\begin{array}{rcl}{\mathrm{IP}}_{{GH}}(x) & = & \displaystyle \sum _{n}{A}_{n}{\left(\displaystyle \frac{2}{\pi {\sigma }^{2}}\right)}^{1/4}\displaystyle \frac{1}{\sqrt{n!{2}^{n}}}{H}_{n}\left(\displaystyle \frac{x\sqrt{2}}{\sigma }\right){e}^{-{\left(\tfrac{x}{\sigma }\right)}^{2}}\\ & & +\ \ {\mathrm{IP}}_{f}(x).\end{array}\end{eqnarray} \tag{ 6 }$$

This function includes the sum of the products between Gaussians of amplitude A_n and Hermite polynomials H_n. The systematic structure seen in Figure 6(c) suggests that we may be able to model the time variable component of the IP with fewer free parameters than a purely time-variable IP, and thus preserve the signal in the spectrum to constrain the Doppler signal we care about rather than the instrumental profile. Unfortunately, the number of additional parameters required to accurately model the time-variable component was comparable to purely time-variable IP described in the following section, and therefore offers no advantage. Additionally, the GH parameterization of the IP is not as well behaved as the sum of Gaussians in our forward modeling procedure.

Our sum-of-Gaussians IP,

$$\begin{eqnarray}&&{\mathrm{IP}}_{G}(x)=\displaystyle \sum _{n}{A}_{n}\exp \left[-{\left(\displaystyle \frac{x-{x}_{0,n}}{{\sigma }_{n}}\right)}^{2}\right],\end{eqnarray} \tag{ 7 }$$

is a time-dependent IP that is described in detail by Valenti et al. (1995). Notice here that we do not include the fixed IP. Also note that A₀ is fixed to 1 so that there is one large central Gaussian while all other Gaussian components act as small satellite Gaussians. To test this IP_G against the fixed IP, we calculated RVs for an RV standard star and a planet-hosting star after forward modeling the data with each of the IPs, as described later in Section 8.

With one telescope, we acquired the radial velocities of HD 122064 during the months of 2016 May and June. As a test, we used both the fixed IP and the sum-of-Gaussians IP in our forward modeling to generate two distinct RV data sets which derive from the same spectra.

For the purposes of the instrumental profile's stability test, we only used one daytime sky exposure to construct the fixed IP. Whenever we perform our forward modeling procedure with the fixed IP, we make the fixed IP more robust by using multiple daytime sky exposures. We perform the same procedure as described in Section 7.5, except the number of crosscuts that we simultaneously fit a cubic spline function to is equal to 128 times the number of daytime sky exposures we use. Each set of 128 crosscuts comes from the same chunk of distinct daytime sky exposures. For the May/June data set, we use ∼5 daytime sky exposures to construct a fixed IP for each chunk and these exposures are somewhat evenly distributed throughout the 1.5 months timescale of the data set. This fixed IP is used to generate the RVs in Figure 7(a) while the sum-of-Gaussians IP is used to produce the RVs of Figure 7(b).

Zoom In Zoom Out Reset image size

After the RVs are extracted, we compute the Allan variance to determine the level of precision MINERVA can achieve. In Figure 8, the line represents the precision if the binned data only contained white noise. We use an error-weighted, overlapping Allan variance to determine the limit for which we can bin down the given data set before it is dominated by systematic errors (Allan 1966; Malkin 2011).

Zoom In Zoom Out Reset image size

We have seventy-five radial velocity measurements tabulated in Table 2 and shown in Figures 7(a) and (b). Figures 8(a) and (b) suggest that a bin size of six roughly marks the limit for which the respective binned data sets begin to deviate from white noise. At this binning, we are sensitive to variations below the 2 m s⁻¹ level for our measurements of this RV standard star. The precision achieved is 1.8 m s⁻¹ for the fixed IP approach and 1.9 m s⁻¹ for the sum-of-Gaussians approach. This could potentially change depending on the standard star or the amount of data we have for a given star. To confirm this, we plan on performing the same test for observations of other RV standard stars. Regardless however, these RVs evaluated through our Doppler pipeline suggest that our fixed IP is doing just as well as our sum-of-Gaussians IP.

MINERVA observations of 51 Peg were taken with one telescope in 2017 October. Again, we use the fixed IP and sum-of-Gaussians IP to extract the radial velocities. It is wise to see if our radial velocities can confirm the existence and characteristics of exoplanet systems. 51 Peg b is the first of such exoplanets to be tested. We use EXOFASTv2 to constrain the properties of this exoplanet system. Our stellar parameters are informed by the broadband photometry summarized in Table 3. The RVs and resultant 51 Peg b properties derived with both IPs are so similar that we only show the results produced with the fixed IP. The unbinned RVs and the EXOFASTv2-generated orbital solution are illustrated in Figure 9(a), tabulated in Table 4, and summarized in Table 5. Figure 9(b) shows the same but with the time series folded to the phase of the planet's orbital period.

Zoom In Zoom Out Reset image size

Table 5.  Properties of 51 Peg

Parameter	Description	EXOFASTv2	Reference (Bu06)	Δでるたσしぐま
Stellar Parameters:
M*	Mass (M⊙)	${1.095}_{-0.07}^{+0.066}$	1.09 ± 0.109	0.039
R*	Radius (R⊙)	${1.139}_{-0.03}^{+0.03}$	1.03 ± 0.07	1.431
L*	Luminosity (L⊙)	${1.39}_{-0.053}^{+0.056}$	⋯	⋯
ρろー*	Density (cgs)	${1.05}_{-0.11}^{+0.12}$	⋯	⋯
log g	Surface Gravity (cm s−2)	${4.365}_{-0.041}^{+0.036}$	4.449 ± 0.06	1.2
Teff	Effective Temperature (K)	${5875}_{-76}^{+73}$	5787 ± 44	1.002
[Fe/H]	Metallicity	${0.201}_{-0.069}^{+0.07}$	0.2 ± 0.03	0.013
[Fe/H]0	Initial Metallicity	${0.204}_{-0.065}^{+0.064}$	⋯	⋯
Age	Age (Gyr)	${4.0}_{-2.5}^{+3.4}$	⋯	⋯
EEP	Equal Evolutionary Point	${367}_{-37}^{+39}$	⋯	⋯
AV	V-band extinction	${0.051}_{-0.036}^{+0.045}$	⋯	⋯
σしぐまSED	SED photometry error scaling	${3.01}_{-0.7}^{+1.1}$	⋯	⋯
d	Distance (pc)	${15.468}_{-0.029}^{+0.03}$	15.36 ± 0.18	0.592
πぱい	Parallax (mas)	${64.65}_{-0.12}^{+0.12}$	⋯	⋯
Planetary Parameters:
P	Period (days)	${4.236}_{-0.025}^{+0.028}$	4.230785 ± 0.000036	0.209
RP	Radius (RJ)	${1.242}_{-0.094}^{+0.093}$	⋯	⋯
TC	Time of Conjunction (BJDTDB)	${2458048.885}_{-0.062}^{+0.081}$	⋯	⋯
T0	Optimal Time of Conjunction (BJDTDB)	${2458044.653}_{-0.068}^{+0.074}$	⋯	⋯
a	Semimajor Axis (au)	${0.0528}_{-0.0012}^{+0.0011}$	0.0527 ± 0.003	0.031
i	Inclination (°)	${61}_{-27}^{+20}$	⋯	⋯
e	Eccentricity	${0.051}_{-0.036}^{+0.062}$	0.013 ± 0.012a	1.001
ωおめが*	Argument of Periastron (°)	${73}_{-79}^{+200}$	58a	⋯
Teq	Equilibrium Temperature (K)	${1315}_{-17}^{+18}$	⋯	⋯
MP	Mass (MJ)	${0.56}_{-0.075}^{+0.3}$	⋯	⋯
K	RV semi-amplitude (m s−1)	${57.6}_{-3.7}^{+3.6}$	55.94 ± 0.69	0.441
log K	Log of RV semi-amplitude	${1.76}_{-0.029}^{+0.026}$	⋯	⋯
RP/R*	Radius of planet in stellar radii	${0.112}_{-0.0089}^{+0.0091}$	⋯	⋯
a/R*	Semimajor axis in stellar radii	${9.98}_{-0.38}^{+0.35}$	⋯	⋯
ρろーP	Density (cgs)	${0.38}_{-0.1}^{+0.25}$	⋯	⋯
log gP	Surface Gravity	${2.97}_{-0.11}^{+0.2}$	⋯	⋯
Θしーた	Safronov Number	${0.0439}_{-0.007}^{+0.025}$	⋯	⋯
$\langle F\rangle$	Incident Flux (109 erg s−1 cm2)	${0.675}_{-0.035}^{+0.038}$	⋯	⋯
TP	Time of Periastron (BJDTDB)	${2458050.24}_{-1.0}^{+0.98}$	⋯	⋯
TS	Time of Eclipse (BJDTDB)	${2458050.941}_{-0.087}^{+0.073}$	⋯	⋯
TA	Time of Ascending Node (BJDTDB)	${2458047.78}_{-0.079}^{+0.065}$	⋯	⋯
TD	Time of Descending Node (BJDTDB)	${2458049.93}_{-0.063}^{+0.07}$	⋯	⋯
e cos ωおめが*		${0.018}_{-0.063}^{+0.032}$	⋯	⋯
e sin ωおめが*		${0.008}_{-0.052}^{+0.032}$	⋯	⋯
MP sin i	Minimum mass (MJ)	${0.484}_{-0.037}^{+0.036}$	0.472 ± 0.039	0.223
MP/M*	Mass ratio	${0.000487}_{-0.000063}^{+0.00028}$	⋯	⋯
d/R*	Separation at mid transit	${10.07}_{-0.55}^{+0.59}$	⋯	⋯
PT	A priori non-gazing transit prob	${0.0881}_{-0.0051}^{+0.0054}$	⋯	⋯
PT,G	A priori transit prob	${0.1104}_{-0.0061}^{+0.0063}$	⋯	⋯
PS	A priori non-gazing eclipse prob	${0.0905}_{-0.005}^{+0.0065}$	⋯	⋯
PS,G	A priori eclipse prob	${0.1133}_{-0.0057}^{+0.0076}$	⋯	⋯
rms	Root Mean Square of residuals (m s−1)	5.6	7.0	⋯
Telescope Parameters:
γがんまrel	Relative RV Offset (m s−1)	$-{5.0}_{-2.5}^{+2.5}$	⋯	⋯
σしぐまJ	RV Jitter (m s−1)	${7.7}_{-1.9}^{+2.8}$ b	3.7b	⋯
${\sigma }_{J}^{2}$	RV Jitter Variance	${59}_{-26}^{+51}$	⋯	⋯
Nobs	Number of observations	15	256	⋯

Notes.

^aThe uncertainties reported by Bu06 for e and ωおめが_* are non-Gaussian because the uncertainty of e is comparable to e, i.e., σしぐま_e ≳ e/2. ^bOur RV jitter is informed by the aforementioned radial velocities while the Bu06 value is informed merely by other stellar parameters.

Download table as:  ASCIITypeset image

The median and 68% confidence intervals determined using EXOFASTv2 with the MINERVA data for all parameters of the 51 Peg system are listed in Table 5. We only employ constraints on three of the stellar parameters. We impose a prior on the stellar metallicity ([Fe/H]) of 0.20 ± 0.07 dex from spectroscopy described in Fuhrmann et al. (1997). We set the V-band extinction's (A_V) upper limit to 0.11811 magnitudes, using the dust maps from Schlafly & Finkbeiner (2011). Lastly, we impose a Gaussian prior on the parallax (πぱい) of 64.65 ± 0.12 mas from Gaia Collaboration et al. (2016, 2018). These priors, coupled with a MIST stellar evolution model (Choi et al. 2016; Dotter 2016) and an SED model, constrain the properties of the host star.

We compare our results with the values in Butler et al. (2006) (hereafter referred to as Bu06). They cite the SPOCS catalog (Valenti & Fischer 2005) as the source for most of their stellar parameters. The Bu06 distance d to the star and its uncertainty are from the Hipparcos³⁷ catalog. Valenti & Fischer (2005) suggest that the typical uncertainties for their stellar parameters among their catalog of stars are 0.06 dex for log g, 44 K for T_eff, and 0.03 dex for [Fe/H]. Bu06 assumes a 10% uncertainty for the stellar mass M_*. Bu06 does not report a value for stellar radius R_*. We therefore calculate this using their stated log g and M_* values. The uncertainty in the stellar radius is found using propagation of error between those two parameters.

The reference planetary and telescope parameters shown in Table 5 were derived solely from Bu06. Their observations were taken at Lick Observatory using the Hamilton spectrograph (Vogt 1987), the 3.9 m Anglo-Australian Telescope using UCLES (Diego et al. 1990), and the Keck Observatory using HIRES (Vogt et al. 1994). Their quoted RV jitter σしぐま_J, however, does not come from their observations. Their jitter comes from the model developed by Wright (2005), which was informed by a sample of 531 stars observed at Keck that had known activity levels, colors, and parallaxes. In general, the jitter depends on the spectral type of the star and the instrument observing it. The model by Wright (2005) uses a stellar activity indicator, B − V color, and difference in magnitude above the main sequence to approximate the stellar jitter.

Values not reported by Bu06 are marked with ellipses (...) in Table 5. For the parameters that have a value and errors reported by Bu06, we state the discrepancy between our values and theirs in terms of 1σしぐま uncertainty. We define this discrepancy as

$$\begin{eqnarray}&&{\rm{\Delta }}\sigma =\displaystyle \frac{{N}_{1}-{N}_{2}}{\sqrt{{\left({\sigma }_{1,L}\right)}^{2}+{\left({\sigma }_{2,U}\right)}^{2}}},\mathrm{when}\quad {N}_{1}\gt {N}_{2},\end{eqnarray} \tag{ 9 }$$

in which σしぐま_1,L is the lower error bar for N₁ and σしぐま_2,U is the upper error bar for N₂. Ideally, the discrepancy should be less than 1σしぐま. Seeing as the discrepancies evaluated for the tabulated parameters—most importantly the planetary parameters—are ≲1.0σしぐま, we find good agreement with results quoted in literature.

We also analyze the residuals of the fit to 51 Peg b to determine our precision. An Allan variance plot of the residuals is shown in Figure 10, demonstrating our per point scatter (5.6 m s⁻¹) and systematic floor (4.2 m s⁻¹) when two observations are binned is considerably worse than achieved for HD 122064, most likely due to stellar jitter of the more active host star. We note that fitting for the period and eccentricity may absorb some of the excess scatter and make our precision appear better than it is. To address this, we also ran a fit fixing the period to the Bu06 value (4.230785 days) and fixing the eccentricity to zero, but found an insignificant increase in the per-point residual rms (5.6 to 6.0 m s⁻¹) and an insignificant decrease in the binned-by-2 residual rms (4.3 to 4.2 m s⁻¹). Therefore, we conclude that fitting the period and eccentricity does not have a major impact on our inferred precision.

Zoom In Zoom Out Reset image size