Aldebaran b's temperate past uncovered in planet search data
Will M. Farr Benjamin J. S. Pope Guy R. Davies Thomas S. H. North Timothy R. White Jim W. Barrett Andrea Miglio Mikkel N. Lund Victoria Antoci Mads Fredslund Andersen Frank Grundahl Daniel Huber




0000-0003-1540-8562]Will M. Farr

w.farr@bham.ac.uk

Birmingham Institute for Gravitational Wave Astronomy,
  University of Birmingham, Birmingham, B15 2TT, United Kingdom

School of Physics and Astronomy, University of
  Birmingham, Birmingham, B15 2TT, United Kingdom

Center for Computational Astrophysics, Flatiron Institute, 162 Fifth Avenue, New York NY 10010, United States

0000-0003-2595-9114]Benjamin J. S. Pope
Sydney Institute for Astronomy (SIfA), School of Physics, University of Sydney, NSW 2006, Australia
Center for Cosmology and Particle Physics, Department of Physics, New York University, 726 Broadway, New York, NY 10003, USA
NASA Sagan Fellow

0000-0002-4290-7351]Guy R. Davies

School of Physics and Astronomy, University of
  Birmingham, Birmingham, B15 2TT, United Kingdom
Stellar Astrophysics Centre, Department of Physics
and Astronomy, Aarhus University, Ny Munkegade 120, DK-8000 Aarhus C, Denmark


School of Physics and Astronomy, University of Birmingham, Birmingham, B15 2TT, United Kingdom
Stellar Astrophysics Centre, Department of Physics
and Astronomy, Aarhus University, Ny Munkegade 120, DK-8000 Aarhus C, Denmark

0000-0002-6980-3392]Timothy R. White
Stellar Astrophysics Centre, Department of Physics
and Astronomy, Aarhus University, Ny Munkegade 120, DK-8000 Aarhus C, Denmark


Birmingham Institute for Gravitational Wave Astronomy,
  University of Birmingham, Birmingham, B15 2TT, United Kingdom

School of Physics and Astronomy, University of
  Birmingham, Birmingham, B15 2TT, United Kingdom



School of Physics and Astronomy, University of
  Birmingham, Birmingham, B15 2TT, United Kingdom
Stellar Astrophysics Centre, Department of Physics
and Astronomy, Aarhus University, Ny Munkegade 120, DK-8000 Aarhus C, Denmark



Stellar Astrophysics Centre, Department of Physics
and Astronomy, Aarhus University, Ny Munkegade 120, DK-8000 Aarhus C, Denmark
School of Physics and Astronomy, University of
  Birmingham, Birmingham, B15 2TT, United Kingdom


Stellar Astrophysics Centre, Department of Physics and Astronomy, Aarhus University, Ny Munkegade 120, DK-8000 Aarhus C, Denmark


Stellar Astrophysics Centre, Department of Physics and Astronomy, Aarhus University, Ny Munkegade 120, DK-8000 Aarhus C, Denmark


Stellar Astrophysics Centre, Department of Physics and Astronomy, Aarhus University, Ny Munkegade 120, DK-8000 Aarhus C, Denmark

