A Decade of Radial-velocity Monitoring of Vega and New Limits on the Presence of Planets

We obtained high-resolution spectra of Vega with the Tillinghast Reflector Echelle Spectrograph (TRES; Fűrész 2008), which is mounted on the 1.5-m Tillinghast Reflector at Fred L. Whipple Observatory on Mount Hopkins, AZ. It has a resolving power of R ∼ 44, 000, and a wavelength coverage of 3850–9100 Å. We obtained a total of 1524 spectra spanning the ten-year period between UT 2009 June 13 and 2019 October 24. Typical exposure times ranged from a few seconds to a few tens of seconds, achieving signal-to-noise ratios (S/Ns) between about 300 and 1000 per resolution element.

We optimally extracted and reduced the spectra following the procedures outlined in Buchhave et al. (2010), and while we begin by following the radial-velocity measurement outlined therein, our final velocity extraction includes a few key differences. We first cross-correlate each spectrum of Vega against the strongest exposure, treating each spectral order separately. The relative RV for each exposure is taken to be the location of the peak of the summed cross-correlation function (CCF; across all orders) from that spectrum. The internal RV uncertainty for each observation is taken to be the standard deviation of the locations of the CCF peak of each order for that spectrum. Next, we shift and median combine the 1524 spectra to generate a master template spectrum. In the case of Vega, the template S/N is 26,000 per resolution element. Though exposures of Vega are typically short enough that cosmic-ray rejection is not very important, we identify outlier pixels and replace them with the median spectrum at that location. Finally, we re-run the order-by-order cross-correlation, this time against the high-S/N template.

TRES was not designed for long-term stability at the level of meters per second, and has at times experienced drifts and jumps in its instrumental zero-point as large as a few tens of meters per second. Changes of this magnitude can mimic or mask the presence of long-period companions. To combat this problem, we track the zero-point with nightly observations of several RV standard stars, allowing us to measure and correct for zero-point changes over time. From the standard deviation of the RV standards, we also estimate the instrumental RV noise floor, to be added in quadrature with the internal error estimates described above. The TRES instrumental precision at the beginning of our data set was ∼50 m s⁻¹, but by the end of 2010 had improved to ∼10–15 m s⁻¹, thanks in large part to hardware upgrades. Though these introduced some zero-point changes, they are corrected for in the same way as other zero-point shifts. Our final, zero-point-corrected, relative RVs are presented in Table 2, and shown in Figure 1.

Zoom In Zoom Out Reset image size

Vega was observed by NASA's TESS mission (Ricker et al. 2015) in Sector 14, between UT 2019 July 18 and 2019 August 15, and in Sector 26, between UT 2020 June 8 and 2020 July 04. The star is bright enough to fill the full well depth of hundreds of pixels, but the TESS detector is designed to preserve the flux by spilling into neighboring pixels. As long as the star is not too close to the edge of the detector, the full frame images (FFIs), returned at 30 minute cadence, can be used to extract photometry from an area encompassing all of the light. The FFIs were processed by the Science Processing Operations Center (SPOC) at NASA Ames (Jenkins et al. 2015, 2016), adapted from the pipeline originally developed for the Kepler mission (Jenkins et al. 2010). We then performed photometry using apertures shaped to trace the distribution of charge on the TESS images, including 3625 pixels in Sector 14 and 4335 pixels in Sector 26. We corrected the flux for sources within the aperture using TESS magnitudes listed in the TESS Input Catalog. To account for systematics caused by the motion of the spacecraft, we followed the procedure outlined in Vanderburg et al. (2019). Namely, we assumed that systematics caused by changes in spacecraft pointing can be corrected by decorrelating against the background scattered light signal and combinations of the spacecraft quaternions, which encode the pointing of the spacecraft every 2 s. We produced quaternion time series by calculating the standard deviation and mean quaternions during each 30 minute exposure. Ultimately, we found the best performance (i.e., lowest resulting light-curve scatter) when decorrelating against the scattered light signal and the first, second, and third order of these quaternion time series. We exclude a few hours of data most strongly affected by scattered light at the beginning of each orbit in Sector 26, and we use a spline to flatten remaining low-frequency systematics that occur on the timescale of the spacecraft orbit.

