First radial velocity results from the MINiature Exoplanet Radial Velocity Array (MINERVA)

TL;DR

MINERVA, a dedicated array of four robotic telescopes, achieves high-precision radial velocity measurements, demonstrating instrument stability and accuracy in detecting exoplanets like super-Earths in nearby star systems.

Contribution

This paper details the improvements to MINERVA's hardware and software, and demonstrates stable, precise radial velocity measurements over nine months using a time-independent instrumental profile.

Findings

Instrumental profile stable for at least nine months

Achieved 1.8 m s$^{-1}$ radial velocity precision

Results consistent with literature for 51 Peg

Abstract

The MINiature Exoplanet Radial Velocity Array (MINERVA) is a dedicated observatory of four 0.7m robotic telescopes fiber-fed to a KiwiSpec spectrograph. The MINERVA mission is to discover super-Earths in the habitable zones of nearby stars. This can be accomplished with MINERVA's unique combination of high precision and high cadence over long time periods. In this work, we detail changes to the MINERVA facility that have occurred since our previous paper. We then describe MINERVA's robotic control software, the process by which we perform 1D spectral extraction, and our forward modeling Doppler pipeline. In the process of improving our forward modeling procedure, we found that our spectrograph's intrinsic instrumental profile is stable for at least nine months. Because of that, we characterized our instrumental profile with a time-independent, cubic spline function based on the profile in the cross dispersion direction, with which we achieved a radial velocity precision similar to using a conventional "sum-of-Gaussians" instrumental profile: 1.8 m s$^{-1}$ over 1.5 months on the RV standard star HD 122064. Therefore, we conclude that the instrumental profile need not be perfectly accurate as long as it is stable. In addition, we observed 51 Peg and our results are consistent with the literature, confirming our spectrograph and Doppler pipeline are producing accurate and precise radial velocities.

Tables

Table 1
Stars of Interest
HD	R.A. (J2000)	Decl. (J2000)	$V$	SpType
122064	13 57 32.059	+61 29 34.30	6.52	K4V
217014	22 57 27.980	+20 46 07.80	5.46	G5V

Table 2
HD 122064 RVs and RV Errors (m s-1)
Date - 2,457,500	IPf		IPG
(BJDTDB)	RV	Error	RV	Error
24.672037	-3.95	2.09	0.51	1.92
24.693414	2.47	2.03	7.07	1.99
24.714667	4.41	2.12	9.56	1.90
28.674272	7.66	1.97	12.03	1.70
…	…	…	…	…
NOTE.—RVs for HD 122064 displayed in Figure 7. The mean error is 2.1 m/s for each data set: the data derived from the fixed IP and data from the sum-of-Gaussians IP. A portion is shown here for guidance regarding its form and content.

Table 3
51 Peg Magnitudes
Band	Mag.	Used	Catalog’s
		Mag. Error	Mag. Error
Tycho-2 Catalog (Høg et al., 2000)
BT	6.249	0.020	0.014
VT	5.526	0.020	0.009
2MASS Catalog (Cutri et al., 2003)
J2M	4.655	0.300	0.300
H2M	4.234	0.270	0.270
K2M	3.911	0.020	0.020
WISE Catalog (Wright et al., 2010)
WISE1	3.909	0.387	0.387
WISE2	3.624	0.246	0.246
WISE3	3.929	0.030	0.016
WISE4	3.904	0.100	0.024

Table 4
RV Results for 51 Peg
Date - 2,458,000	RV	Error	Residuals
(BJDTDB)	(m s-1)	(m s-1)	(m s-1)
38.897636	42.28	3.11	0.41
40.748225	-29.54	2.60	5.05
40.883873	-49.87	3.19	-5.30
43.703150	56.90	2.81	9.99
43.888337	44.76	3.13	3.24
44.892530	-22.88	2.99	3.94
45.888237	-68.62	3.27	-3.70
47.681088	38.81	2.72	-9.92
48.686860	7.81	2.70	-1.24
48.844736	-10.45	3.04	-7.23
50.718624	-24.11	2.30	2.54
52.787278	25.94	3.05	6.67
53.797358	-47.75	2.59	8.28
55.696227	37.76	2.81	4.24
63.696364	-1.04	3.99	2.76
NOTE.—The mean formal error derived from the third column is 3.0 m s-1. The fourth column represents the residuals between the RV data (of the second column) and the best-fit orbital solution from EXOFASTv2. The RMS of the residuals is 5.6 m s-1, and the systematic error floor, achieved when binning by 2, is 4.2 m s-1.