0000-0001-8832-4488]Daniel Huber
Institute for Astronomy, University of Hawai`i, 2680 Woodlawn Drive, Honolulu, HI 96822, USA
Sydney Institute for Astronomy (SIfA), School of Physics, University of Sydney, NSW 2006, Australia
SETI Institute, 189 Bernardo Avenue, Mountain View, CA 94043, USA
Stellar Astrophysics Centre, Department of Physics and Astronomy, Aarhus University, Ny Munkegade 120, DK-8000 Aarhus C, Denmark


The nearby red giant Aldebaran is known to host a gas giant planetary companion from decades of ground-based spectroscopic radial velocity measurements. Using Gaussian Process-based Continuous Auto-Regressive Moving Average (CARMA) models, we show that these historic data also contain evidence of acoustic oscillations in the star itself, and verify this result with further dedicated ground-based spectroscopy and space-based photometry with the Space Telescope. From the frequency of these oscillations we determine the mass of Aldebaran to be 1.16 ± 0.07 , and note that this implies its planet will have been subject to insolation comparable to the Earth for some of the star's main sequence lifetime. Our approach to sparse, irregularly sampled time series astronomical observations has the potential to unlock asteroseismic measurements for thousands of stars in archival data, and push to lower-mass planets around red giant stars.




§ INTRODUCTION


Aldebaran (α Tauri) is a well-known first-magnitude red giant star, and has long been the subject of astronomical investigations. It was one of the first stars around which an extrasolar planet candidate was identified, by looking for Doppler shifts from the star's reflex motion around the common centre of mass with its companion <cit.>. While the hot Jupiter 51 Peg b <cit.> was the first exoplanet to be recognized as such, before this, <cit.> had noted RV variations in Pollux <cit.>, Arcturus (α Boo, unconfirmed), and Aldebaran. After further investigation by <cit.>, <cit.> have now claimed a firm RV detection of a planetary-mass companion Aldebaran b, with a period of 629.96 ± 0.90 d.

In this paper, we present a re-analysis of these original RV data in which we not only confirm the planetary signal, but detect acoustic oscillations in Aldebaran for the first time. <cit.> noted night-to-night RV variability: this is the noise floor limiting the sensitivity to sub-M_J planets around giants in RV surveys <cit.>. We pull out the asteroseismic signal in this planet hunting noise. We validate this method and its result with new independent RV observations with the Hertzsprung SONG Telescope, and photometry from the Mission. By measuring the frequency of maximum power of these p-mode oscillations, , we asteroseismically determine the mass of Aldebaran to be 1.16 ± 0.07 . This precise stellar mass allows us to calculate that Aldebaran b and any satellites it may have, although they are now likely to be very hot, would have had equilibrium temperatures comparable to that of the Earth when Aldebaran was on the main sequence. It is possible that they may have once been habitable, billions of years ago.

Our new approach to asteroseismic data analysis, based on Continuous Auto-Regressive Moving Average (CARMA) models, can extract exoplanet signals together with measures of from sparse and irregularly-sampled time series. An all-sky survey to find planetary companions and to precisely measure the masses of all nearby red giant stars is feasible with this new approach, and the required data either already exist in large radial velocity exoplanet surveys, or are easy to obtain with ground-based telescopes.



§ ASTEROSEISMOLOGY OF RED GIANTS

Asteroseismology is a powerful tool for the charaterisation of red giant stars <cit.>.  Red giants exhibit oscillations that are excited and damped by stellar convection.  In the power spectrum of either radial velocity or photometric observations, these modes of oscillation form a characteristic pattern of peaks which can be used to infer intrinsic stellar properties <cit.>.  The easiest properties of the pattern to determine are the frequency of maximum amplitude and the so-called large separation <cit.>, which are often referred to as global asteroseismic parameters.


These parameters can be used to estimate the radius R, mass M, and surface gravity g of stars when combined with an estimate of the effective temperature  through scaling relations:


    ( R/R_⊙) ≃   ( ν_ max/ν_ max, ⊙)  ( Δν/Δν_⊙)^-2 ( T_ eff/T_ eff, ⊙)^0.5, 
    ( M/M_⊙) ≃   ( ν_ max/ν_ max, ⊙)^3 ( Δν/Δν_⊙)^-4 ( T_ eff/T_ eff, ⊙)^1.5, 
    ( g/g_⊙) ≃   ( ν_ max/ν_ max, ⊙)  ( T_ eff/T_ eff, ⊙)^0.5.


While the accuracy of these scaling relations is still a matter of ongoing work <cit.>, stellar parameters can be estimated by comparing observables to parameters from models of stellar evolution.


The ability to measure and depends on the length and sampling rate of a data set.
In practice, for the higher luminosity red giants it is more straightforward to measure than .  Typical values for  range from 0.1-20  μ Hz for luminous giants and 20-50  μ Hz for stars near the red clump.  For the evolutionary stages, typically ranges from 0.02-3  μ Hz and 3-7  μ Hz depending on stellar mass and 
<cit.>.  Hence, requires less frequency resolution than to establish a good measurement.





§ TIME-DOMAIN MODELS


Radial velocity variations intrinsic to the star, such as those caused by stellar oscillations, have previously limited the precision with which red giant planets have been studied <cit.>. We aim to model these noise processes and use them to recover asteroseismic information, as well as to improve our estimates of the orbital parameters of the planet. It is easy to observe Aldebaran and similarly-bright stars with ground-based spectroscopic instruments, typically only requiring short exposures that can be obtained even under adverse observing conditions. There is indeed a considerable archive of such observations already, as a legacy of radial velocity (RV) surveys conducted to find exoplanets. In most cases, however, these have not so far been useful for asteroseismology because these RV data are sparsely and irregularly sampled. Because we have to pause observations during the day, during poor weather conditions, or simply when targets of higher priority are being observed, we get time series which may have significant and uneven gaps. This introduces a window-function effect: the power spectrum as constructed for example by a Fourier transform, or a Lomb-Scargle Periodogram <cit.>, is convolved with the Fourier transform of the window function, introducing strong sidelobes adjacent to real frequency peaks and causing crosstalk between adjacent frequency channels. This imposes significant limitations both on the signal-to-noise and frequency resolution of power spectra derived from linear methods such as the Lomb-Scargle periodogram, and in practice makes asteroseismology with these conventional approaches difficult or impossible from the ground for stars with oscillation periods ranging from ∼ 12 h to ∼ a few days.

If we apply nonlinear statistical inference methods this situation can be
improved. <cit.> show that a system of driven, damped harmonic
oscillators such as we encounter in asteroseismology can be statistically
modelled as a Gaussian Process (GP), with a covariance kernel consisting of a sum of
damped sinusoids. The hyperparameters of this process encode features of the
power spectral density distributions and are insensitive to the window function. GPs with quasi-periodic kernels have previously been used to re-analyse RV data for main sequence stars <cit.>, where the noise is from stellar activity, but have not previously been applied to red giant asteroseismology.
Unfortunately, these methods are impractical for long time series, as the computational
cost of evaluating the standard GP likelihood function scales as
𝒪(n^3).

Fortunately, the damped and driven harmonic oscillator GP can be written as the
solution of a class of stochastic ordinary differential equations, the
Continuous Auto-Regressive Moving Average (CARMA) models, whose likelihoods can
be evaluated in linear time. For problems such as this, we therefore have access
to these powerful computational tools for inferring their power spectral
densities.  Our treatment of CARMA models is described more fully in Appendix
<ref>.   <cit.> also give more details on CARMA
processes in an astronomical context, and show how a state-space model of the
process y(t) can be tracked through a time-series of uncertain observations
y_k = y( t_k ) + ϵ_k, with Gaussian noise ϵ_k of
zero mean and known (heteroskedastic) variance, using a Kalman filter to produce
a likelihood function p( y_k |σ, r_i, b_j ) depending on
the amplitude of the stochastic forcing, σ, the (possibly complex)
eigenfrequencies of the ODE, r_i, and parameters describing the power spectrum
of the stochastic forcing, b_j in 𝒪( n ) computational
time. <cit.> has demonstrated that the same model can also be implemented via a novel matrix factorization, and demonstrated its asteroseismic potential modelling the light curve of KIC 11615890.

Here we have adopted the <cit.> approach.  The 
package[<https://github.com/farr/CARMA.jl>] implements Kalman
filters that can compute the likelihood for a set of observations parameterised
by either the r_i and b_j or the RMS amplitudes A_i of each eigenmode of
the ODE and the corresponding roots r_i. This likelihood is computed for the
residuals of a deterministic Keplerian RV mean model for the planetary
companion.  Because we have a set of observations taken with different
instruments at different sites (see Section <ref>), we include
RV offset (i.e. mean value) parameters that are instrument- and site-dependent
and also an instrument- and site-dependent uncertainty scaling parameter.  Thus,
the ith RV measured at site and intstrument combination j, y_ij, is
assumed to be


    y_ij = v_r(t_i) + μ_j + y( t_i ) + ν_j
    ϵ_ij,


where t_i is the time of the measurement, v_r is the planetary RV signal,
μ_j is the mean RV offset for site and instrument combination j, y(t) is
the CARMA process representing stellar activity, ν_j is the uncertainty
scaling parameter for site and instrument combination j, and ϵ_ij
is an independent Gaussian random variable with mean zero and standard deviation
equal to the quoted uncertainty of the observation.  We impose broad priors (see
Table <ref>) on the parameters of the planetary RV signal, the
CARMA process, the μ_j offsets, and the ν_j scale factors.  We treat the
frequency of maximum asteroseismic power, ν_max, as the frequency
of an appropriate eigenmode in the CARMA process representing the stellar
activity; this is equivalent to fitting a Lorentzian profile to a set of
unresolved asteroseismic modes.  We used the nested sampling algorithm of the
 package[<https://github.com/farr/Ensemble.jl>;
this package implements several stochastic sampling algorithms based on the
“stretch move” proposal used in  <cit.> and
introduced in <cit.>] to calculate the marginal likelihood
(evidences) and draw samples from the posterior distribution over the parameters
of our various models.



§ ALDEBARAN

Aldebaran is a red giant star with spectral type K5, one of the nearest such stars at a distance of only 19.96 ± 0.38 pc as determined by <cit.>. Its position near the Ecliptic permits the determination of its angular diameter by lunar occultations and by interferometry <cit.>. These tight constraints are valuable in breaking degeneracies in stellar modelling and make this an ideal star for asteroseismic characterization (Appendix <ref>).



 §.§ Archival Observations



<cit.> reported on RV observations of Aldebaran from the coudé
échelle spectrograph of the Thüringer Landessternwarte Tautenburg (TLS),
the Tull Spectrograph of the McDonald 2.7 m telescope, and the Bohyunsan
Observatory Échelle Spectrograph (BOES) spectrograph of the Bohyunsan
Optical Astronomy Observatory (BOAO) spanning the period from 2000.01 to
2013.92.  Combined with earlier observations <cit.>, we construct a
full time series spanning more than two decades from 1980.80 to 2013.71
(Figure <ref>) and fit this with a Keplerian model and a CARMA process.
We find evidence for a low-quality (Q =) oscillatory mode with
ν_max =  μHz (here and throughout
we quote the posterior median and the range from the 0.16 to the 0.84 posterior
quantile).  This is consistent with the presence of a number of un-resolved
asteroseismic modes in the RV data.  The low Q factor in this case refers to the
inverse fractional width of the band over which the modes are excited, not the
quality of any individual mode.  Our posterior for  from this data set
is shown in Figure <ref>.







 §.§ SONG Observations
    
 In order to confirm the oscillations detected in
the archival data, we used the Hertzsprung SONG telescope
<cit.> to conduct high time cadence follow-up observations.
These were carried out in the highest resolution mode (R≈110000) with an
integration time of 30 s, between 2016 September 27 and 2017 January 12. During
this time we attempted to obtain at least one visit per available night, for a
total of 254 epochs over the campaign. The spectra were extracted using the SONG
pipeline <cit.>. SONG employs an iodine cell for
precise wavelength calibration and the velocities were determined using the
iSONG software following the same procedures as described in
<cit.>. The typical velocity precision per visit is
2.5m s^-1 allowing us to easily detect the oscillations which display
a peak-to-peak amplitude of ≈  170   m   s^-1.
Figure <ref> displays the radial velocities obtained—oscillations
are easily visible as well as a long-term trend, due to the planetary companion.
From this data set alone we also find evidence for a low-quality (Q =) oscillatory mode with ν_max =  μHz which is consistent with the inference from the archival data.
This posterior is shown in Figure <ref>.










 §.§ K2 Observations


In order to verify the results of the novel analysis presented above, we sought to obtain an independent detection of the oscillations of Aldebaran and compare the frequencies determined with the two methods.
Aldebaran was observed with <cit.>, a two-wheeled revival of the Mission <cit.>, under Guest Observer Program 130471 in Campaign 13, from 2017 March 8 to 2017 May 27. Aldebaran and similar high-luminosity K giants are well known to show photometric variability due to solar-like oscillations <cit.>, and hence the continuous, high-precision K2 light curves allow an independent confirmation of our CARMA results.

As Aldebaran is extremely bright, it saturates the detector and it is therefore not possible to use standard photometry pipelines to extract a K2 lightcurve. We instead use halo photometry <cit.>, whereby unsaturated pixels from the outer part of the large and complicated halo of scattered light around bright stars are used to reconstruct a light curve as a weighted linear combination of individual pixel-level time series, as described in Appendix <ref>.

Figure <ref> shows the light curve and power spectrum of the K2 observations of Aldebaran. We detect clear variability on an average timescale of ∼ 5 days, consistent with the expected timescale for solar-like oscillations for a high-luminosity red giant. To measure global asteroseismic parameters we model the background variability in the Fourier domain using the methodology described in <cit.>, yielding a frequency of maximum power of = 2.2 ± 0.25 and amplitude per radial mode of 1850 ± 500 ppm. Due to the limited frequency resolution of the 70-day K2 time series we were unable to measure the large frequency separation , which is expected to be ≈ 0.5.

Using the k2ps planet-search code <cit.> to examine this light curve, we search for transits across a wide range of periods, and find no evidence either of short-period planetary transits or an eclipsing stellar companion.





 §.§ Planetary and Stellar Parameters



By jointly fitting a Keplerian for the planet with a CARMA model for the stellar
oscillations to the combined archival and SONG data sets (Figure <ref>),
we obtain a precise estimate of both and the planetary orbital
parameters. The planetary orbital parameters we obtain are similar to
<cit.>, but with larger uncertainty that is likely due to our
more-flexible and correlated model for the stellar component of the RV signal:
period P =  d,  eccentricity e =; and radial velocity semi-amplitude K =
 m s^-1.

The stellar oscillations are best fit with a single, low-quality (Q =) oscillatory mode with = 2.24^+ 0.09_- 0.08
.  We find no improvement in the marginal likelihood (evidence) or other
model-selection information criteria <cit.> from models with
additional oscillatory modes, so we conclude that the large spacing Δν
is not observable in this data set.

We have used the additional constraint of determined from this analysis
to update the stellar properties of Aldebaran.  The details of the stellar
modelling are presented in Appendix <ref>.  We have considered the
impact of our additional constraint and run our analysis for different
sets of spectroscopic estimates.  For final values we adopt the
<cit.> spectroscopic solution as being `middle of the road'
estimates.  For the analysis without constraint on we find estimates of
the stellar properties as M = 1.27^+0.24_-0.20 M_⊙ and
age 4.9^+3.6_-2.0  Gyr.  With the inclusion of we find M =
1.16 ± 0.07 M_⊙ and age 6.4^+1.4_-1.1 
Gyr.  It is clear that the addition of provides substantially more
precise estimates of mass and age. With this stellar mass, the radial velocity
signal translates to an inclination-dependent planetary mass of M sin i =
5.8 ± 0.7 M_J <cit.>, a massive giant planet; for
sufficiently high inclinations (≳ 45^∘), the mass of Aldebaran b may
exceed 13 M_J, making it a possible low-mass brown dwarf candidate. A recent
study by <cit.> has cast doubt on the validity of the
detection of the massive planet γ Draconis b, and by extension similar
planets around similar K giants such as Aldebaran. We remain convinced that the
γ Draconis problem is not an issue in this case, and Aldebaran b is a
bona fide planet, for reasons outlined in Appendix <ref>.