The resulting light curve is shown in Figure 2, and is remarkably quiet. The median standard deviation in 6.5 hr bins is only 5.0 ppm. We run Transit Least Squares (TLS; Hippke & Heller 2019) in search of transiting planets, but we find no evidence for transit-like features. In Section 4.2, we produce detection limits using transit injection recovery. The maximum power in a Lomb–Scargle periodogram occurs for a period of 0.68 day (Figure 2, lower left panel), which is consistent with previous measurements of the rotation period of Vega. We do detect the signal in both sectors individually at lower significance, but because of the year-long gap between sectors, we cannot confirm whether the signal is stable in phase for the full time span. It is interesting to note that while the period is still significant (>3σ), if we exclude photometric outliers, the significance is highest when they are included. This suggests that the apparent outliers vary in phase with the rest of the data and may be associated with stellar activity rather than systematics related to the instrument or data processing. This could indicate that some surface features are evolving in brightness on very short timescales, while others appear more stable over the course of a month, and perhaps over the full year spanned by the TESS data.

Zoom In Zoom Out Reset image size

3.1.1. Detection of the Activity Signal

One of the strongest peaks in the RV periodogram is located at 0.677 day, which corresponds to the previously reported rotation period (Alina et al. 2012; Böhm et al. 2015). Table 3 lists all of the significant peaks in the periodogram, and shows how they relate to the rotation period. Nearly every signal is an alias of the 0.677 day periodicity, with harmonics also present but at lower significance. While we have good reason to suspect that this value corresponds to the true signal, we conduct a more rigorous test of the aliases using the methodology of Dawson & Fabrycky (2010). This procedure simulates a sinusoidal signal at the candidate period using the time stamps of the observed data. The amplitudes and phases of each peak in the periodogram are then compared between the synthetic and actual data. The simulated signal that best reproduces the results from the observations is chosen as the true period. While the signal we observe is likely caused by activity, and therefore is somewhat irregular, this test does confirm 0.677 day as the best match to the observed periodogram.

Table 3. Significant Periodogram Signals (FAP < 0.1%)

Period	Frequency	Notes	Amplitude	FAP
(days)	(days−1)		( m s−1)
0.105	9.48	Frot + 8	7.16	7.71e-6
0.111	9.00	⋯	7.12	2.49e-5
0.118	8.48	Frot + 7	8.21	3.80e-8
0.134	7.48	Frot + 6	7.60	3.37e-6
0.154	6.48	Frot + 5	7.59	7.63e-6
0.182	5.48	Frot + 4	7.58	1.06e-6
0.221	4.52	$\left|{F}_{\mathrm{rot}}-6\right|$	8.31	2.35e-8
0.223	4.48	Frot + 3	9.15	8.23e-11
0.263	3.80	⋯	6.88	9.05e-5
0.284	3.52	$\left|{F}_{\mathrm{rot}}-5\right|$	8.17	5.79e-9
0.288	3.48	FFrot + 2	9.43	1.06e-10
0.397	2.51	$\left|{F}_{\mathrm{rot}}-4\right|$	7.40	4.36e-7
0.403	2.48	Frot + 1	8.72	8.48e-10
0.677	1.48	Frot	9.06	9.86e-10
2.10	0.478	Frot − 1	7.29	2.31e-5

Download table as:  ASCIITypeset image

A useful diagnostic for distinguishing between activity and orbiting companions is the stacked periodogram, which is described by Mortier & Collier Cameron (2017). As the number of observations included in a periodogram increases, the power of a peak corresponding to a planetary signal should monotonically increase, within the limits of the noise, while the power of a peak corresponding to rotation of features on the stellar surface may not, due to changes in the phase of the signal driven by the evolution of active regions. In Figure 4, we present a stacked periodogram for our TRES radial velocities. For the first ∼1000 observations, the power at 0.677 day increases, but afterwards, the power falls off. This supports the conclusion that the periodicity originates from activity, and its variation over the course of hundreds of observations is consistent with a long timescale of evolution.

Zoom In Zoom Out Reset image size

3.1.2. Modeling the Activity Signal

We first attempt to model the stellar activity with a Keplerian orbit using the radvel Python package (Fulton et al. 2018). However, Markov Chain Monte Carlo (MCMC) samples failed to converge. Additionally, a periodogram of the residuals to the best fit showed significant power remaining near 0.677 day. This suggests that a Keplerian model does not adequately fit the signal, which is expected given the structure seen in the stacked periodogram. It would have been surprising for the stellar surface to be described well by a single Keplerian (corresponding to a single active region) over the span of a decade.