Table 5
Properties of 51 Peg
Parameter	Description	EXOFASTv2	Reference (Bu06)	$\Delta \sigma$
Stellar Parameters:
$M_{\ast }$	Mass ($M_{\odot }$)	${1.095}_{-0.07}^{+0.066}$	$1.09\pm 0.109$	0.039
$R_{\ast }$	Radius ($R_{\odot }$)	${1.139}_{-0.03}^{+0.03}$	$1.03\pm 0.07$	1.431
$L_{\ast }$	Luminosity ($L_{\odot }$)	${1.39}_{-0.053}^{+0.056}$
${\rho }_{\ast }$	Density (cgs)	${1.05}_{-0.11}^{+0.12}$
$\log g$	Surface Gravity (cm s-2)	${4.365}_{-0.041}^{+0.036}$	$4.449\pm 0.06$	1.2
$T_{\mathrm{eff}}$	Effective Temperature (K)	${5875}_{-76}^{+73}$	$5787\pm 44$	1.002
$\left[\mathrm{Fe}/H\right]$	Metallicity	${0.201}_{-0.069}^{+0.07}$	$0.2\pm 0.03$	0.013
${\left[\mathrm{Fe}/H\right]}_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}$
$A_V$	V-band extinction	${0.051}_{-0.036}^{+0.045}$
${\sigma }_{SED}$	SED photometry error scaling	${3.01}_{-0.7}^{+1.1}$
$d$	Distance (pc)	${15.468}_{-0.029}^{+0.03}$	$15.36\pm 0.18$	0.592
$\pi$	Parallax (mas)	${64.65}_{-0.12}^{+0.12}$
Planetary Parameters:
$P$	Period (days)	${4.236}_{-0.025}^{+0.028}$	$4.230785\pm 0.000036$	0.209
$R_P$	Radius ($R_J$)	${1.242}_{-0.094}^{+0.093}$
$T_C$	Time of Conjunction (BJDTDB)	${2458048.885}_{-0.062}^{+0.081}$
$T_0$	Optimal Time of Conjunction (BJDTDB)	${2458044.653}_{-0.068}^{+0.074}$
$a$	Semi-major Axis (AU)	${0.0528}_{-0.0012}^{+0.0011}$	$0.0527\pm 0.003$	0.031
$i$	Inclination (Degrees)	${61}_{-27}^{+20}$
$e$	Eccentricity	${0.051}_{-0.036}^{+0.062}$	$0.013\pm 0.012$a	1.001
${\omega }_{\ast }$	Argument of Periastron (Degrees)	$-{73}_{-79}^{+200}$	58a
$T_{eq}$	Equilibrium Temperature (K)	${1315}_{-17}^{+18}$
$M_P$	Mass ($M_J$)	${0.56}_{-0.075}^{+0.3}$
$K$	RV semi-amplitude (m/s)	${57.6}_{-3.7}^{+3.6}$	$55.94\pm 0.69$	0.441
$\log K$	Log of RV semi-amplitude	${1.76}_{-0.029}^{+0.026}$
$R_P/R_{\ast }$	Radius of planet in stellar radii	${0.112}_{-0.0089}^{+0.0091}$
$a/R_{\ast }$	Semi-major axis in stellar radii	${9.98}_{-0.38}^{+0.35}$
${\rho }_P$	Density (cgs)	${0.38}_{-0.1}^{+0.25}$
$\log g_P$	Surface Gravity	${2.97}_{-0.11}^{+0.2}$
$\Theta$	Safronov Number	${0.0439}_{-0.007}^{+0.025}$
$⟨F⟩$	Incident Flux (109 erg s-1 cm2)	${0.675}_{-0.035}^{+0.038}$
$T_P$	Time of Periastron (BJDTDB)	${2458050.24}_{-1.0}^{+0.98}$
$T_S$	Time of Eclipse (BJDTDB)	${2458050.941}_{-0.087}^{+0.073}$
$T_A$	Time of Ascending Node (BJDTDB)	${2458047.78}_{-0.079}^{+0.065}$
$T_D$	Time of Descending Node (BJDTDB)	${2458049.93}_{-0.063}^{+0.07}$
$e\cos {\omega }_{\ast }$		$-{0.018}_{-0.063}^{+0.032}$
$e\sin {\omega }_{\ast }$		$-{0.008}_{-0.052}^{+0.032}$
$M_P\sin i$	Minimum mass ($M_J$)	${0.484}_{-0.037}^{+0.036}$	$0.472\pm 0.039$	0.223
$M_P/M_{\ast }$	Mass ratio	${0.000487}_{-0.000063}^{+0.00028}$
$d/R_{\ast }$	Separation at mid transit	${10.07}_{-0.55}^{+0.59}$
$P_T$	A priori non-gazing transit prob	${0.0881}_{-0.0051}^{+0.0054}$
$P_{T,G}$	A priori transit prob	${0.1104}_{-0.0061}^{+0.0063}$
$P_S$	A priori non-gazing eclipse prob	${0.0905}_{-0.005}^{+0.0065}$
$P_{S,G}$	A priori eclipse prob	${0.1133}_{-0.0057}^{+0.0076}$
RMS	Root Mean Square of residuals (m/s)	5.6	7.0
Telescope Parameters:
${\gamma }_{\mathrm{rel}}$	Relative RV Offset (m/s)	$-{5.0}_{-2.5}^{+2.5}$
${\sigma }_J$	RV Jitter (m/s)	${7.7}_{-1.9}^{+2.8}$b	3.7b
${\sigma }_J^2$	RV Jitter Variance	${59}_{-26}^{+51}$
$N_{\mathrm{obs}}$	Number of observations	15	256
NOTE.
aThe uncertainties reported by Bu06 for $e$ and ${\omega }_{\ast }$ are non-Gaussian because the uncertainty of $e$ is comparable to $e$, i.e. ${\sigma }_e\gtrsim e/2$.
bOur RV jitter is informed by the aforementioned radial velocities while the Bu06 value is informed merely by other stellar parameters.