With the determination of the stellar mass it is possible to infer the parameters Aldebaran had while it was on the main sequence. We conduct a Monte Carlo simulation, drawing masses randomly from a Gaussian distribution M ∼𝒩(1.16,0.07) and metallicities [Fe/H] ∼𝒩(-0.15,0.2) <cit.> to predict the luminosity of the main-sequence progenitor, and the semi-major axis of the planet's orbit.
Using the Mesa Isochrones and Stellar Tracks <cit.> models, we compute the evolution of the luminosity L of the star along the main sequence; L evolves from 2.0 ± 0.7 at 0.5 Gyr to 3.5 ± 2.3 at 4.5 Gyr.
We note that this implies that the planet at a semi-major axis of 1.50 ±  0.03 AU would have been subject to an insolation comparable to that of the Earth, evolving from 0.5 Gyr to 4.5 Gyr from 0.86 ±  0.27 to 1.55  ± 1.0 times that of the Earth today (Figure <ref>). Subject as well as this to the great uncertainties of the planet's orbital evolution and albedo, Aldebaran b and any of its moons (or its S-type planets if it is a brown dwarf) may well have hosted temperate environments for some of their history, now long-since destroyed by their star's evolution away from the main sequence.






§ CONCLUSIONS


Using sophisticated time-domain CARMA models for the stellar activity in
Aldebaran, we have confirmed the previously-suspected planet and additionally
detected acoustic oscillations that permit a mass determination with a precision
of ∼ 6%. We have confirmed both of these results with new ground-based
observations, and space-based photometry with .