We next consider a model in which active regions are very long-lived, but multiple regions may exist with slightly different periods due to differential rotation. For this exercise, we follow an iterative whitening procedure in which we fit a Keplerian with initial conditions determined from the periodogram following Delisle et al. (2016). After refining that fit, we repeat the process with the strongest peak in the periodogram of the residuals, until there are no signals remaining with FAP < 0.1%. This procedure returns only two signals, 0.676085 and 0.676684 day, both coincident with the rotation period. This makes sense given that the same hemisphere of the star is almost always visible. Very low-latitude (equatorial) active regions, visible for only a portion of the rotational phase, may appear at harmonics of the rotation period. Otherwise, variability caused by static surface features on Vega should be well described by a sum of Keplerians at the rotation period, with small differences in period caused by differential rotation. In Table 4 we present the parameters of the Keplerians used to model these two dominant signals. Some power remains in the periodogram of residuals near the rotation period, its aliases, and its harmonics, which may indicate the presence of additional surface features. However, it may also indicate that the surface of Vega is not stable on decade timescales, in which case we need to apply a model that can better address its time variability.

Table 4. Stellar Activity RV Models

Parameter	Units	Value
Multiple Keplerians
P1	days	0.676682 ± 0.000006
Tc1	BJD	2457187.123 ± 0.013
e1	⋯	0.438 ± 0.052
ω1	deg	344 ± 12
K1	m s−1	11.6 ± 1.2
P2	days	0.676086 ± 0.000006
Tc2	BJD	2457187.195 ± 0.051
e2	⋯	0.46 ± 0.22
ω2	deg	179 ± 15
K2	m s−1	11.1 ± 3.8
γ	m s−1	-33.69 ± 0.69
σjit	m s−1	13.1 ± 1.0
Quasiperiodic GP
η1	m s−1	${14.4}_{-1.8}^{+2.2}$
η2	days	${179}_{-58}^{+65}$
η3	days	${0.6764}_{-0.00021}^{+0.13}$
η4	⋯	${0.44}_{-0.15}^{+0.16}$
γ	m s−1	$-{33.9}_{-2.8}^{+2.9}$
σjit	m s−1	12.9 ± 1.4

Download table as:  ASCIITypeset image

Quasiperiodic Gaussian processes (GPs; Rasmussen & Williams 2006; Roberts et al. 2012) are known to represent the effects of active regions rotating in and out of view (Haywood et al. 2014; Rajpaul et al. 2015; Angus et al. 2018). Using radvel, we model the correlated stellar noise with a kernel matrix whose elements are defined as

$$\begin{eqnarray}&&{C}_{{ij}}={\eta }_{1}^{2}\exp \left(-\displaystyle \frac{\left|{t}_{i}-{t}_{j}\right|}{{\eta }_{2}^{2}}-\displaystyle \frac{\sin \left(\tfrac{\pi \left|{t}_{i}-{t}_{j}\right|}{{\eta }_{3}}\right)}{2{\eta }_{4}^{2}}\right),\end{eqnarray} \tag{ 1 }$$

where t_i and t_j are observations made at any two times; η₁ is the amplitude hyperparamter; η₂ is the exponential decay timescale, which is physically related to the lifetime of the active regions; η₃ is the period of the variability, corresponding to the rotational period; and η₄ is the length scale of the periodic component, describing the high-frequency variation in the stellar rotation. It is important to note, however, that these physical interpretations are often not straightforward, particularly given the degeneracies between the hyperparameters.

No trend is apparent in our data or the residuals of early fits, and we choose to exclude a linear or quadratic term in our model. However, we do include offset and jitter terms to account for instrumental offset and noise. Priors are only placed on parameters to keep them within physically possible limits. The complete log-likelihood of this model is

$$\begin{eqnarray}&&\mathrm{ln}{ \mathcal L }=-\displaystyle \frac{1}{2}{r}^{T}{K}^{-1}r-\displaystyle \frac{1}{2}\mathrm{ln}\left(\det K\right)-\displaystyle \frac{n}{2}\mathrm{ln}\left(2\pi \right),\end{eqnarray} \tag{ 2 }$$