Equations

$$F(y)=\mathbf{p}(x,y)F(x)+n(y).$$

$$\mathbf{p}(x,y)=N(x)e^{\left(-0.5{\left(\frac{|y-y_c(x)|}{\sigma (x)}\right)}^{p(x)}\right)}.$$

$$F_{obs}(x)=[F_{I_2}(\lambda (x))\times F_{\star }({\lambda }^{{}^{\prime }}(x))]\ast \text{IP}(x),$$

$$G(x)=A\textrm{ }\textrm{exp}\bigg{[}-\bigg{(}\frac{x-x_{0}}{\sigma}\bigg{)}^{2}\bigg{]}+b,$$

$$M(x)=N\times {\text{IP}}_f(x-\Delta x)+b,$$

$$\begin{split}\text{IP}_{GH}(x)=&\sum_{n}A_{n}\bigg{(}\frac{2}{\pi\sigma^{2}}\bigg{)}^{1/4}\frac{1}{\sqrt{n!2^{n}}}H_{n}\bigg{(}\frac{x\sqrt{2}}{\sigma}\bigg{)}e^{-\big{(}\frac{x}{\sigma}\big{)}^{2}}\\ &+\text{ IP}_{f}(x).\end{split}$$

$$\textrm{IP}_{G}(x)=\sum_{n}A_{n}\textrm{exp}\bigg{[}-\bigg{(}\frac{x-x_{0,n}}{\sigma_{n}}\bigg{)}^{2}\bigg{]},$$

$$M_{p}\sin i\approx K\sqrt{1-e^{2}}\bigg{(}\frac{PM^{2}_{*}}{2\pi G}\bigg{)}^{1/3}$$

$$\Delta \sigma =\frac{N_1-N_2}{\sqrt{({\sigma }_{1,L}{)}^2+({\sigma }_{2,U}{)}^2}},\text{when }{N}_1>N_2,$$

Full text

instrumentation: spectrographs — methods: data analysis — methods: observational — planets and satellites: detection — planets and satellites: general — techniques: radial velocities — techniques: spectroscopic

††journal: PASP

1 Introduction

The discovery of the first planets orbiting solar-type stars was achieved using Doppler spectroscopy (Campbell et al., 1988; Latham et al., 1989; Mayor & Queloz, 1995). As the first exoplanet detections and confirmations were made, Doppler spectroscopy instruments gradually improved from attaining a radial velocity (RV) precision of $\sim$15 m s*-1* (Campbell et al., 1988) to $\sim$3 m s*-1* (Butler et al., 1996) thanks to the advent of the iodine absorption cell technique. Two decades later, the next generation of precision RV instruments aims for instrumental stability at the 30 cm s*-1* level (Wright & Robertson, 2017). However, our sensitivity to exoplanets is likely limited by stellar activity at the $\sim$1 m s*-1* level for most stars (e.g., Saar & Donahue, 1997; Haywood et al., 2016). Detections below this level will not be achieved until astrophysical noise sources are understood as well as sources of instrumental noise. Observing with high cadence throughout a planet’s full orbit may allow us to understand and correct for non-planetary RV signals induced by stellar activity (O’Toole et al., 2008; Pepe et al., 2011; Dumusque, 2012).