From this pilot study with Aldebaran, we have shown that with limited quantities of data of limited quality, we can nevertheless do asteroseismology from the ground and obtain very precise estimates of stellar parameters. Furthermore, improvements in the algorithm may also permit the detection of not only but also from similar measurements, permitting more sophisticated asteroseismic analysis. Sufficient radial velocity data either already exist, or are trivial to obtain, in order to do this for essentially all bright giants; while for asteroseismology of solar-like stars, these new methods will allow significantly relaxed observing requirements. One rich archive is the series of observations with the Hamilton Échelle Spectrograph at Lick Observatory, studying 373 G and K giants <cit.>. The method may also be useful in detecting rotational or other periodic variability in the LCES HIRES/Keck Precision Radial Velocity Exoplanet Survey <cit.> of 1,624 FGKM dwarfs, although the time sampling is by design too coarse to detect p-mode oscillations in these stars. Furthermore, CARMA models and related methods will enhance deep, all-sky, sparse photometric surveys: an immediate future test of this will be from ; while only 58 epochs of photometry are available for Aldebaran, at a sampling we find to be insufficient for our purposes, stars at higher latitudes may often have 150–200 epochs and more even sampling <cit.>, and the 14 bright K giants in noted by <cit.> are an ideal first test case. This can naturally be extended to Gaia in space <cit.>, or LSST from the ground <cit.>, from which many thousands of new asteroseismic determinations will be possible. In future RV surveys, it may be possible to beat the <cit.> ∼ 30 m/s RV precision limit for red giants by taking many closely spaced observations and modelling-out the effects of stellar oscillations: this will allow us to dig deeper into this intrinsic stellar noise to detect less-massive planets around these stars.