where r is the vector of residuals, K is the covariance matrix, and n is the number of data points.

We perform an affine-invariant MCMC exploration of the parameter space using the ensemble sampler emcee (Foreman-Mackey et al. 2013, 2019). Our MCMC analysis used eight ensembles of 50 walkers and converged after 4,950,000 steps, achieving a maximum Gelman-Rubin statistic (Gelman & Rubin 1992) of 1.006. The resulting posterior distributions are shown in Table 4. We do note that this model is over-fitting the data, with a reduced χ² statistic of 0.869; either the model is fitting noise or our cadence is not sufficient to constrain high-frequency variations in the rotation signal. In either case, the GP is likely too complex of a model for our data. We check our results with a second GP fit of the radial velocities using a quasiperiodic kernel with the PyMC3 Python package (Salvatier et al. 2016), which returns results well within 1σ of those in Table 4.

The activity signal at the period of stellar rotation is by far the strongest we observe, but other candidate signals have been reported in previous Vega data sets, and we also search for signals of lower significance in the TRES RVs.

3.2.1. Previously Suggested Signals

3.2.2. Signals at Long Periods

3.2.3. A Candidate Planetary Companion

Interestingly, after removing two Keplerian signals associated with the stellar rotation, the strongest remaining peak has a period of 2.43 days and a formal FAP of only <0.01. Unlike the 194 day signal, it is robust to the inclusion of additional Keplerians near the rotation period and its harmonics; the FAP of the 2.43 day signal remains below 1% whether we account for stellar noise or not and even if we include more complex static models of the stellar surface with as many as five Keplerian components. We do not find any relationship between this period and the rotation period and its aliases listed in Table 3, and the uncertainties on the periods are much smaller than the differences between them. A circular orbit with a period of 2.43 days has an RV semi-amplitude of about 6 m s⁻¹, corresponding to a minimum mass of about 20 M_⊕. According to the difference in the Bayesian Information Criterion (ΔBIC), a circular orbit is statistically preferred, but we present both eccentric and circular solutions in Table 5. Figure 5 shows the phase-folded RVs of the candidate planetary signal. While we cannot conclusively rule out false-positive scenarios, we discuss in Section 5 ways in which we might confirm the candidate with future observations and analyses.

Zoom In Zoom Out Reset image size

Table 5. Candidate Planetary Companion to Vega

Parameter	Units	Value
		Eccentric	Circular
P	days	2.42977 ± 0.00016	2.42971 ± 0.00018
Tc	BJD	2457186.51 ± 0.18	2457186.631 ± 0.067
e	⋯	0.25 ± 0.15	0
ω	deg	304 ± 40	⋯
K	m s−1	6.4 ± 1.1	6.1 ± 1.1
γ	m s−1	−34.30 ± 0.75	−34.34 ± 0.75
σ	m s−1	9.7 ± 1.0	9.7 ± 1.0
${m}_{p}\sin i$	M⊕	21.9 ± 5.1	21.2 ± 3.8

Download table as:  ASCIITypeset image

4.1.1. Isotropic Orbits

We begin the calculation of our detection limits by randomly generating ∼10⁶ models corresponding to planets with semimajor axes between 0 and 15 au and masses ranging from 0 to 100 M_J. Each orbit is assigned an inclination (drawn from a uniform distribution in $\cos i$), eccentricity (drawn from a beta distribution described by Kipping 2013, with parameters a = 0.867 and b = 3.03), an argument of periastron (drawn from a uniform distribution), and a time of periastron passage (drawn from a uniform distribution). With radvel, we calculate the expected radial velocities for each orbit at our times of observation. We then add noise scaled to the observed uncertainties at each time stamp. Using radvel, we fit a flat line to the synthesized RVs and calculate the χ² statistic and its associated p-value for this model. Low p-values indicate that the synthesized signal is unlikely to arise from white noise—i.e., it is distinguishable from a flat line, and we therefore consider it detected. To set a threshold for detection, we follow Latham et al. (2002), who demonstrated that p < 0.001 is a conservative threshold below which signals are almost entirely real. This makes sense, since data comprising only white noise will exhibit p < 0.001 only 0.1% of the time. While we find that other methods of injection recovery—such as a requirement that ΔBIC between the detected signal and a flat line model be greater than 10—yield more sensitive limits, we adopt the conservative p(χ²) < 0.001 because we do not address correlated noise in our simulations. In the likely event that we cannot perfectly model the stellar activity, it may hinder our ability to detect some signals. The top panel of Figure 6 shows the distribution of detection probabilities for isotropically oriented orbits, from which it can be seen that we are sensitive to sub-Saturn masses orbiting at 0.1 au and masses as low as about 2 M_J at 6 au, which corresponds to an orbital period of 10 yr around Vega. Beyond 6 au, our detection limits degrade more quickly, as our data no longer cover an entire orbit. Nonetheless, we are sensitive to the most massive giant planets, all the way out to 15 au. Given the roughly isotropic distribution of stellar obliquities for hot stars hosting transiting hot Jupiters, the decision to draw inclinations from an isotropic distribution may be the most realistic assumption for short periods. However, it is unclear whether the same is true at long periods. We therefore also explore our detection limits for well-aligned orbits in the following section.