The MINiature Exoplanet Radial Velocity Array (MINERVA) is a dedicated observatory aiming for both high cadence and high precision RV measurements (Swift et al., 2015). It is a robotic array of four 0.7 m telescopes located on Mt. Hopkins in Arizona. The MINERVA mission ultimately has two objectives.

The primary science objective is to detect and characterize super-Earths in the habitable zones of nearby stars. Our RV target list is a subset of the targets monitored during the NASA/UC $\eta_{\oplus}$ Survey performed by the California Planet Search (CPS) group at the Keck Observatory using the HIRES spectrograph (Howard et al., 2009). Out of the 230 GKM stars they surveyed, 166 are considered chromospherically quiet (Wright et al., 2004; Isaacson & Fischer, 2010). The MINERVA RV target list consists of 125 of the brightest ($V\lesssim 8$) chromospherically quiet stars from their survey that can be observed from southern Arizona. With MINERVA’s effective aperture of 1.4 m and use of the NASA/UC $\eta_{\oplus}$ targets, the RV precision goal of the MINERVA mission was set to detect planets at the 80 cm s*-1* level (Swift et al., 2015). At this level, we plan to characterize super-Earths while providing insight into the importance of cadence as a tool for understanding the problem of stellar activity. We show that we are about a factor of 2 of that goal in §8.1, which is already within the top tier of the current generation of precision RV instruments. Coupled with our unmatched observational cadence, we are already operating in a unique parameter space that will enable us to detect new planets and provide valuable insight about the importance of cadence in understanding stellar jitter. We can do this with our cost-efficient, four-telescope, robotic array observing at an unprecedented cadence. The high cadence is attributed to the autonomous, flexible target scheduling, and quick slewing of the CDK-700 telescopes. Most importantly, the majority of the robotic array’s time is not split between multiple teams or science goals.

The secondary science objective is to search for transits of the super-Earths we find. This requires a broadband photometric precision of $<$1 mmag in the optical: a goal that has already been demonstrated by Swift et al. (2015). Multiband light curves provide information that otherwise cannot be deduced from Doppler spectroscopy alone. For example, the minimum mass of the planet can be found from radial velocities, but if the planet happens to transit, the transit photometry can determine the radius and inclination of the planet (see, e.g. Winn, 2010). Therefore, both exoplanet detection methods used together can indicate the true planetary mass and bulk density.

MINERVA’s secondary objective has already contributed to a variety of exoplanet science endeavours (Swift et al., 2015; Vanderburg et al., 2015; Croll et al., 2017; Lund et al., 2017; Pepper et al., 2017; Rodriguez et al., 2017; Siverd et al., 2018; Labadie-Bartz et al., 2019). Thus, in this work we focus on the commencement of MINERVA’s primary objective. We report our survey performance in §2, the changes to our hardware since our last paper in §3, the environmental stability of the spectrograph in §4, our revised telescope control software in §5, our one-dimensional extraction in §6, our Doppler pipeline in §7, our first RV results in §8, and our final remarks in §10.

2 Survey Performance

Observing at Mt Hopkins is divided naturally into seasons by the July/August monsoon shutdown. The first full-season MINERVA observing campaign in radial velocity survey mode began 2017 September 14 and ran through 2018 June 29. Spectra were obtained on 196 of 293 nights. Weather prevented observations on 44 nights, and 53 nights were spent on engineering or lost to system malfunctions. We obtained 1936 exposures with 4 spectra each of 28 survey target stars, with a maximum of 222 exposures of a single (high decl.) target. Fourteen targets had at least 60 exposures. In addition, we obtained 199 exposures of hot stars, at least one per night, used for spectral calibration. A typical night full of observing led to 12 to 19 exposures (most in December, less on shorter nights) of 8 to 10 target stars. The open shutter fraction was highly variable at the beginning of the season, but stabilized at $\sim 69$% after implementation of the autonomous scheduler in late October. Given the rapid slewing and settling time of our telescopes, the majority of the overhead per spectrum was the result of robotic target acquisition on the fiber tip.

The 2018-19 observing campaign began 2018 October 15 and is in progress at the time of writing this manuscript. Through 2019 March 31, spectra have been obtained on 107 of 168 nights, with 35 nights lost to weather and 26 spent on engineering or lost to system malfunctions. We have obtained 1455 exposures with 4 spectra each of 19 survey target stars, a 32% increase in spectra obtained over the same period from the previous season. Twelve targets have at least 60 exposures thus far. In addition we have obtained 137 exposures of hot calibration stars. Changes in our acquisition algorithm have reduced the overhead, resulting in an average open shutter fraction of 86% since 2018 November. Historically, April through June provide very reliable weather at the site (in 2018-2019 we lost only 3 of these 91 nights to weather), and we anticipate that this observing season will lead to a larger set of RV data than the previous season — 43% of our 2017-18 spectra were obtained in April–June.