We would like to thank Tim Bedding, David Hogg and Dan Foreman-Mackey for their helpful discussions.

We would like to commend the authors of <cit.> for making their data
publicly available on SIMBAD.

The authors would like to acknowledge the Gadigal people of the Eora Nation on whose ancestral lands the University of Sydney is built.

BP is grateful for funding from the Breakthrough Prize Foundation. This work was performed in part under contract with the Jet Propulsion Laboratory (JPL) funded by NASA through the Sagan Fellowship Program executed by the NASA Exoplanet Science Institute.
TRW acknowledges the support of the Villum Foundation (research grant 10118).
D.H. acknowledges support by the Australian Research Council's Discovery Projects funding scheme (project number DE140101364) and support by the National Aeronautics and Space Administration under Grants NNX14AB92G and NNX16AH45G issued through the Kepler Participating Scientist Program and the K2 Guest Observer Program.

Funding for the Stellar Astrophysics Centre is provided by the Danish National Research Foundation.
This research has made use of the SIMBAD database, operated at CDS, Strasbourg, France, and NASA's Astrophysics Data System. The data presented in this paper were obtained from the Mikulski Archive for Space Telescopes (MAST). STScI is operated by the Association of Universities for Research in Astronomy, Inc., under NASA contract NAS5-26555. Support for MAST for non-HST data is provided by the NASA Office of Space Science via grant NNX13AC07G and by other grants and contracts.

This work has made use of observations made with the Hertzsprung SONG telescope operated on the Spanish Observatorio
  del Teide on the island of Tenerife by Aarhus and Copenhagen Universities and by the Instituto de
  Astrofìsica de Canarias.

The Villum Foundation, the Danish Council for Independent Research | Natural Sciences and
 the Carlsberg Foundation are thanked for their support on building the SONG prototype on
 Tenerife.
 We thank Jens Jessen-Hansen, Rasmus Handberg and René Tronsgaard Rasmussen for their contribution
 to the development of the SONG data reduction pipeline and archive system.

 WMF, GRD, TN, JWB, AM, and MNL are supported in part by the STFC.