Zoom In Zoom Out Reset image size

4.1.2. Well-aligned Orbits

4.1.3. Including Direct Imaging Limits

Using the Python package batman (Kreidberg 2015), we randomly generate ∼10⁵ transit models for planets with periods ranging from 0.5 to 30 days and radii between 1 and 8 R_⊕. Each orbit is assumed to be circular, and semimajor axes are calculated using Kepler's law, assuming that the planet's mass is negligible compared to Vega. Transit times are randomly assigned from a uniform distribution, and inclinations are drawn from a uniform distribution in $\cos i$ such that

$$\begin{eqnarray}&&0\lt \cos i\lt \displaystyle \frac{{R}_{\star }-{R}_{p}}{a},\end{eqnarray} \tag{ 3 }$$

where a is the semimajor axis and R_⋆ and R_p are, respectively, the stellar and planetary radii, ensuring that the planet transits. Additionally, each synthetic transit follows quadratic limb-darkening laws for an A0 star in the TESS bandpass with u₁ = 0.1554 and u₂ = 0.2537 (Claret 2018).

We then inject each model into our light curve and conduct a transit search using TLS. We consider a transit to be detectable if the best period recovered by TLS is within 1% of the injected period and if at least one of the transit times is within the same margin of an injected transit. Additionally, because a transit could appear within the gaps of our light curve, any recovered integer multiples of the injected period are considered detectable.

We also require a transit to be distinguishable from a false positive in order for it to be considered a true recovery, ensuring that the signal would not be dismissed due to a low signal in a real search. Using synthetic light curves and transits, Hippke & Heller (2019) find that a signal detection efficiency (SDE) threshold of 7 corresponds to a false-positive rate of 1% for TLS. Therefore, for a transit to be detectable, we also require that the best period have an SDE of 7 or greater.

The results of our injection recovery test are shown in Figure 7. At the shortest periods, we are sensitive to transiting planets as small as 2 R_⊕, and as small as 3 R_⊕ for orbital periods similar to the duration of a TESS orbit (≲14 days). For periods longer than a couple weeks, our formal sensitivity drops, as fewer transits are expected; some planets may exhibit only one or two transits total, and some might not even be observed once by TESS due to data gaps between orbits. Many of these will not pass our automated TLS criteria for detection. On the other hand, the transits of even relatively small planets orbiting Vega would be easily visible by eye in the quiet TESS light curve. A 4 R_⊕ transit, for example, would be about 180 ppm. It is clear that there are no Neptune-sized planets transiting in the TESS data even once. While the formal detection limits shown in Figure 7 are illustrative of the types of transiting planets we could detect most easily, they should be taken as a conservative estimate.

Zoom In Zoom Out Reset image size

Rotational modulation has consistently been found in spectroscopic observations of Vega (Böhm et al. 2015; Petit et al. 2017), providing evidence for active regions on the surface of the star. While chemically peculiar A stars are known to have star spots, normal A stars conventionally do not exhibit the same traits. Vega is the first normal A-type star observed to have a weak surface magnetic field (Lignières et al. 2009; Petit et al. 2010), accounting for the rotational modulation. However, the origins of these magnetic fields are uncertain. One possible mechanism is dynamo action, where thin convective layers host a dynamo driven by convective motion, generating a magnetic field. On the other hand, the fields could be "failed fossils," which are generated in the early phases of a star's life and dynamically evolve toward fossil equilibrium (e.g., Braithwaite & Cantiello 2013). These models can be distinguished by the time variation of the magnetic field. Dynamos are intrinsically variable, leading to spot lifetimes similar to the rotation period. In contrast, failed fossil fields would evolve much more slowly.