3 Hardware

The overall hardware design for the MINERVA facility has remained largely unchanged from that described in Swift et al. (2015). However, we have made several changes to improve our science capability, which we discuss in detail below.

3.1 Fiber Acquisition Unit (FAU) Cameras

The SBIG ST-i cameras originally used in the FAUs had two major problems. First, their small field of view (3${}^{{}^{\prime}}$.6 x 2${}^{{}^{\prime}}$.7) 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 x 1216, 5.86 $\mu m$ pixels, and is similarly priced. This provides us an 8.5’ x 5.4’ 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.

3.2 Fiber

The science fiber was originally purchased from CeramOptec and has a 50 $\mu$m octagonal core and 330 $\mu$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 $\mu$m circular core and a 125 $\mu$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 $\mu$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 $\mu$m octagonal core and 110 $\mu$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 ($\sim 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.

3.3 Spectrograph

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 Optics111http://www.primeoptics.com.au/. 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 $\lesssim 0.1\%$.

We determined our total system throughput to be $\sim 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-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 $\mu$m octagonal fibers at $f$/6.6 in our FAU (see §3.1). Three meters before the spectrograph entrance, each of the fibers are butt-coupled to circular 50 $\mu$m fibers, which are then arranged into a V-groove at the entrance to the spectrograph, separated by 220 $\mu$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+(\lambda-5500$Å$)\times 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.

3.4 Exposure Meter

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 $\sim$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.

3.5 Backlight

The FAU design described in §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 $\sim 10$ $\mu$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.

3.6 Slit flat

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 §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 2a and 2b 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.

4 Spectrograph environmental performance

Here we show the pressure stability of the spectrograph (Figure 3a) from 2017 March through 2018 July at 3 s intervals. It should be noted however that during 2017 mid-April, a power outage occurred at the MINERVA facility. Typically, the spectrograph is continuously pumped, but the outage caused the valve from the pump to the spectrograph to close, which went unnoticed for an extended period of time (mostly during the monsoon when we were not operational). We have since implemented a watchdog that sends an email notification if the pressure rises above 10 $\mu$bar (see §5). This power outage resulted in a swift increase in pressure as it leaked up toward atmosphere. By early October, the pressure once again became sufficiently stable for the collection of good quality data, and remained so through 2018 July aside from some minor fluctuations due to maintenance. Figure 3b shows this stability over the month of 2018 March with a RMS of 0.065 $\mu$bar (dramatically exceeding our requirement of 12 $\mu$bar).

Temperature readings meanwhile, were recorded for MINERVA over the period of 2018 January through July (Figure 4a). These measurements were taken at the side of the cell holding the Echelle grating, and so it is most relevant for RV stability. Throughout that time period we manage to remain fairly stable from January to May, one such example being Figure 4b which depicts the stability over the month of March, with an RMS of 0.0052 K (two times better than our requirement of 0.01 K). Slightly larger fluctuations were seen to occur from May onward, where we removed and reinstalled the outer thermal enclosure for maintenance.

Moving forward, we intend to investigate the causes of some of the more minor fluctuations present in the environmental data. These could be a result of events such as the backlight being turned on and off, moving of the iodine stage in and out (which holds the cell which must be heated to 55 ∘C), fluctuations of the room HVAC, or other events related to the operation of MINERVA itself.

The scatter in the empirically determined wavelength solution of a single chunk for all our targets on sky (0.003 Å, or about 0.1 pixels, for the zero-point and 0.07% for the dispersion) dominates any trends seen on nightly or monthly timescales.

5 Telescope Control Software

5.1 Architecture

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 Github222https://github.com/MinervaCollaboration.

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 §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.

5.2 Operations

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 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.

6 Spectroscopic Data Reduction

The first step in our spectroscopic data reduction is to calibrate the science exposures of our spectrograph’s CCD. We collect and median stack eleven frames each night for the bias, dark current, and slit flats. For each science exposure, we subtract the overscan from the raw exposure. We experimented with dark current subtraction but omit this in our present pipeline because the corrections are negligible and it only serves to increase the noise. We re-normalize each bias-corrected exposure by the stacked slit flats, similar to the procedure developed in Bernstein et al. (2015), although we retain the blaze function. Finally, we interpolate between fiber bundles to estimate and subtract scattered light.

With our calibrated science frames, we are prepared to extract the one-dimensional spectrum from the two-dimensional CCD exposure. We wrote a custom pipeline using the optimal extraction algorithm (Horne, 1986; Piskunov & Valenti, 2002; Zechmeister et al., 2014; Bernstein et al., 2015). Optimal extraction requires that flux is a separable function of $x$ and $y$ so that $F(x,y)=F(x)F(y)$, a condition that is very nearly satisfied in our instrument. This allows us to independently find the observed flux in each row, $x$, through

[TABLE]

Here $F(x)$ is the underlying spectrum we wish to extract. We determine this from the observed flux in the cross-dispersion direction, $F(y)$, a model for the noise $n(y)$, and a normalized cross-dispersion profile, $\bm{p}(x,y)$.

Our pipeline presently uses a modified Gaussian for $\bm{p}(x,y)$. This gives us the form

[TABLE]

The free exponent $p(x)$ is slightly broader than a typical Gaussian with $p\approx 2.2$. The value $N(x)$ is a numerically determined normalization coefficient and $y_{c}(x)$ indicates the trace centroid, determined during calibration from archival fiber flats. We model $\sigma(x)$ and $p(x)$ as slowly varying polynomials along the dispersion direction

We simultaneously extract all fibers within each column, accounting for any cross-talk. During extraction, we include a slowly varying background term to account for any additional scattered light. We also apply a cosmic ray rejection algorithm and mask any hits. Although our precise wavelength solution $\lambda(x)$ is found with the Doppler pipeline, we generate an initial solution from archived thorium argon exposures we took during the installation and maintenance of the spectrograph. This allows the subsequent code to quickly lock on to the correct solution.

The spectra from individual telescopes are extracted and modeled independently all the way through to the orbital modeling. This gives us four data points per exposure and a unique insight into systematic errors having to do with the telescope and the position of the trace on the detector (cosmic rays, scattered light, and flat-fielding).

7 Doppler Pipeline

The one-dimensional spectrum is the primary input for our Doppler code333https://github.com/MinervaCollaboration/minerva-pipeline. The architecture and general principles of our Doppler code are inspired by the code that is comprehensively described by Wang (2016), although the algorithm is originally introduced by Butler et al. (1996). Our code implements a forward modeling procedure on this spectrum that can be summarized mathematically as

[TABLE]

where $x$ is the pixel position in the dispersion direction, $\lambda(x)$ is the wavelength solution, $\lambda^{{}^{\prime}}(x)$ is the Doppler-shifted wavelength solution, $F_{\rm obs}$ is the one-dimensional spectrum extracted from our observations, $F_{I_{2}}$ is the normalized absorption spectrum of our iodine cell, $F_{\star}$ is the stellar flux, and IP$(x)$ is a model of the spectrograph’s intrinsic instrumental profile (which is sometimes referred to as the spectrograph response function or the one-dimensional spectral point spread function). After determining the product of the iodine absorption spectrum and stellar spectrum, the observed spectrum is modeled as this product convolved with the instrumental profile.

7.1 Iodine Absorption Spectrum

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 $\lambda(x)$ because the resolving power of the FTS ($R\approx 300,000$) is about a factor of 4 greater than that of our KiwiSpec spectrograph ($R\approx 80,000$). Because the molecular iodine lines span from 500 to 630 nm, the wavelength solution is determined solely within this range.

7.2 Choice of IP$(x)$

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 §7.5 and the other is a time-varying summation of satellite Gaussian profiles stacked on one central Gaussian profile that is described in §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.

7.3 Reference Stellar Spectrum

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_{\star}(\lambda_{\rm ref}(x))$ is solved for via $F_{\rm obs,ref}(x)=F_{\star}(\lambda_{\rm ref}(x))*\textrm{IP}(x)$. This IP$(x)$ and $\lambda_{\rm 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 $\lambda_{\rm 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_{\rm obs,ref}(x)$ to get $F_{\star}(\lambda_{\rm ref}(x))$, whose wavelength solution is assumed to be the averaged $\lambda_{\rm 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_{\rm 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.

7.4 Doppler-shifted Wavelength Solution

The forward modeling procedure finds the best fit to $F_{\rm 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_{\star}(\lambda^{{}^{\prime}}(x))$. The Doppler-shifted wavelength solution can be decomposed as $\lambda^{{}^{\prime}}=\lambda_{\rm ref}\cdot(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 ($\sim 0.5$ km s*-1*) and orbital motion around the Sun ($\sim 30$ km s*-1*). 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.

7.5 Fixed IP and IP Stability

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 $\textrm{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 time scales at which our instrumental profile is stable. To find this time scale, 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.

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,

[TABLE]

where $A$ is the amplitude, $x$ is the pixel position in the cross-dispersion direction, $x_{0}$ is the centroid, $\sigma$ 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_{\rm norm}(x)=(D_{\rm raw}(x)-b)/N$, where $D_{\rm norm}$ is the normalized background-subtracted crosscut, $D_{\rm raw}$ is the original crosscut, and $N$ is the normalization factor calculated by the aforementioned integration.

Once $D_{\rm 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_{\rm norm}(x))$. To find the optimal fixed IP, we perform an iterative process of modeling the $D_{\rm raw}(x)$ and subsequently evaluating a new $D_{\rm norm}(x)$ and IP${}_{f}(x)$.

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

[TABLE]

where $M(x)$ is the model for one crosscut and $\Delta x$ is a translational shift parameter. Now, the least-squares optimizer has only three parameters to evaluate: $\Delta 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 $\chi^{2}$ of all the crosscuts’ data and $M(x)$ models collectively, and lastly evaluate the difference between the previous iteration’s reduced $\chi^{2}$ and the current iteration’s reduced $\chi^{2}$. The most important distinction between iterations is the differing fixed IPs; when the $\chi^{2}$ 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 $\chi^{2}$ values lessens. Once this difference is less than $10^{-4}$, 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 6a is an example of the final fixed IP. The blue data points represent the $D_{\rm 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_{\rm norm}(x)-\textrm{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.

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 §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 6b and it has the same general pattern of noise in its residuals as Figure 6a. This implies that the instrumental profile has not changed within that time period. Note that for Figure 6b the 2018 fixed IP was used to model the 2017 spectra (via Eq. 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 September 14 and October 2 in 2017, 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 6c 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 6a and 6b. 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 §4 and shown in Figure 3a, 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 6c is likely to occur.

7.6 Gauss-Hermite IP

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 6c and explained in Figure 3. We can formulate such an IP as

[TABLE]

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 6c 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.

7.7 Sum-of-Gaussians IP

Our sum-of-Gaussians IP,

[TABLE]

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_{0}$ 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 IPG 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 §8.

8 RV Performance

The radial velocity we measure is the reflex motion of the star induced by the gravitational pull of its planetary companion. This motion is accounted for in Kepler’s laws. Kepler’s laws suggest that the lower limit of the planetary mass can be described as

[TABLE]

(see parameter symbols with Table LABEL:table:_51_Peg descriptions). To model the stellar system, we use EXOFASTv2 (Eastman et al., 2013; Eastman, 2017).

Before the minimum mass is determined, the RV semi-amplitude must be extracted from the Doppler-shifted spectra via a forward modeling procedure. The results of the instrumental profile work discussed in §7.5 provided fruitful information that paved the way for the successful extraction of MINERVA’s first radial velocity results.

We present RV measurements of two target stars to demonstrate MINERVA’s precision. These stars are HD 122064 and HD 217014 (51 Peg), which have their coordinates, $V$ magnitude, and spectral type reported in Table LABEL:table:_stars. HD 122064 is chromospherically inactive, has no known companions, and serves as a convenient RV standard star. For the hot Jupiter 51 Peg b, we compare the planetary properties derived from MINERVA data with results from the literature.

8.1 HD 122064

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 §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 $\sim$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 7a while the sum-of-Gaussians IP is used to produce the RVs of Figure 7b.

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).

We have seventy-five radial velocity measurements tabulated in Table LABEL:table:_RV_standard_star_RVs and shown in Figures 7a and 7b. Figures 8a and 8b 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*-1* level for our measurements of this RV standard star. The precision achieved is 1.8 m s*-1* for the fixed IP approach and 1.9 m s*-1* 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.

8.2 HD 217014

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 broad band photometry summarized in Table LABEL:table:_mags. 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 9a, tabulated in Table LABEL:table:_51_Peg_RVs, and summarized in Table LABEL:table:_51_Peg. Figure 9b shows the same but with the time series folded to the phase of the planet’s orbital period.

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 LABEL:table:_51_Peg. We only employ constraints on three of the stellar parameters. We impose a prior on the stellar metallicity ([Fe/H]) of 0.20 $\pm$ 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 ($\pi$) of 64.65 $\pm$ 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 Hipparcos444Vizier Online Data Catalog, I/239 (ESA 1997) catalog. Valenti & Fischer (2005) suggest that the typical uncertainties for their stellar parameters amongst their catalog of stars are 0.06 dex for log $g$, 44 K for $T_{\rm 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 LABEL:table:_51_Peg 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 $\sigma_{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 LABEL:table:_51_Peg. 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$\sigma$ uncertainty. We define this discrepancy as

[TABLE]

in which $\sigma_{1,L}$ is the lower error bar for $N_{1}$ and $\sigma_{2,U}$ is the upper error bar for $N_{2}$. Ideally, the discrepancy should be less than 1$\sigma$. Seeing as the discrepancies evaluated for the tabulated parameters—most importantly the planetary parameters—are $\lesssim 1.0\sigma$, 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*-1*) and systematic floor (4.2 m s*-1*) 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*-1*) and an insignificant decrease in the binned-by-2 residual RMS (4.3 to 4.2 m s*-1*). Therefore, we conclude that fitting the period and eccentricity does not have a major impact on our inferred precision.

9 Future Work

While we are about a factor of 2 of our original goal and already operating with an unmatched combination of cadence and precision that enables us to detect new planets and provide valuable insight into stellar jitter, there are several areas for improvement that may help us achieve our original goal of 80 cm s*-1*.

The relatively long exposures we typically take increases the uncertainty in the flux-weighted midpoint time and therefore introduces additional error in the corresponding barycentric correction (Wright & Eastman, 2014). Either shortening our exposures or improving the determination of the flux-weighted midpoint may improve our ultimate RV precision.

The stability of our spectrograph demonstrated here allows us several avenues to improve our Doppler pipeline. The templates we use, derived from Keck/HIRES observations, is likely to contain systematic differences due to the instrument and atmospheric conditions (e.g., increased water column) from our instrument, and may be dominating our RV error. We will investigate generating our own templates, but we will also investigate fitting for the template from dozens of spectra, known as the “grand solution” (Gao et al., 2016; Czekala, 2017). Further, modeling many spectra at once with the same instrumental profile and/or fitting for the iodine cell removes many sources of potential systematic error and unnecessary nuisance parameters that may be covariant with the radial velocity.

While we have octagonal fibers coupled to a circular fiber for optimal near-field scrambling, we have opted not to introduce a double scrambler to improve our far field scrambling, which typically reduces the throughput by 10-20% (Halverson et al., 2015). We may revisit this trade-off in the future, as well as explore the possibility of introducing an agitator to improve modal noise. We are also actively exploring the improvement achievable by using spectro-perfectionism to improve the extraction of our 1D spectra using a 2D instrumental profile (Bolton & Schlegel, 2010).

Finally, having four simultaneous spectra gives MINERVA a unique insight into telescope and detector level systematics which we have yet to fully capitalize on. In particular, we will explore the possibility of adapting the “vanking” stage of the Doppler pipeline, which is a sophisticated outlier rejection algorithm, to incorporate the knowledge that the four simultaneous spectra should produce identical radial velocities.

10 Conclusion

Since the commissioning of MINERVA, we have substantially modified our telescope control software and our Doppler pipeline. The MINERVA mission’s secondary goal is accomplished much more efficiently with the control software changes. This work marks our first achievement toward our primary science goal of obtaining precise RVs. We have confirmed which of the IPs we had at our disposal would yield reliable results from our pipeline. These are the aforementioned time-invariant cubic spline function (the fixed IP) and the sum-of-Gaussians function. While testing the fixed IP, we have also confirmed that our spectrograph’s intrinsic instrumental profile remains stable for months. When there is significant fluctuation in the intrinsic instrumental profile, it is likely due to disturbances to the instrument, as opposed to any natural and gradual perturbations within the instrument. The agreement between both IPs implies that using an instrumental profile from the cross dispersion direction, and therefore has systematic errors, is sufficient if the instrumental profile is stable. Consequently, we precisely characterized our spectrograph’s instrumental profile from the cross-dispersion direction of the echellogram.

MINERVA is a collaboration among the Harvard-Smithsonian Center for Astrophysics, The Pennsylvania State University, the University of Montana, and the University of Southern Queensland. MINERVA is made possible by generous contributions from its collaborating institutions and Mt. Cuba Astronomical Foundation, The David & Lucile Packard Foundation, National Aeronautics and Space Administration (EPSCOR grant NNX13AM97A), The Australian Research Council (LIEF grant LE140100050), and the National Science Foundation (grants 1516242 and 1608203). Any opinions, findings, and conclusions or recommendations expressed are those of the author and do not necessarily reflect the views of the National Science Foundation. Funding for MINERVA data-analysis software development is provided through a subaward under NASA award MT-13-EPSCoR-0011. This work was partially supported by funding from the Center for Exoplanets and Habitable Worlds, which is supported by the Pennsylvania State University, the Eberly College of Science, and the Pennsylvania Space Grant Consortium. We are grateful to Dr. Gillian Nave and R. Paul Butler for providing FTS measurements of our iodine gas cell.