IPython <cit.>, SciPy <cit.>, matplotlib <cit.>, Astropy <cit.>, NumPy <cit.>, isochrones <cit.>, Julia <cit.>.

The data, code, and  source used to produce this document can be found
at <https://github.com/farr/Aldebaran> under an open-source “MIT” license.









§ CARMA MODELS



A CARMA process corresponds to the solution, y(t), of the stochastic ODE


    ∏_i=1^p[ t - r_i ] y(t) =
      ∏_j = 1^q[ t - b_i ] η(t),


where p, q ∈ℤ, p > q ≥ 0 define the order of the process;
r_i ∈ℂ, r_i < 0 are the roots of the characteristic equation
of the ODE; b_j ∈ℂ, b_j < 0 are corresponding roots of the
inverse process; and η(t) is a white-noise GP with
⟨η(t) ⟩ = 0 and ⟨η(t) η(t')
⟩ = σ^2 δ(t - t' ).  The constraint that r_i < 0 and b_j < 0 ensures that the linear operator defining the process
and its inverse are invertable.  Under the assumption that the solution, y(t),
is real, roots r_i and b_j must either be real or occur in complex conjugate
pairs.  The power spectrum of the process is


    P_y(f) = σ^2 | ∏_j=1^q [ 2 π i f - b_j ]/∏_i=1^p [ 2 π i f - r_i ]|^2.


The power spectrum is a rational function of frequency with poles at 2π i f =
r_i and zeros at 2π i f = b_i; because of the invertable constraints, these
poles and zeros all occur when f ≠ 0.  The autocorrelation function
corresponding to the power spectrum in Eq. (<ref>) is


    ρ_y (τ) = ⟨ y( t ) y ( t + τ) ⟩ = ∑_i=1^p A_i^2 e^r_i τ,


where the RMS amplitudes of the p different modes, A_i, are functions of
σ and the values of the roots r_i and b_j.  If q = p-1, then the
A_i are independent, while if q < p-1 then there are interdependencies among
the A_i.  From either the power spectrum in Eq. (<ref>) or the
autocorrelation function in Eq. (<ref>), it is apparent that a
real root r_i represents an exponentially-decaying mode with r_i the rate
constant, while a complex conjugate pair of r_i represent an oscillatory mode
with decay rate r_i and angular frequency | r_i|.
Oscillatory modes can also be described by their frequency f and quality
factor Q via r = f/(2Q) ± 2π i f.  Decaying modes generate a damped
random walk, also known as a Ornstein-Uhlenbeck process.  More discussion of CARMA processes, including formulae for A_i in terms of σ, r_i, and b_j, can be found in <cit.> and references therein.

Our preferred model for the data uses one decaying mode and one oscillatory mode
in the CARMA process to capture the un-resolved superposition of many
asteroseismic modes (see Section <ref>); we have fit
models with more than one mode of each type, but find no improvement in
model-selection criteria from the expanded parameter space.  The priors we
impose in all cases on the parameters of the CARMA model, planetary RV signal,
and instrumental and site effects are given in Table <ref>.  A
draw from the posterior for our canonical model in shown in data space in Figure
<ref>.



The Kalman filter implementation of the CARMA likelihood is not unique.
<cit.> show how a dimensional expansion can be used to
reduce the standard GP covariance matrix for this same process,  with
autocorrelation given by Eq. (<ref>) under observations at times
t_k, to a banded form that can then be Cholesky-decomposed in the standard GP
likelihood function in 𝒪( n ) time. These
`celerite' models are formally equivalent to a Kalman CARMA model with q
= p-1 and differ only in their computational implementation.

lll
3
 The parameters in our model Eq. (<ref>) and their priors.  Site and instrument combinations are indexed by j.  Modes in the stochastic process are indexed by k ∈ [1,p] and are either real (parameterised by a rate parameter r_k) or complex (parameterised by a frequency f_k and quality factor Q_k).
Parameter     Prior     Description
Site/Instrument RV Offset and Uncertainty Scaling
μ_j     N[ ⟨ y_j⟩, 10 σ_y_j] ( μ _j )a     The site/instrument RV offset. 

ν_j     1/ν_j,   ν_j ∈ [1/2, 2]     The site/instrument uncertainty scale factor. 

Keplerian Parameters
K     1/K,    K ∈[ min_j σ_y_j/100, max_j 10 σ_y_j]     RV semi-amplitude. 

P     1/P,    P ∈ [600, 700]  d     Keplerian period 

e     2e,   0 ≤ e < 1b     Keplerian eccentricity 

ω     const,    0 ≤ω < 2π     Longitude of pericentre 

χ     const,    0 ≤χ < 1     Pericentre passage is at t = χ P

CARMA Parameters
A_k     1/A_k,    A_k ∈[min_j σ_y_j/100, max_j 10 σ_y_j]     RMS amplitude of stochastic mode 

r_k     1/r_k,    r_k ∈[ 1/ (2T), 0.05  μHz],  r_k = 0c     The damping rates for real modes 

f_k     1/f_k,    f_k ∈[ 0.1, 10 ] μHz     The frequency of oscillatory modes 

Q_k     1/Q_k,    Q_k ∈[ 1, 1000 ]     The quality factor of oscillatory modes 

a⟨ y_j ⟩ and σ_y_j are the mean and standard deviation of the RV time series from site and instrument combination j.
bWe actually sample in and impose a (flat) prior on x ≡ e cosω and y ≡ e sinω.
cHere T is the time span of the measurements over all sites.




§ HALO PHOTOMETRY



The Space Telescope <cit.>
suffered a critical reaction wheel failure in May 2013, which made it impossible to maintain a stable pointing and therefore continue its nominal mission. It was revived as <cit.>, balanced by orienting perpendicular to the Sun. This requires that observes fields in the Ecliptic in ∼80 d Campaigns; Aldebaran was observed in Campaign 13.

The detector saturates for stars brighter than the ∼11^th magnitude. Nevertheless, the excess flux deposited in a saturated pixel spills conservatively up and down the pixel column, such that it is possible to sum this `bleed column' for bright stars and still obtain precise photometry, such as was done for the brightest star in the nominal mission, θ Cyg <cit.>. There are two main reasons why this is not possible for all bright stars in general. First, because the on-board data storage and downlink bandwidth from are limited, it is often not desirable to store and download the large number of pixels that are required for such bright stars. Second, if the bleed column for a sufficiently bright star reaches the edge of the chip, flux spills over and is not conserved, imposing a hard brightness limit that depends on the distance to the detector edge.
Collateral `smear' data, which are collected to help calibrate the photometric bias from stars sharing the same column as a target, can be used to reconstruct light curves for un-downloaded bright stars and thereby avoid bandwidth constraints <cit.>, but these data are still rendered unusable if the bleed column falls off the edge of the chip and contaminates the smear rows.

Bright stars have a wide, complicated, position-dependent point spread function (PSF) arising from diffraction and scattering from secondary and higher-order reflections inside the instrument, with the result that they may contaminate thousands of nearby pixels with significant flux. We can therefore use this `halo' of unsaturated pixels for photometry. The brightness of this halo varies in the same way as that of the primary star, and we therefore obtain data in a region of 20 pixel radius around the mean position of Aldebaran, and discard saturated pixels.
In this paper we proceed as in <cit.>, in which the method was demonstrated on the seven bright Pleiades, with only minor changes.

The flux f_i at each cadence i is constructed as a weighted sum of pixel values p_ij:


    f_i = ∑_j=1^M w_j p_ij.


We choose the weights w_j such that they lie between 0 and 1, add to unity, and minimize the Total Variation (TV) of the weighted light curve. In the continuous case, n-th order TV is defined as the integral of the absolute value of the n-th derivative of a function; in the discrete case, replacing the derivative with finite differences, first-order TV becomes


    TV = ∑_i=1^N |f_i - f_i-1|∑_j=1^Nf_j


and likewise second-order TV the equivalent expression in second-order finite differences. The efficacy of this method was recently confirmed by <cit.>, comparing BRITE-Constellation observations of Atlas to halo photometry and finding excellent agreement in the frequency and amplitude of the reported oscillations.

In an improvement since <cit.>, we use the autograd library <cit.> to calculate analytic derivatives for the TV objective function, which reduces the computational time for a single halo light curve on a commercial laptop from tens of minutes to a few seconds.

As a final step to reduce residual uncorrected systematics, we apply the k2sc <cit.> GP-based systematics correction code to the initial halo light curve, but the effect of this in the present case for Aldebaran is minimal.
There is somewhat higher than usual residual noise at harmonics of 46 μHz (4 d^-1, the satellite thruster firing frequency), but this is nevertheless very small in comparison to the signal from Aldebaran, and may be ascribed to the large fraction of the pixel mask occupied by the bleed column from this extremely bright star.



§ STELLAR MODELLING



 §.§ Stellar Models

We used our determination of and several combinations of the asteroseismic and spectroscopic parameters, along with luminosity, to estimate the fundamental stellar parameters, via fitting to stellar models. We used MESA models <cit.> in conjunction with the Bayesian code PARAM <cit.>.  A summary of our selected “benchmark” options is as follows;



  * Heavy element partitioning from <cit.>.

  * OPAL equation of state <cit.> along with OPAL opacities <cit.>, with complementary values at low temperatures from <cit.>.

  * Nuclear reaction rates from NACRE <cit.>.

  * The atmosphere model according to <cit.>.

  * The mixing length theory to describe convection (we adopt a solar-calibrated parameter α_MLT =1.9657).

  * Convective overshooting on the main sequence is set to α_ov=0.2H_p, with H_p the pressure scale height at the border of the convective core. Overshooting was applied according to the <cit.> step function
scheme.

  * No rotational mixing or diffusion is included.


We do not need to correct for the line-of-sight Doppler shift at the frequency precision available in our data <cit.>.


 §.§ Additional modelling inputs

In addition to the asteroseismic parameters, spectroscopically-determined temperature and metallicity values are needed. There exist multiple literature values for Aldebaran. We chose to compare a range of literature values to investigate what uncertainty these systematically-differing models produce in inferred stellar properties.

To ensure the values are self-consistent, when a literature value was chosen for temperature, we took the stellar metallicity from the same source i.e. matched pairs of temperature and metallicity. The final constraint is the stellar luminosity, which may be estimated as follows (e.g. see ):

    log_10L/L_⊙ = 4.0+
    0.4 M_bol,⊙ -2.0 log_10π [mas] -0.4(V-A_V + BC(V)).

The solar bolometric magnitude M_bol,⊙=4.73 is taken from <cit.>, from which we also take the polynomial expression for the bolometric correction BC(V). We assume extinction A_V to be zero.

The final constraint available for Aldebaran is the angular diameter of the star as measured by long baseline interferometry and lunar occultations <cit.>, combined with the parallax of 19.96 ± 0.38 pc to produce a physical radius constraint of R_int=44.2±0.9R_⊙.




As Table <ref> shows, the spectroscopic parameters of Aldebaran are somewhat unclear, particularly logg and [FeH], which may have an impact on the recovered stellar properties when fitting to models. To explore what impact each parameter is having on the final stellar properties, multiple PARAM runs were performed, using different constraints. Two constraints potentially in tension were ν_max and the spectroscopically-determined logg. ν_max has been shown to scale with the stellar logg <cit.>,

    ν_max/ν_max,⊙=logg/logg,⊙(T_eff/T_eff,⊙)^-1/2.


Using Eq <ref> with the values in Table <ref> predicts ν_max in the range 0.5-6μHz. Reversing the equation to produce a predicted logg from the observed ν_max,obs=2.23±0.1μHz results in a predicted logg∼1.2 dex, using an assumed temperature of 3900K. The solar calibration values used here are log g_⊙ / ( m s^-1)=2.44  dex, ν_max,⊙=3150 μHz and T_eff,⊙=5777K.

Table <ref> shows the results, for all modelling variations, both different inputs and different constraints. It shows that results with the addition of ν_max as a constraint exhibit in general smaller uncertainties, with or without the addition of logg as a constraint.

Recovering the mass without the use of asteroseismic constraints produces considerable scatter on the results (0.96-1.5M_⊙), whilst the use of asteroseismology brings the mass estimates into closer agreement with one another, with the exception of the very low metallicity solution of <cit.>. Any systematic offset between asteroseismic masses and spectroscopic masses is sensitive to chosen reference mass, in agreement with <cit.>.

We have also investigated running PARAM and relaxing one or more constraints: between runs with and without the luminosity constraint, the average absolute mass offset between the two sets of results was <0.04M_⊙ if we also ignore the ν_max constraint, and <0.02M_⊙ including the ν_max.






§ THE Γ DRACONIS PROBLEM IS NOT A PROBLEM



γ Draconis (also known as Eltanin, or γ Dra) is a second-magnitude K5III giant which had been observed by <cit.> as part of the same campaign that led to the discovery of Aldebaran b and Pollux b. <cit.> have recently shown that the 702 d period RV variations of the putative Eltanin b disappeared from 2011-2013, returning in 2014 with a different amplitude and phase. They ascribe this to a previously-unidentified kind of stellar variability, and warn that “Given that the periods found in α Tau [Aldebaran] are comparable to those in γ Dra and both stars are evolved with large radii, a closer scrutiny of the RV variability of α Tau is warranted.”

For several reasons we are not convinced that this poses an issue for Aldebaran or for K giant RV planets more generally. <cit.> suggest that the new type of stellar variability may be oscillatory convective modes, but both Aldebaran and Eltanin have much lower luminosities and longer periods than would seem to be allowed by the period-luminosity relation predicted for these otherwise unobserved modes <cit.>. If these are identified with the long secondary periods (LSPs) observed in some bright red giants (L ≈ 1000), then Aldebaran and Eltanin are both too faint, and lack the mid-IR excess typical of LSP stars <cit.>. It would also be surprising if the shape of the RV curve could reproduce the harmonic structure of an eccentric Keplerian such as in Aldebaran b (e =), let alone the much higher eccentricities observed in other giants <cit.>.

It would be a cruel conspiracy of nature if red giants support a type of oscillation which is common and closely resembles a planetary signal. We believe this cannot be the case: the populations of planets around subgiants and giants evolved from intermediate mass stars are similar <cit.>, as expected if subgiants evolve into giants and retain their planetary systems, and with very different stellar structures subgiants and giants are unlikely to share modes of long-period pulsation. The similarity of distributions of systems hosted by subgiant and giants could not be reproduced if a large fraction of the giants' planets were false positives. Moreover giant planets are expected to be common around intermediate-mass stars <cit.>, and some are definitely known to be bona fide planets because they transit their star <cit.>.

We therefore believe that either Eltanin b is not a false positive, or if it is, that it is not a common type of false positive, and is unlikely to affect our certainty that Aldebaran b is real. <cit.> offer an alternative to the pulsation hypothesis, namely beating between the stellar rotation and the planetary signal; this does not seem implausible.