Vega appears to have a very complex magnetic field. Böhm et al. (2015) found no evidence for rapid variability consistent with dynamo action in their data, suggesting the presence of star spots that last for over 5 days. However, applying Doppler imaging techniques to the same data set, Petit et al. (2017) did find rapidly evolving surface features in combination with stable structures. Cantiello & Braithwaite (2019) suggested that this rapid variation is a sign of dynamo-generated magnetic fields.

While our spectroscopic observations are likely not high-cadence enough to characterize the short-term evolution of corotating structures, they stretch over a span of time much longer than the five consecutive nights used by Böhm et al. (2015) and Petit et al. (2017), which may provide insight regarding the long-term evolution. In Section 3, we described that our RVs can be modeled by either a GP with a quasiperiodic kernel or the combination of Keplerian signals very close to the rotation period. In both cases, the implied timescale for evolution is long: the GP fit suggests an exponential decay timescale of half a year, while the Keplerian model that adequately reproduces the signal seen in our ten-year data set suggests that some surface features may be even more long-lived. It is hard to distinguish between the two models, but both imply the presence of surface features with lifetimes longer than those expected for dynamo fields. Even so, this is not necessarily at odds with the complex and time-varying reconstructed surface map presented by Petit et al. (2017), who identified both variable and stable surface features in their data. We speculate that their stable features may be the ones responsible for our long-lived RV signal, while their rapidly varying features contribute high-frequency noise to our data that we cannot characterize due to our limited observing cadence. We do, however, measure stellar jitter with a magnitude comparable to that of the coherent modulation, which may arise from rapidly varying surface features. If both types of features exist, this may imply magnetic structures influenced both by a failed fossil field and a subsurface dynamo. However, long-term spectropolarimetric observations of Vega are necessary to conclusively determine the behavior of its active regions.

Another interesting feature of Vega's activity is its low photometric variation, particularly given its maximum radial-velocity amplitude of ${14.4}_{-1.8}^{+2.2}$ m s⁻¹. A dark spot inducing a 10 m s⁻¹ variation on a star rotating with $v\sin {i}_{\star }=20$ km s⁻¹ should also induce a photometric signal of ${A}_{\mathrm{RV}}/v\sin {i}_{\star }\,\sim 500\,\mathrm{ppm}$ (Vanderburg et al. 2016). The observed amplitude of only 10 ppm provides independent evidence that dark spots do not dominate the surface. To explore this further, we use SOAP 2.0 (Dumusque et al. 2015) to compare the effects of star spots and plages at different latitudes on a Vega-like star (Figure 8). We vary the size of the active region to reproduce an RV signal with an amplitude of ∼15 m s⁻¹ and compare the expected photometric signals. We find that a dark spot would induce photometric variations ranging from 100 to over 700 ppm. On the other hand, a high-latitude plage can reproduce photometric variations less than 10 ppm, consistent with the observed TESS light curve. This is because plages are only marginally hotter than the rest of the stellar surface, resulting in little flux variation, but are still surrounded by local magnetic fields that suppress convective blueshift and cause regions to appear redshifted, creating an RV signal (Dumusque et al. 2015). Given that Petit et al. (2017) detected many bright and dark regions on Vega's surface, the picture is clearly not as simple as in our simulations, but it suggests that the observed variations are primarily driven by plages.

Zoom In Zoom Out Reset image size

Any planetary signals in our radial velocities may be difficult to detect, since the dominant signal arises from stellar activity and its removal is not straightforward. Nevertheless, we are able to rule out the presence of a candidate signal near 0.53 or 0.56 day that was previously reported by Böhm et al. (2015), and we do identify two signals worthy of further investigation. The first, with a period of ${196.4}_{-1.9}^{+1.6}$ days, would have a semimajor axis of ${0.853}_{-0.022}^{+0.020}$ au and a minimum mass of 0.252 ± 0.066 M_J, falling well below the detection limits of previous direct imaging surveys. The power of the signal monotonically increases with the addition of new data, which is what one would expect for a real orbiting companion. However, it is not formally significant, is located near a small peak in the window function, and its significance is reduced further when the activity signal is modeled. We conclude that there is not good evidence for a planet at this period.

The second interesting signal emerged from the velocity residuals when modeling the activity with multiple Keplerians near the rotation period. We identify a short-period signal with a formal false-alarm probability less than 1%, a best-fit period of 2.43 days, and a semi-amplitude of 6 m s⁻¹, implying a minimum mass of ∼20 M_⊕. A true mass this low would require a polar orbit, but A stars hosting short-period planets do display a wide range of stellar obliquities, and even a highly inclined orbit that is well aligned with the stellar spin (i = 62) would correspond to a planetary mass companion ( ∼0.6 M_J). Continued radial-velocity observations could provide further insight into the presence of a planet orbiting Vega, but more work needs to be done to account for the activity signal. Further spectroscopic observations should be planned carefully, for example, to achieve a cadence that resolves the rotational and orbital periods, allowing for better simultaneous modeling of planet and stellar activity.

If there is a planet with a period of 2.43 days orbiting Vega, there are a few ways that one might confirm its existence. TESS data rule out transits of objects larger than about 2 R_⊕ at this period, so unless it is very dense, this candidate planet does not transit. It has taken 1500 spectra obtained over 10 yr to detect the RV signal at low significance, so it may require a large investment of spectroscopic resources to increase the S/N of the detection even if it is real. A high-cadence campaign spanning many orbital periods (ideally from multiple longitudes to achieve uninterrupted data) might be best suited for mitigating stellar activity and limiting its evolution during observations. Alternatively, it may be possible to directly detect the planet using high-resolution spectroscopy and cross-correlation against a template of elemental or molecular species in the planetary atmospheres (e.g., Snellen et al. 2010). With such a short period around such a hot star, the planet would have an equilibrium temperature of ∼3250 K (assuming a Bond albedo of 0.25), and would be the second hottest known exoplanet after KELT-9b (Gaudi et al. 2017). KELT-9b and several other hot Jupiters have had atmospheric lines successfully detected with high-resolution spectroscopy (e.g., Brogi et al. 2012; Birkby et al. 2013; Yan et al. 2019). Importantly, these detections are not reliant on a transiting geometry. To predict the expected signal, we need to estimate the star-to-planet contrast ratio, and therefore the planetary radius. Nearly all highly irradiated giant planets have inflated radii (some as large as ∼2 R_J), so 1 R_J would be a conservative estimate. As outlined in Birkby (2018), the S/N of the planetary signal is described by

$$\begin{eqnarray}&&{{\rm{S}}/{\rm{N}}}_{{\rm{planet}}}=\left(\displaystyle \frac{{S}_{{\rm{p}}}}{{S}_{\star }}\right){{\rm{S}}/{\rm{N}}}_{{\rm{star}}}\sqrt{{N}_{{\rm{lines}}}},\end{eqnarray} \tag{ 4 }$$

where S_p /S_⋆ is the planet-to-star contrast ratio. Adopting R_p = 1 R_J and assuming blackbody radiation, the planet-to-star contrast ratio for the candidate Vega planet would be 2 × 10⁻⁵. As a reference, assuming 1 R_J for tau Boötis b (the first non-transiting planet with a detection using this technique; Brogi et al. 2012), we would estimate a contrast of 2.6 × 10⁻⁵. The total S/N of our Vega spectra is about 26,000, implying ${{\rm{S}}/{\rm{N}}}_{{\rm{planet}}}\sim 0.5\sqrt{{N}_{{\rm{lines}}}}$. If we can resolve ∼100 lines in the planetary spectrum, it is possible that we could detect such a planet with our TRES data. However, a high confidence detection may have to wait for a data set better suited to this purpose. The SOPHIE data obtained by Böhm et al. (2015) have higher resolving power and total S/N than our TRES data but cover a narrower wavelength range. TRES covers the red optical, where contrasts are not as extreme and where lines from molecules such as TiO, CO, and H₂O may be present in the planetary atmosphere. On the other hand, Fe i has been detected on the daysides of ultra-hot Jupiters like KELT-9b and WASP-33b, so direct detection of the planet is possible even in the blue optical. Ultimately, it may be necessary to move to the near-infrared where contrasts are further reduced and the presence of strong molecular bands of CO and H₂O can boost the S/N.

Using 1524 TRES RVs, we are able to place new detection limits on planets orbiting within 15 au of Vega, a region in which direct imaging surveys have only ruled out brown dwarfs. Assuming no preference for the orientation of the orbital plane—consistent with the observed distribution for short-period planets orbiting hot stars—we can rule out nearly all hot Jupiters. We are sensitive to Saturn-mass objects at 0.1 au, Jupiters at 1 au, 5 M_J at 10 au, and 13 M_J at 15 au. For a planet well aligned with the stellar spin, these masses increase by a factor of about eight and even massive planets would therefore be difficult to detect beyond ∼3 au. Nevertheless, these limits indicate that brown dwarf companions in any orientation are exceedingly unlikely within 3 au. While the majority of planets that have been hypothesized from the observed architecture of Vega's disk would reside beyond 15 au, Raymond & Bonsor (2014) found that systems of planets with masses ranging from ∼1 M_J at 5–10 au to Neptune masses at tens of astronomical units could replenish the hot dust in the inner belt. Our data are not quite sensitive enough to constrain planets this small at these distances, with our detection probabilities falling off steeply below about 2 M_J at 5 au, even for planets drawn from an isotropic inclination distribution.

With TESS data, we are able to place constraints on the size of any transiting planets. While most planets would not transit, requiring both a very high stellar obliquity and a semimajor axis ≲0.2 au, the extremely precise photometry obtained by TESS results in very sensitive detection limits, despite the large stellar radius. We can rule out any transiting hot Jupiters, along with most other planets with radii greater than 3 R_⊕. Further observations could place even more extensive limits on transiting planets while helping us to better understand the active regions on the star. The TESS extended mission will return to the northern hemisphere in its second year of operations, which may be the next opportunity for additional uninterrupted measurements at a similar precision. The 10 minute FFIs used in the extended mission will also improve the time resolution by a factor of three, opening the possibility of further characterization of high-frequency stellar variations.

Current and future space telescopes offer additional opportunities to search for planets around Vega. Meshkat et al. (2018) showed that planned James Webb Space Telescope (JWST) NIRCam GTO observations of Vega (Beichman et al. 2010) will have greater sensitivity than previous surveys, perhaps extending to Saturn-mass planets. Additionally, Mid-Infrared Instrument (MIRI) Guaranteed Time Observations program (GTO) observations are expected to resolve the potential asteroid belt analog (Beichman et al. 2017), providing further insight into the disk structure of Vega and its implications for a planetary system. However, JWST will only be able to search the region beyond 15 (11 au) from Vega for planets; any closer, and the star will saturate the instrument. The Nancy Grace Roman Space Telescope coronagraph instrument (CGI) is more promising for close companions, as it is intended to observe small fields around bright stars. Because Vega is so bright, Roman would not need to observe a bright point-spread function reference star, potentially allowing better CGI stability and improving contrast for point sources. On the other hand, the CGI is optimized for stars with angular diameters less than ∼ 2 mas; with an angular diameter over 3 mas, there may be some light leakage when observing Vega. While it is uncertain how these factors would balance out, if the conventional Roman CGI sensitivity limits apply to the star, the mission could observe Jupiter-sized planets within 02 (close to 1 au; Nemati et al. 2017; Krist et al. 2018). Combining Roman and JWST observations, future space missions promise to place new direct imaging limits on planets throughout the entire Vega system.

As a nearby star, Vega is also an interesting candidate for astrometric detection of companions, though the best precision—i.e., with Gaia—may not be possible. Sahlmann et al. (2016, 2018) described protocols to observe very bright stars with Gaia, but Vega is far beyond the bright limit for standard Gaia processing, and they do not quantify what the uncertainties for such measurements may be. If astrometric precision for Vega were to rival the typical performance (e.g., 35 μas per measurement), we might expect Gaia to be sensitive to Jupiter masses beyond about 1.5 au (Ranalli et al. 2018). However, we do not know the precision with which Gaia can observe a star like Vega, and it is unlikely to be close to the instrumental floor.

Radial velocities remain the most sensitive technique for companions within 1 au for the foreseeable future.