Three years of Sun-as-a-star radial-velocity observations on the approach to solar minimum

TL;DR

This study presents three years of high-precision solar radial-velocity measurements using a dedicated telescope and spectrometer, analyzing solar activity effects and noise sources to improve exoplanet detection methods.

Contribution

It introduces a Bayesian data-quality approach and provides detailed analysis of solar activity impacts on radial velocities, advancing techniques for detecting Earth-like planets.

Findings

Photon-noise limited precision better than 0.43 m/s per 5-minute observation

Granulation noise dominates within a day with 0.4 m/s amplitude

Active regions cause 5 m/s radial-velocity signals and line asymmetry changes

Abstract

The time-variable velocity fields of solar-type stars limit the precision of radial-velocity determinations of their planets' masses, obstructing detection of Earth twins. Since 2015 July we have been monitoring disc-integrated sunlight in daytime using a purpose-built solar telescope and fibre feed to the HARPS-N stellar radial-velocity spectrometer. We present and analyse the solar radial-velocity measurements and cross-correlation function (CCF) parameters obtained in the first 3 years of observation, interpreting them in the context of spatially-resolved solar observations. We describe a Bayesian mixture-model approach to automated data-quality monitoring. We provide dynamical and daily differential-extinction corrections to place the radial velocities in the heliocentric reference frame, and the CCF shape parameters in the sidereal frame. We achieve a photon-noise limited radial-velocity precision better than 0.43 m s$^{-1}$ per 5-minute observation. The day-to-day precision is limited by zero-point calibration uncertainty with an RMS scatter of about 0.4 m s$^{-1}$. We find significant signals from granulation and solar activity. Within a day, granulation noise dominates, with an amplitude of about 0.4 m s$^{-1}$ and an autocorrelation half-life of 15 minutes. On longer timescales, activity dominates. Sunspot groups broaden the CCF as they cross the solar disc. Facular regions temporarily reduce the intrinsic asymmetry of the CCF. The radial-velocity increase that accompanies an active-region passage has a typical amplitude of 5 m s$^{-1}$ and is correlated with the line asymmetry, but leads it by 3 days. Spectral line-shape variability thus shows promise as a proxy for recovering the true radial velocity.

Tables

Table 1. Table 1: Bayesian prior probability distributions for the parameters of the Gaussian mixture model used to determine daily extinction coefficients and data quality.

Variable	Parameter	Prior
${\theta }_1$	Fiducial magnitude $m_{60}\left(\hat{x}\right)$ at airmass $x=\hat{x}$	$𝒰[-0.4,0.4]$
${\theta }_2$	Extinction coefficient $k_{60}$	$𝒰[0.075,0.4]$
${\theta }_3$	Transparency fluctuation amplitude ${\sigma }_{\mathrm{jit}}$	$𝒰[0.001,0.08]$
${\theta }_4$	Outlier zero point offset $b_{\mathrm{bg}}$ relative to $\hat{y}$	$𝒰[0.0,0.8]$
${\theta }_5$	Outlier standard deviation ${\sigma }_{\mathrm{bg}}$	$𝒰[0.1,0.7]$
$Q$	Fraction of data in foreground population	$𝒰[0.0,1.0]$

Table 2. Table 2: Dates and daily extinction and mixture-model parameters for each day of observation.

JD-2400000.0	$k_{60}$	$m_{60}\left(x=0\right)$	$b_{\mathrm{bg}}$	${\sigma }_{\mathrm{bg}}$	${\sigma }_{\mathrm{jit}}$	Q
Day	Mag/airmass	Mag	Mag	Mag
57222.1914	0.1415	-13.280	0.509	0.394	0.0030	0.931
57223.9179	0.1213	-13.369	0.456	0.391	0.0085	0.884
57233.0057	0.0965	-13.323	0.734	0.386	0.0035	0.976
57234.0432	0.0891	-13.309	0.462	0.418	0.0042	0.991
57235.0145	0.0995	-13.326	0.465	0.425	0.0044	0.993
57235.9870	0.1074	-13.335	0.471	0.391	0.0059	0.991
57237.0092	0.1089	-13.346	0.489	0.422	0.0058	0.993
57237.9919	0.0902	-13.262	0.506	0.414	0.0103	0.991
57239.0619	0.1351	-13.339	0.479	0.383	0.0016	0.979
57239.9333	0.1090	-13.321	0.623	0.354	0.0065	0.965
57242.0245	0.2982	-13.463	0.432	0.417	0.0505	0.982
57242.8892	0.3109	-13.355	0.335	0.564	0.0561	0.330
57244.9912	0.2314	-13.353	0.181	0.290	0.0058	0.670
57252.0397	0.0963	-13.320	0.634	0.534	0.0079	0.976
57252.9576	0.0841	-13.297	0.437	0.547	0.0079	0.742
57253.9311	0.0829	-13.297	0.882	0.305	0.0071	0.944
.	.	.	.	.	.	.
.	.	.	.	.	.	.
.	.	.	.	.	.	.

Table 3. Table 3: Radial velocities and CCF profile parameters (provided online-only).

BJD	JD	RVbary	RVhel	RVfinal	RVerr	FWHMobs	FWHMsid	Contrast	BISobs	H	Airmass	$m_{60}$	$\sigma \left(m_{60}\right)$	Dec	Diam	PA	Qualflag
57222.1780	57222.1838	0.1122	0.1040	0.1043	0.0004	7.0064	7.0485	45.990	-0.0186	0.820	1.366	-13.0887	0.0026	0.36708	0.009155	0.0886	0.9992
57222.1817	57222.1876	0.1123	0.1041	0.1044	0.0004	7.0081	7.0502	45.987	-0.0193	0.843	1.393	-13.0840	0.0026	0.36706	0.009155	0.0886	0.9993
57222.1855	57222.1913	0.1121	0.1040	0.1043	0.0004	7.0066	7.0487	45.992	-0.0184	0.867	1.422	-13.0750	0.0026	0.36705	0.009155	0.0886	0.9990
57222.1892	57222.1951	0.1115	0.1034	0.1037	0.0004	7.0066	7.0488	45.999	-0.0174	0.890	1.453	-13.0713	0.0026	0.36704	0.009155	0.0887	0.9991
57222.1930	57222.1989	0.1119	0.1038	0.1041	0.0004	7.0083	7.0504	45.989	-0.0194	0.914	1.485	-13.0745	0.0026	0.36703	0.009155	0.0887	0.9985
.	.	.	.	.	.	.	.	.	.	.	.	.	.	.	.	.	.
.	.	.	.	.	.	.	.	.	.	.	.	.	.	.	.	.	.
.	.	.	.	.	.	.	.	.	.	.	.	.	.	.	.	.	.

Equations

$$W=\frac{CF}{2}\sqrt{\frac{\pi }{\ln 2}}.$$

$$y_i\equiv m_{60,i}=-5{\log }_{10}S{N}_{60,i},$$

$${\sigma }_i=\frac{2.5}{\ln 10}\frac{1}{S{N}_{60,i}}.$$

$$m_{60}(x_i)$$

$$\Rightarrow m_{60}(x_i)-m_{60}(\hat{x})$$

$$p(y_i|x_i,{\sigma }_i,\theta ,q_i=0)=\frac{1}{\sqrt{2\pi ({\sigma }_i^2+{\theta }_3^2)}}\exp \left(-\frac{[y_i-{\theta }_1-{\theta }_2(x_i-\hat{x}){]}^2}{2({\sigma }_i^2+{\theta }_3^2)}\right),$$

$$p(y_i|x_i,{\sigma }_i,\theta ,q_i=1)=\frac{1}{\sqrt{2\pi ({\sigma }_i^2+{\theta }_5^2)}}\exp \left(-\frac{[y_i-\hat{y}-{\theta }_4{]}^2}{2({\sigma }_i^2+{\theta }_5^2)}\right).$$

$$p(\mbox{\boldmath${y}$}|\mbox{\boldmath${x}$},\mbox{\boldmath${\sigma}$},\mbox{\boldmath${\theta}$})=\prod_{i=1}^{N}p(y_{i}|x_{i},\sigma_{i},\mbox{\boldmath${\theta}$})$$

$$p(y_{i}|x_{i},\sigma_{i},\mbox{\boldmath${\theta}$})=\sum_{q_{i}}p(q_{i})p(y_{i}|x_{i},\sigma_{i},\mbox{\boldmath${\theta}$},q_{i}).$$

$$\displaystyle p(\mbox{\boldmath${y}$}|\mbox{\boldmath${x}$},\mbox{\boldmath${\sigma}$},\mbox{\boldmath${\theta}$})=\prod_{i=1}^{N}$$

$$p(q_{i}|\mbox{\boldmath${y}$})\approx\frac{1}{M}\sum_{m=1}^{M}p(q_{i}|\mbox{\boldmath${y}$},\mbox{\boldmath${\theta}$}^{(m)}),$$

$$p(q_{i}|\mbox{\boldmath${y}$},\mbox{\boldmath${\theta}$})=\frac{p(q_{i})p(\mbox{\boldmath${y}$}|\mbox{\boldmath${x}$},\mbox{\boldmath${\sigma}$},\mbox{\boldmath${\theta}$},q_{i})}{Qp(\mbox{\boldmath${y}$}|\mbox{\boldmath${x}$},\mbox{\boldmath${\sigma}$},\mbox{\boldmath${\theta}$},q_{i}=0)+(1-Q)p(\mbox{\boldmath${y}$}|\mbox{\boldmath${x}$},\mbox{\boldmath${\sigma}$},\mbox{\boldmath${\theta}$},q_{i}=1)}.$$

$$p(\mbox{\boldmath${y}$}|\mbox{\boldmath${x}$},\mbox{\boldmath${\sigma}$},\mbox{\boldmath${\theta}$},q_{i}=0)=\prod_{i=1}^{N}p(y_{i}|x_{i},\sigma_{i},q_{i}=0)$$

$$\delta X=X_{\mathrm{upper}}-X_{\mathrm{lower}}$$

$$\sin q$$

$$\cos q$$

$$f(x)={\int }_{-\sqrt{1-x^2}}^{\sqrt{1-x^2}}(1-u+u\mu (x,y))dy,$$

$$\delta v_{\mathrm{r}}=v_{\mathrm{eq}}\frac{{\int }_{-1}^{1}xf(x)dx}{{\int }_{-1}^{1}f(x)dx}.$$

$$g(x,\varphi )={\int }_{-\sqrt{1-x^2}}^{\sqrt{1-x^2}}I(x,y)(1+ϵ(-x\sin \varphi +y\cos \varphi ))dy.$$

$$\delta v_{\mathrm{r}}$$

$${\omega }_{\mathrm{rot}}(\sin b)=A+B{\sin }^{2}b+C{\sin }^{4}b,$$

$$\sin b=y\sin i+\sqrt{1-x^2-y^2}\cos i.$$

$$\delta v_{\mathrm{rot}}=\frac{{\int }_{-1}^{1}xh(x,\varphi )dx}{{\int }_{-1}^{1}g(x,\varphi )dx}{R}_{\odot }\sin i$$

$$h(x,\varphi )={\int }_{-\sqrt{1-x^2}}^{\sqrt{1-x^2}}I(x,y){\omega }_{\mathrm{rot}}(x,y)(1+ϵ(-x\sin \varphi +y\cos \varphi ))dy.$$

$$\delta v_{\mathrm{r}}=\frac{{\int }_{-1}^1(xh(x,\varphi )R_{\odot }\sin i+c(x,\varphi ))dx}{{\int }_{-1}^{1}g(x,\varphi )dx}-\frac{{\int }_{-1}^{1}xg(x,\psi )dx}{{\int }_{-1}^{1}g(x,\psi )dx}{R}_{\odot }{\omega }_{\mathrm{orb}}.$$

$$\sin \mathrm{P}{\mathrm{A}}_{\mathrm{ecl}}=\sin \epsilon \cos \alpha$$

$$c(x,\varphi )={\int }_{-\sqrt{1-x^2}}^{\sqrt{1-x^2}}I(x,y)v_{\mathrm{conv}}(\mu (x,y))(1+ϵ(-x\sin \varphi +y\cos \varphi ))dy.$$

$$v_{\mathrm{conv}}(x,y)=0.4-1.2{\mu }^2+0.6{\mu }^4{\mathrm{km}\mathrm{s}}^{-1},$$

$$\delta v_{\mathrm{c}}={\int }_{-1}^{1}c(x,\varphi )dx$$

$$\delta v_{\mathrm{r}}$$

$${\omega }_{\oplus }(t)=\frac{2\pi }{P_{\oplus }}\frac{a_{\oplus }^2}{r^2(t)}\sqrt{1-e^2}$$

Full text

††pubyear: 2019††pagerange: Three years of Sun-as-a-star radial-velocity observations on the approach to solar minimum.–B

1 Introduction

The Sun is the only star that can be observed with spatial resolution fine enough to discern the finest convective and magnetic elements that decorate its surface. Tiny though they are, these small surface elements have a profound impact on the integrated solar spectrum. The contrast between the hot, upwelling cores of granules only a few hundred km across and the cooler surrounding downflow lanes gives rise to global spectral-line asymmetries (Dravins, Lindegren, & Nordlund, 1981). These asymmetries are strongly suppressed by small-scale magnetic fields in the faculae that surround large sunspot groups (Cegla et al., 2013). The finite number of granules, and their finite lifetimes, give rise to statistical fluctuations in the global solar radial velocity (Ludwig, 2006). The changing filling factor of sunspots and faculae as the Sun rotates gives rise to larger perturbations in global radial velocity and line asymmetries (Meunier, Desort, & Lagrange, 2010).

Early campaigns to monitor the solar radial velocity in integrated sunlight (e.g. Deming et al. 1987) were motivated by the need to understand the intrinsic variability of solar-type stars at a time when high-precision studies of stellar radial velocities (e.g. Campbell & Walker 1979) were in their infancy. At about the same time, the Birmingham Solar Oscillations Network (BiSON, Chaplin et al. 1996) started a 39-year campaign of radial-velocity monitoring of integrated sunlight, with the goal of using low-order $p$-modes in the solar oscillation spectrum to probe the Sun’s deep interior.

In the era of ultra-high precision radial velocity instruments such as ESO’s High-Accuracy Radial-velocity Planet Searcher (HARPS, Pepe et al. 2004) and its northern counterpart HARPS-N (Cosentino et al., 2012), Keck’s High Resolution Échelle Spectrometer (HIRES, Vogt et al. 1994), the Echelle SPectrograph for Rocky Exoplanets and Stable Spectroscopic Observations (ESPRESSO, Pepe et al. 2014) at the VLT, the EXtreme PREcision Spectrometer (EXPRES, Jurgenson et al. 2016) at Lowell Observatory and others following up transiting terrestrial-sized planet candidates from the CoRoT, Kepler/K2 and TESS space photometry missions, the need has become acute to understand the frequency spectrum and origins of stellar radial-velocity variability. Studies such as the recent “radial-velocity fitting challenge" of Dumusque et al. (2017) are designed to test the efficacy of novel analysis tools for modeling star-induced RV variability as part of the measurement process, and have highlighted the need to understand the underlying physics.

Using an approach developed by Fligge, Solanki, & Unruh (2000) for modelling solar irradiance variations, Meunier, Desort, & Lagrange (2010) pioneered the study of the effects of different types of solar activity on the global solar radial velocity. By partitioning solar images from the Michaelson Doppler Imager (MDI) aboard SoHo into quiet-sun, sunspot and facular regions according to their continuum brightness and magnetic flux density, and isolating the relative velocities of the different components in the Dopplergrams, Meunier et al. concluded that the dominant contributor to solar radial-velocity (RV) variability is convective suppression of the granular blueshift. This prediction was confirmed observationally by Haywood et al. (2016), who used the same technique on images from the Helioseismic and Magnetic Imager (HMI) on the Solar Dynamics Observatory (SDO) to model radial-velocity variations in integrated sunlight reflected from asteroid 4/Vesta.

The HARPS-N solar telescope (Phillips et al., 2016; Dumusque et al., 2015) was conceived with a longer-term goal in mind: to monitor precisely the solar radial velocity during the day using the same instrument as is being used at night for measuring the reflex orbital motions of exoplanet host stars. By using the stellar data-reduction pipeline for the solar spectra, we aim to characterise the impact of both solar variability and instrumental and data-reduction systematics on the radial velocities delivered by the instrument.

The analysis of the solar data is complicated by two considerations that do not apply to stellar targets. The observer and the target are both participants in solar-system gravitational dynamics, and the Sun is not a point source. The purpose of the present paper is to describe the methods used to correct for these non-stellar effects, transforming the data to the equivalent of the sidereal, heliocentric frame. The goal is to separate the effects of solar photospheric physics from those of solar-system dynamics and differential atmospheric effects.

In Section 2 we describe briefly the instrument, the observing strategy and the data-reduction pipeline. We develop a Gaussian mixture-model approach to determine daily extinction coefficients and to quantify the reliability of data points affected by short-term obscuration of parts of the solar disc. We correct the radial velocities for differential extinction across the solar rotation profile and use the good data within single days to assess the level of residual $p$-mode and granulation noise on minutes-to-hours timescales. In Section 3 we transform the radial velocities delivered by the pipeline from the barycentric frame (which is dominated by the synodic radial-velocity signal of Jupiter) to the heliocentric frame. We provide algorithms for correcting line-profile moments of the HARPS-N cross-correlation function (CCF) for the Earth’s changing orbital velocity and the obliquity of the solar rotation axis to the ecliptic. In Section 4 we analyse the behaviour of the radial velocity and CCF profile parameters on timescales from days to years. We identify fluctuations in the width and asymmetry of the CCF profile with the passages of sunspot groups and large facular regions across the solar disc, respectively. We identify correlations and temporal offsets between these parameters and the radial velocity, and discuss their viability as proxy indicators for disentangling exoplanetary orbital reflex motion from host-star activity.

2 HARPS-N solar telescope observations

2.1 Instrument and observing strategy

The HARPS-N solar telescope comprises a small guided telescope on an amateur mount, housed in a perspex dome on the exterior of the enclosure of the 3.58-m Telescopio Nazionale Galileo (TNG) at the Observatorio del Roque de los Muchachos, Spain. Its 7.6-cm achromatic lens of 200mm focal length feeds sunlight via an integrating sphere and an optical fibre into the calibration unit of the HARPS-N spectrograph. The telescope, fibre feed and control systems are described by Phillips et al. (2016). Early results from the instrument were published by Dumusque et al. (2015), confirming that it achieves uniformity of throughput better than 1 part in $10^{4}$over the solar disc. This performance is essential to achieve velocity precision of 10 cm s*-1* across the width of the solar rotation profile.

The instrument observes the Sun continuously each day, with 5-minute integration times designed to average out solar $p$-mode oscillations. The seasonal variations in hour angle are shown in Fig. 1. Observing begins each clear day when the Sun rises above an altitude of $\sim 20$ degrees, and ends either at the same altitude or earlier, if the instrument is required for afternoon setup by the night-time observer.

2.2 HARPS-N data reduction pipeline

The data are reduced using the same HARPS-N Data-Reduction System (DRS) employed for night-time stellar radial velocimetry. Calibration exposures taken after the end of solar observations and before the start of night-time observing each day provide order-by-order information on the locations of the echelle orders and the wavelength calibration scale. Fabry-Perot exposures are recorded simultaneously with the solar exposures to monitor instrumental drift, via a second optical fibre from the spectrograph calibration unit. Following optimal extraction (Horne, 1986; Marsh, 1989) to obtain one-dimensional background-subtracted spectra in each order, the data are calibrated in wavelength. They are then cross-correlated with a digital mask (Baranne et al., 1996; Pepe et al., 2002) derived from a typical G2 stellar spectrum, and corrected for instrumental drift derived from the Fabry-Perot spectrum. The radial velocity is computed as the mean velocity of a Gaussian fit to the CCF profile. The (dimensionless) contrast of the CCF is expressed as the maximum depth of the fitted Gaussian, expressed as a percentage of the surrounding pseudo-continuum. The area $W$ of the fitted Gaussian is then proportional to the product of the contrast $C$ and the full width at half maximum depth ($F\equiv{\rm FWHM}$):

[TABLE]

The units of $W$ and $F$ are the same, in this case velocities in km s*-1*, because the pseudo-continuum of the CCF is normalised to unity. Although $C$ is normally expressed as a percentage, the expression above treats it as a fraction $0<C<1$. The “area” is therefore defined in a manner similar to the equivalent width of a single spectral line.

The CCF itself is not perfectly symmetric, mainly because of the convective asymmetry in solar and stellar line profiles resulting from the brightness and velocity structure of the photospheric granulation pattern. This asymmetry is quantified as the difference in the line-bisector velocity in the upper and lower parts of the CCF, commonly referred to as the Bisector Inverse Slope (BIS; Queloz et al. 2001).

The DRS computes the formal uncertainty in each radial velocity by propagating the photon-noise error of the spectrum through all stages of extraction and cross-correlation, into the CCF. The error of the RV precision is then measured using the derivative of the CCF (Bouchy, Pepe, & Queloz, 2001). The median radial-velocity precision achieved with a standard 5-minute exposure is 0.43 ms -1 in observing conditions of high transparency; the upper 95th percentile is 0.77 m s*-1*.

2.3 Data quality assessment

Under normal circumstances the SNR of the radial velocity is determined primarily by photon shot noise and wavelength calibration error. Large systematic errors in the radial velocity can, however, arise when part of the solar disc is obscured, for example by cloud or by temporary tracking errors in the telescope mount. This happens because the solar disc has a finite angular diameter and a projected equatorial rotation speed of order 2 km s*-1*. Partial obscuration of the solar disc distorts the rotation profile and corrupts the measured radial velocity. It is therefore mandatory to identify and flag data points with anomalously low fluxes relative to neighbouring points.

The worst outliers are eliminated by rejecting the 5 percent of the data with the poorest estimated radial-velocity errors, followed by an iterative 6-sigma clip of the radial velocities in the heliocentric frame.

It is relatively straightforward to perform more rigorous automated data-quality assessment because the spectral fluxes observed in cloud-free conditions should follow an exponential extinction law as a function of airmass. To determine the apparent magnitude of the Sun, we use the SNR estimate for pipeline order 60 (echelle order 98, central wavelength 6245 Å) recorded in the HARPS-N data headers. The SNR is proportional to the square root of the recorded photon count. We therefore construct a sequence of instrumental magnitudes of which the $i$th measurement is

[TABLE]

with corresponding magnitude uncertainty

[TABLE]

On a day of near-optimal observing conditions, such as the clear spring day affected by occasional scattered clouds illustrated in Fig. 2, we expect the relation between the instrumental magnitude $m_{60}(x_{i})$ in order 60 and the airmass $x_{i}$ of the centre of the solar disc to be represented by a linear (Bouguer’s Law, Bouguer 1729) model

[TABLE]

The slope $k_{60}$ of the linear extinction law is the primary extinction coefficient, while $m_{60}(x=0)$ is the instrumental magnitude extrapolated to zero airmass using the extinction law. We use the inverse variance-weighted mean airmass $\hat{x}$ and the corresponding magnitude $m_{60}(x=\hat{x})$ as the origin for the linear regression, to eliminate correlation between the slope and the zero-point of the extinction law.

In addition to the fiducial magnitude $m_{60}(\hat{x})$ at the mean airmass and the extinction coefficient $k_{60}$ we include a white-noise variance parameter $\sigma^{2}_{\rm jit}$ representing low-level transparency fluctuations. This mitigates the risk of rejecting usable data in observing conditions where the transparency is imperfect but slowly-varying, as is the case shown in Fig. 3 for a day affected by the calima conditions often encountered in summer when fine dust is blown westward from the Sahara desert over the observatory. Even on a clear day, the extinction coefficient may vary subtly with time and/or hour angle (Poretti & Zerbi, 1993).

It is not clear that simple methods such as iterative sigma-clipping can make a reliable distinction between good and bad data under calima conditions in particular. Instead, we adopt a more rigorous and robust Bayesian mixture-model approach (e.g. Hogg, Bovy, & Lang, 2010) to identify those data points that lie close enough to the daily linear extinction law to be considered reliable, and to distinguish them from a separate background population of outliers that do not. For this purpose we introduce a discrete classifier $q_{i}$ which takes the value 0 for a good (“foreground”) observation, and 1 for a “background” outlier.

The extinction model yields the likelihood that a single measurement $y_{i}=m_{60}(x_{i})$ was obtained in good (i.e. $q_{i}=0$) conditions:

[TABLE]

where we have defined a vector of extinction-model parameters ${\theta}$ such that $\theta_{1}\equiv m_{60}(\hat{x})$, $\theta_{2}\equiv k_{60}$, and $\theta_{3}\equiv\sigma_{\rm jit}$.

Outliers produced by clouds or guiding errors generally lie below the linear relation, and are treated as having been drawn from a distinct “background" statistical population with mean $b_{\rm bg}$ (relative to the inverse variance-weighted mean $\hat{y}$) and variance $\sigma^{2}_{\rm bg}$.

For a point drawn from the background (i.e. $q_{i}=1$) outlier population, and an extended vector of model parameters \mbox{\boldmath{\theta}}\equiv\{m_{60}(\hat{x}),k_{60},\sigma_{\rm jit},b_{\rm bg},\sigma_{\rm bg}\} such that $\theta_{4}=b_{\rm bg}$ and $\theta_{5}=\sigma_{\rm bg}$, the likelihood is

[TABLE]

As Hogg, Bovy, & Lang (2010) and Foreman-Mackey (2014) point out, the marginal likelihood for the full dataset can be expressed as

[TABLE]

where $N$ is the number of observations on the day in question, and

[TABLE]

The simple prior $p(q_{i}=0)=Q$ and $p(q_{i}=1)=1-Q$, yields the likelihood function

[TABLE]

(Hogg, Bovy, & Lang, 2010).

In the present application, the prior $Q$ can be thought of as the fraction of all observations obtained on a given day that belong to the foreground (“good”) population. It is a quantity that must be derived from the data themselves, because weather conditions change from day to day.

We sample from the distribution described by Eq. 9 using a simple Markov-chain Monte Carlo scheme, saving the values of $p(y_{i}|x_{i},\sigma_{i},\theta,q_{i}=0)$ and $p(y_{i}|x_{i},\sigma_{i},\theta,q_{i}=1)$ for each sample for later use. This allows us to estimate the conditional foreground probabilities marginalised over the model parameters via the method described by Foreman-Mackey (2014):

[TABLE]

where $M$ is the number of samples in the Markov chain. In this expression, the conditional probability that a given data point $i$ is good ($q_{i}=0$) or bad ($q_{i}=1$) for a specified parameter set $\theta$ is given by

[TABLE]

As before, the priors $p(q_{i})$ take values $Q$ and $1-Q$ for $q_{i}=0$ and 1 respectively, while

[TABLE]

using Eq. 5 (and similarly for $q_{i}=1$ using Eq. 6).

The marginalised mixture membership probabilities described by Eq. 10 (referred to hereafter as “mixture probabilities”) provide a robust way of quantifying data quality. Fig. 4 illustrates the performance of the model on a day of generally poor and variable transparency, probably caused by cirrus cloud. Points lying along the upper envelope of the airmass plot in the upper-right panel follow a linear extinction law, indicating that they were obtained in clear intervals between cloud passages. As in Figs. 2 and 3, the radial velocities of the observations classified as having mixture probabilities greater than 0.9 of being “good” (blue) show a much lower point-to-point scatter than those with lower mixture probabilities (cyan, green, yellow, red). This confirms that the mixture-model approach is successful, and that a threshold probability of about 0.9 is appropriate for detecting and filtering out observations with anomalous extinction high enough to bias the radial-velocity measurements.

The prior probability distributions for the parameters are listed in Table 1. To ensure orthogonality between the fiducial magnitude $\theta_{1}=m_{60}(\hat{x})$ and the extinction coefficient $\theta_{2}=k_{60}$, we subtract the inverse variance-weighted mean airmass $\hat{x}$ from all airmass values, and the corresponding fiducial magnitude $m_{60}(\hat{x})$ from all magnitude values. The extinction coefficient was allowed to vary between 0.075 mag/airmass, which represents a transparency slightly better than is encountered in any of the 904 daily datasets examined, and 0.4 mag/airmass, representing the worst transparency encountered in calima conditions that can still yield a meaningful airmass plot. Short-term transparency fluctuations with a standard deviation as great as 0.08 mag/airmass are allowed in the foreground population; anything worse than this is likely to give differential extinction across the solar disc bad enough to distort the radial velocity measurements. The mean of the background outlier population was determined to be up to 0.8 magnitude fainter than the daily mean, with a standard deviation between 0.1 and 0.7 magnitudes. The quality of a given daily sequence is not known a priori, so the prior probability $Q$ that any given data point belongs to the foreground population is allowed to vary between 0.0 and 1.0.

Because the parameters are very close to being mutually orthogonal, we used an MCMC scheme employing a simple Metropolis-Hastings sampler. Two initial burn-in phases of $10^{4}$ samples were employed. At the end of each, the standard deviation of the parameter values was determined from the last $10^{3}$ steps to ensure better mixing in the next stage. A third production phase of $10^{4}$ steps yielded the posterior samples used to determine the daily extinction coefficients and to calculate the foreground membership probabilities for the individual data points. Median values of these samples are listed in machine-readable form in Appendix A, Table 2. The mixture probabilities from Eq. 10 provide clear thresholds against which the user can partition usable from unusable data, and are used as quality flags in Appendix B, Table 3.

2.4 Differential extinction corrections

The daily mean extinction coefficient defines the vertical extinction gradient across the diameter of the solar disk as a function of airmass. This gradient affects the measured solar radial velocity whenever the solar rotation axis is inclined relative to the vertical direction. This generally leads to a spurious redshift in the morning, when the approaching hemisphere of the Sun is lower in the sky than the receding hemisphere. The effect was first identified and corrected by Deming et al. (1987), and has subsequently been modelled by the BiSON team (Davies et al., 2014). The analytic approach we adopt here was developed independently of, and is slightly different from, these previous studies, but the results obtained are essentially identical.

The difference in airmass $X$ between the upper and lower limbs of the Sun is

[TABLE]

where $X=\sec{z}[1-0.0012(\sec^{2}{z}-1)]$ (Young & Irvine, 1967). The zenith angles $z_{\rm lower,upper}=z_{\odot}\pm\sin^{-1}(R_{\odot}/d_{\odot})$ depend on the zenith angle $z_{\odot}$ of the solar centre, the solar radius $R_{\odot}$ and the Earth-Sun distance $d_{\odot}$ at the time of observation.

The corresponding difference in magnitudes of atmospheric extinction is $\delta m_{60}=k_{60}\delta X/2$ across the solar radius. We define a fractional difference in brightness $\epsilon\simeq\delta m_{60}\ln{10}/2.5$, with $\epsilon=0$ at the centre of the solar disc.

The angle between the solar rotation pole and the local vertical is $\phi=PA-q$, where $PA$ is the position angle of the solar rotation axis relative to the north celestial pole at the time of observation. The parallactic angle $q$ is the angle subtended at the Sun by the celestial pole and the zenith, giving

[TABLE]

where $H$ is the solar hour angle and $\beta$ is the latitude of the TNG.

For a linear limb-darkening law with limb darkening coefficient $u$, the form of the rotation profile in the absence of differential extinction is (Gray, 1992):

[TABLE]

where the brightness is $I(x,y)=1-u+u\mu(x,y)$ and the foreshortening cosine $\mu(x,y)\equiv\sqrt{1-x^{2}-y^{2}}$.

The orthogonal Cartesian coordinates $x$ and $y$ are in the directions of increasing solar rotation velocity and the solar north pole respectively in the plane of the sky. For solid-body rotation with equatorial velocity $v_{\rm eq}$, the brightness-weighted correction to the solar radial velocity is

[TABLE]

In the absence of differential extinction across the rotation profile, the integrand in the numerator of Eq. 17 is an odd function, so the brightness-weighted velocity offset is zero.

We modify the calculation of the limb-darkened rotation profile to include such a density gradient, parameterised by the fractional extinction difference $\epsilon$ across the solar radius. A brightness gradient of uniform slope $\epsilon$ in $x$ and $y$ at an angle $\phi$ to the local vertical modifies the rotation profile in Eq. 17 to a new form:

[TABLE]

The sign of the $x\sin{\phi}$ term should be negative if $\phi$ is positive west of the meridian.

The velocity bias $\delta v_{\rm r}$ introduced by the extinction gradient thus depends on $v_{\rm eq}$, $u$, $\phi$, and $\epsilon$ according to

[TABLE]

So far we have assumed that the solar photosphere rotates as a solid body, ignoring the effects of differential rotation, axial inclination and centre-to-limb variations in convective blueshift. Reiners et al. (2016) modelled the effects of these modifications to the solar rotation profile on the solar radial velocity in the presence of an extinction gradient. Following their example, we can include differential rotation in our model by allowing the apparent rotational angular frequency $\omega_{\rm obs}$ to vary as a function of heliographic latitude $b$. This alters eq. 17 such that the intrinsic solar equatorial rotation speed $v_{\rm eq}$ is replaced by $R_{\odot}\omega_{\rm rot}(x,y)$. We adopt the sidereal differential rotation law of Snodgrass & Ulrich (1990) as a function of heliographic latitude $b$, defining

[TABLE]

with $A=2.972\times 10^{-6}\ {\rm rad\ s}^{-1}$, $B=-0.484\times 10^{-6}\ {\rm rad\ s}^{-1}$ and $C=-0.361\times 10^{-6}\ {\rm rad\ s}^{-1}$. When the solar rotation axis is inclined at angle $i$ to the line of sight, a point at location $(x,y)$ on the visible solar hemisphere has

[TABLE]

The integration for $v_{\rm rot}$ becomes

[TABLE]

where

[TABLE]

To take the Earth’s instantaneous orbital motion and the centre-to-limb variation in convective radial velocity into account we must include additional terms:

[TABLE]

Here $\psi$ is the position angle between ecliptic north and the local vertical at the Sun’s location. It always differs from $\phi$ by less than the 7-degree tilt of the solar rotation axis relative to the ecliptic. The position angle of ecliptic north relative to celestial north at the Sun is given by

[TABLE]

where $\varepsilon$ is the 23.4-degree obliquity of the ecliptic and $\alpha$ is the solar RA at the time of observation. By analogy with the position angle of the solar rotation pole, $\psi={\rm PA_{ecl}}-q$ where $q$ is the parallactic angle. The integration is specified separately in Eq. 24 in a Cartesian coordinate system in which $y$ is oriented toward the ecliptic north pole and $x$ is in the direction of Earth’s orbital motion. The solution has the same form as Eq. 19 with $v_{\rm eq}$ replaced by $-R_{\odot}\,\omega_{\rm orb}$.

The convective contribution is

[TABLE]

Both $I(x,y)$ and $v_{\rm conv}(\mu(x,y))$ are azimuthally-symmetric functions of radius only. The contribution of $c(x,\phi)$ to $\delta v_{\rm r}$ in Eq. 24 is therefore independent of $\epsilon$ and $\phi$, being a constant that depends on the radial variations in limb darkening and convective velocity only. To illustrate this point by way of example, we approximate the centre-to-limb dependence of the line-of-sight convective velocity as a sigmoidal function with the form

[TABLE]

the coefficients being chosen to approximate Fig. 18 of Beeck et al. (2013) and to give a derivative with respect to $\mu$ equal to zero at the centre and limb. The double integral

[TABLE]

has the analytic solution $\delta v_{\rm c}=0.376991-0.273631u$ for any linear extinction gradient across the solar disc. Since there is no dependence on either $\epsilon$ or $\phi$, the extinction gradient across the convective velocity pattern does not produce any change in radial velocity with either airmass or position angle. We therefore treat the effect of the centre-to-limb variation in convective velocity as a time-invariant offset to the velocity zero-point, and evaluate the integrals in Eqs. 22 and 24 to obtain

[TABLE]

The velocity correction for differential extinction varies almost linearly with airmass, and has an amplitude of 2 to 3 m s*-1* in an uninterrupted day’s observations. As $\cos{i}$ varies through the year, the values of $\delta v_{\rm r}$ depart by only 1 or 2 mm s*-1* from the solid-body rotation values obtained using a sidereal rotation rate $\omega_{\rm sid}=2.875\times 10^{-6}\ {\rm rad\ s}^{-1}$. The inclination dependence involves only even powers of $\cos i$, and does not change sign when the Earth crosses the solar equator.

The radial-velocity plots in the lower-left panels of Figs. 2, 3 and 4 have all been corrected for this effect, using Eq. 29 with $u=0.6$. The full set of extinction-corrected instrumental radial velocities are plotted in Fig. 7 for the duration of the observing campaign to date.

2.5 Intra-day solar radial-velocity variability

Phillips et al. (2016) have shown that the 5-minute solar $p$-mode oscillations are clearly detected in sequences of 20-s exposures taken as part of the science verification programme for the HARPS-N solar telescope. Although the 5-minute exposures used in the longer-term study reported here are designed to suppress these oscillations, the daily radial-velocity traces in the lower-left panels of Figs. 2 and 3 show that short-term correlated noise is still present. The daily RMS scatter shown in the upper panel of Fig. 5 is nonetheless only slightly greater than the 0.43 m s*-1* estimated uncertainty on a single radial-velocity measurement in conditions of good transparency. This agrees well with early results studied by Dumusque et al. (2015).

To ensure that the correlated noise is not caused by variable atmospheric transparency in calima conditions, in the lower panel of Fig. 5 we plot the same values against the daily extinction coefficient. Calima is characterised by values of $k_{60}$ in the approximate range from 0.15 to 0.40 magnitudes per airmass. No dependence of the daily RMS scatter on extinction coefficient is apparent.

To obtain the autocorrelation function of the Sun on minutes-to-hours timescales in the quietest possible conditions, we selected the period from 2018 July 21-October 12. Daily SDO magnetograms show that six weak or decaying bipolar magnetic regions crossed the centre of the solar disc during this period. These disc-centre crossings are marked by blue arrows in the upper panel of Fig. 6. We estimated the RMS scatter of the individual 5-minute radial-velocity measurements under these quiet-Sun conditions by taking all days on which we obtained more than 4 hours’ data (48 observations) with mixture probabilities (cf. Eq. 10) greater than 99 percent. Within each day we applied an iterative 4-sigma clip to remove very occasional outlying RV measurements with causes other than transparency fluctuations, which had not been rejected by the mixture modelling.

Even during periods when no discernible magnetic regions are present, the day-to-day scatter of the radial velocities greatly exceeds the intra-day scatter. The BIS and FWHM, however, display day-to-day scatter no greater than their daily averages. HARPS-N observations entail comprehensive daily calibration sequences. These establish the geometry for the extraction of the spectral orders from the CCD frames, and the wavelength dispersion relation for the extracted spectra. The accuracy of the zero point of the radial velocities therefore depends on the accuracy of the daily calibration. Our results indicate that the zero-point accuracy of the daily calibrations is very similar to the 0.4 m s*-1* RMS scatter measured by Dumusque (2018) in HARPS data reduced with the same version of the DRS used to extract and calibrate the solar spectra. We conclude that the daily jumps in the RV are likely to be a data-reduction artifact.

To eliminate the effects of day-to-day calibration zero-point errors and any residual rotationally-modulated stellar activity, we subtracted the daily mean from each day’s RV measurements. The resulting set of velocities is plotted in the middle panel of Fig. 6. The autocorrelation function (ACF) of these data was computed using the method of Edelson & Krolik (1988), and is plotted in the lower panel of the same figure. For comparison, the ACF of white noise with the same variance, sampled at the same times, is also shown.

The ACF confirms clearly that the radial velocities are strongly correlated, with a correlation half-life close to 15 minutes. This is somewhat longer than the average lifetime of a solar granule. The radial velocities become effectively uncorrelated for time lags longer than about 0.07 day (1.7 hours).

The 5-minute cadence of the observations is designed to suppress short-term variability caused by the solar 5-minute $p$-mode oscillations. For stellar observations, the HARPS and HARPS-N Rocky Planet-Search projects use 2 or 3 exposures of at least 15 minutes, spaced at least 2 or 3 hours apart. The 15-minute exposure mitigates the effects of $p$-mode oscillations, and the 2 to 3-hour sample spacing ensures that granulation noise is effectively uncorrelated (Dumusque et al., 2011; Meunier et al., 2015). The autocorrelation properties of the HARPS-N solar data illustrate the need for such a strategy.

3 Corrections for solar-system dynamics

3.1 Solar barycentric motion

The radial velocities derived from the solar spectra are processed using the same HARPS-N Data Reduction System that is employed for night-time observations. The exposure-weighted mid-time of each exposure is converted from UTC to Barycentric Julian Date in the Barycentric Dynamical Time standard, and the radial velocity is transformed to the reference frame of the solar-system barycentre. For both calculations, the apparent J2000 position of the centre of the Sun replaces conventional J2000 stellar coordinates.

When the target is the Sun, the use of the Barycentric Julian Date remains valid to within 2 or 3 seconds. The radial velocity, however, represents the component of the Sun’s barycentric motion in the direction of the observer, in this case HARPS-N at the TNG. The raw radial velocities thus comprise the radial component of the total solar reflex motion induced by the planets. The sum of these motions is clearly seen in Fig. 7 as a sinusoidal variation with an amplitude of 12 m s*-1* and a period close to 13 months. This is the synodic period of Jupiter observed from Earth.

To eliminate this reflex motion, we transform the radial velocities to the heliocentric system, using the JPL Horizons ephemeris (Giorgini et al., 1996) to obtain the component of the Sun’s barycentric motion in the direction of the TNG. This is subtracted from the raw velocities, reducing the observed velocities to a reference frame in which the rate of change of distance of the solar centre from the TNG is zero. The resulting radial velocities are shown in the lower panel of Fig. 7. The average observed velocity of 102 m s*-1* includes the instrumental zero point, the solar gravitational redshift and the granular convective blueshift, and agrees with the mean values determined by Molaro et al. (2013) and Haywood et al. (2016) using HARPS on the ESO 3.6-m telescope within 2 m s*-1*.

3.2 Earth orbital eccentricity and obliquity

The FWHM of the Gaussian fit to the CCF is known to be rotationally modulated in exoplanet host stars with moderate to high levels of stellar activity (Queloz et al., 2009). The periodogram of the observed FWHM of the solar CCF, constructed with the Bayesian Generalised Lomb-Scargle (BGLS) formalism of Mortier et al. (2015), however, displays a clear additional modulation on timescales of 6 months and 1 year (Fig. 8).

These additional modulations are of dynamical origin. The solar apparent equatorial velocity is $v\sin{i}=R_{\odot}\omega_{\rm obs}\sin{i}$, where $i$ is the inclination of the solar rotation axis to the line of sight. The solar apparent angular velocity $\omega_{\rm obs}$ is the difference between the solar sidereal angular velocity $\omega_{\odot}$ and the instantaneous orbital angular velocity $\omega_{\oplus}(t)$ of the observer. The Earth’s orbital eccentricity imposes an annual variation on the latter,

[TABLE]

where $e$ is the Earth’s orbital eccentricity, $P_{\oplus}$ is the orbital period and $r(t)/a_{\oplus}$ is the Earth-Sun separation at time $t$ expressed in au.

The variation in $v\sin i$ with time is straightforward to calculate given the Cartesian ecliptic direction vector $\hat{\boldmath{\omega}}_{\odot}$ of the solar rotation pole and the instantaneous Sun-Earth direction vector $\hat{\boldmath{r}}_{\oplus}$ in ecliptic coordinates, since $\sin i=\sqrt{1-(\hat{\boldmath{\omega}}_{\odot}.\hat{\boldmath{r}}_{\oplus}})^{2}$. The correction needed to recover the value of the FWHM that would be observed in the sidereal reference frame is, however, less straightforward. The observed CCF can be considered as the convolution of a limb-darkened rotation profile with the CCF that would be obtained for a hypothetical non-rotating star. For the purposes of this study we treat the observed FWHM $F_{\rm obs}$ as the quadrature sum of the intrinsic linewidth $F_{0}$ and an unknown fraction $\gamma$ of the observed $v\sin i$, such that

[TABLE]

An observer at interstellar distance in the Sun’s equatorial plane would observe the sidereal FWHM

[TABLE]

We find that the RMS scatter in the resulting values of $F_{\rm sid}$ is minimised for $\gamma=1.04$. Using this value we plot the corrected FWHM in the sidereal frame in Fig. 9. To first order, the CCF profile is the convolution of a local line profile with a stellar rotation profile of unit area, re-normalised to the resulting continuum level. So long as the area of the local line profile does not change strongly between facular and quiet-Sun areas, and the centre-to-limb variation in the area of the local line profile is small, the area $W$ of the CCF will remain approximately constant even when parts of the stellar disc are occupied by dark spots or bright faculae. The area of the fitted profile is thus proportional to the product of the FWHM and the contrast, so $C_{\rm sid}F_{\rm sid}=C_{\rm obs}F_{\rm obs}$ gives the corresponding correction to the contrast.

4 Line-profile behaviour on the approach to solar minimum

When corrected to the sidereal reference frame, the FWHM of the solar CCF shows a barely significant increase over the duration of the study. Dumusque (2018) found a similar slow change in the FWHM for HARPS, which they attributed to long-term changes in the focus of the instrument. The FWHM shows strong short-term variability, however, and peaks strongly at the times when one or more major spot groups are near the centre of the disk (Fig. 9, top two panels; SILSO World Data Centre 2015–2018; Clette et al. 2016).

The Bisector Inverse Slope, plotted in the third panel of Fig. 9, exhibits a secular decline over the three years of observation amassed to date, with a shorter modulation clearly seen at the solar rotation period and its first harmonic. Although it is likely that the BIS will be affected to some extent by the modulation of the solar $v\sin i$, Fig. 9 shows that it is not the dominant signal.

This long-term trend is mirrored in the area of the Gaussian fit to the CCF, plotted in the bottom panel of Fig. 9.

4.1 Long-term trends

Early in the sequence the BIS shows rapid oscillations, as might be expected if multiple active regions cross the disc during every solar rotation. Underlying the peaks associated with disc-centre passages of large magnetic regions, the BIS exhibits a long-term downward trend. This trend appears to mirror the general decline in sunspot number over the same period.

The CCF area shows a corresponding secular increase over the three years of observation amassed to date. This increase is not seen in the FWHM, so cannot be attributed to, for example, a long-term change in spectrograph focus. On shorter timescales of order weeks to months, the FWHM is anticorrelated with the CCF contrast. The CCF area, being the product of the contrast and the FWHM, does not therefore contain a strong rotational signal. The secular increase in the CCF area thus appears to be driven by a slow increase in the CCF contrast. This might happen if, for instance, the relative contributions of the blue and red echelle orders were changing slowly with time, as might happen through uncorrected ageing of a calibration lamp’s colour temperature. We can rule out this possibility because the HARPS-N data-reduction system corrects the blaze function and normalises the continuum level in all orders before computing the CCF, with the explicit purpose of preventing such chromatic imbalances.

The trends in both the area and the BIS are interrupted temporarily at approximately BJD 2457850.0, and resume from their new levels from about BJD 2457900.0 (Fig. 9). This interruption coincides with the appearance of the large, persistent bipolar active regions mentioned above (AR12653 and AR12654 in late 2017 March, then AR12674 and AR12675 early in 2017 September), following a lengthy period of low activity. From late 2017 onward, the decline toward solar minimum resumed, with few spot groups visible.

At face value, a secular increase in the equivalent widths of the thousands of spectral lines that contribute to the CCF area gives the appearance of a global decrease in the stellar effective temperature. The area of the CCF represents a reliable proxy of the temperature and metallicity of a star, as demonstrated by Malavolta et al. (2017). Under the reasonable assumption that the Sun did not change its metallicity during the three years spanning our observations, we can use the CCF area to track global temperature changes in the photosphere of the Sun. The correlation with the BIS suggests that the trends are also associated with a change in line asymmetry, and hence with changes in either the magnetic field (Brandt & Solanki, 1990; Meunier et al., 2017) or the flux effect due to dark spots (Meunier, Desort, & Lagrange, 2010; Haywood et al., 2016).

Unlike the BIS, however, the CCF area shows little evidence of rotational modulation. This lack of rotational modulation suggests that the structures responsible for the change in CCF area are axisymmetrically distributed over the solar surface. We speculate that the CCF area may be tracing the evolution of the solar magnetic network, which largely comprises dispersed magnetic flux elements from decayed active regions being swept up towards the poles by meridional circulation. As the activity cycle declines, the supply of flux elements to the network decreases, reducing localised facular brightening and giving the appearance of a decreased effective temperature.

A decline in network flux also increases the dominance of the asymmetric quiet-Sun pattern of granular line asymmetry (e.g. Cegla et al. 2013). The changes in line asymmetry measured by the BIS show both the long-term trend seen in the CCF area and a rotationally-modulated component that responds to the visibility of the facular magnetic fields of large bipolar regions.

The connection between the effective temperature and the coverage of activity regions and network is not yet fully explored, but the solar data already offer some useful insights. The visibility of the more strongly RV-suppressive facular magnetic-field concentrations in SDO images has been found to be a stronger predictor of the solar radial velocity than the total visible-hemisphere coverage of active-region and network field elements (Milbourne et al., 2019). Milbourne et al. (2019) find the filling factor of the small, weakly RV-suppressive magnetic areas that they label “network” to be greater (and around solar minimum, much greater) than the filling factor of the more strongly RV-suppressive faculae. If the CCF area is sensitive to the temperature contrast and filling factor of network/faculae in the optical, changes in the CCF area will mostly be driven by the slow, small, long-term variation of network rather than the much less dominant (by area) faculae. The BIS, on the other hand is sensitive to velocity suppression, so it would see both the long-term changes (due to changes in network and facular area) and the short term effects of rotationally modulated faculae.

4.2 Short-term line-profile variability

The FWHM and the BIS both peak strongly at times when major sunspot groups and their accompanying bipolar magnetic regions cross the solar disc for the first time. Unlike the FWHM, however, the BIS shows a persistent periodic signal at the solar rotation period for about 3 rotations following the first passage of each major new spot group.

This behaviour is particularly apparent following the appearance of the pair of large bipolar solar active regions AR12654 and AR12655, which crossed the solar disc centre near JD 2457845.5 (2017 March 30). One rotation later, these regions had decayed to the extent that only remnants of the leading spots are visible in continuum images from SDO/HMI 111https://www.solarmonitor.org/. The SDO/HMI magnetograms, however, show that the bipolar fields from these active regions persisted for a further 3 or 4 rotations.

The same behaviour is seen following disc centre passage of the large bipolar regions AR12674 and AR12675 near JD 2458002.5 (2017 September 5). These same two spot groups survived to cross disc centre one Carrington rotation later on 2017 October 1, when they produced secondary peaks in both the FWHM and the BIS. A third “echo" of this event appears strongly in the BIS plot (Fig. 9) after a further Carrington rotation, even though the spots themselves had almost disappeared. SDO magnetograms taken on 2017 October 28 show that the bipolar magnetic regions associated with these groups are still very much present.

The coincidence of strong, isolated peaks in the FWHM time series with similarly isolated peaks in the sunspot number, and the greater observed persistence of the BIS signal for several rotations after the emergence of each new active region suggests strongly that the FWHM responds to dark spots while the line asymmetry measured by the BIS traces inhibition of convection in active-region faculae. To place this observation on a more quantitative footing we computed the autocorrelation functions of the timeseries data in RV, FWHM, BIS and CCF area. For each time series, we took the daily median value and time of observation, and subtracted a linear trend in time. We used the method of Edelson & Krolik (1988) to obtain the ACF of the resulting irregularly-sampled timeseries, with a bin width of 2 days.

The results are presented in Fig. 10, for time lags from 0 to 150 days. The BIS shows the most persistent periodicity, with peaks occurring at multiples of 27 days, out to 110 days (4 solar rotations). The RV and FWHM show clear recurrence peaks out to 55 days (2 solar rotations). The autocorrelation function of the CCF area shows a steady decline with lag. There is little evidence of periodicity.

4.3 Active-region BIS signal and radial velocity

As discussed above, the BIS shows a secular trend mirroring the smooth long-term change in the CCF area. To isolate the rotational signal in the BIS we decorrelated the BIS against the CCF area. The residual BIS signal exhibits annual and six-month periodicities similar to those seen in the FWHM (Fig. 11). We fitted out these periodicities using sinusoids with periods of 6 months and 1 year, whose phases and amplitudes are to be treated as free parameters. The corrected BIS signal is shown in the lower panel of Fig. 11.

Periodograms of the daily median RV and de-trended BIS time series bear a close resemblance to each other (Fig. 12). Both periodograms are dominated by the solar rotation period and its first harmonic.

4.4 RV proxies and time lags

Given the close resemblance of the RV and BIS periodograms, the similarity of the patterns of peaks in the RV and BIS in fig. 9 and the evidence described above that the oscillatory part of the BIS tracks the active-region facular filling factor, it is tempting to infer that the BIS will be a good proxy indicator of the radial velocity. Fig. 13, however, shows that the two are almost entirely uncorrelated.

This lack of correlation arises from the temporal behaviour of the RV and the line-profile parameters. Fig. 14 shows the sunspot number, corrected FWHM, radial velocity and uncorrected BIS during the period between 16 June and 3 November 2017, when a number of large, relatively isolated active regions crossed an otherwise quiet Sun. Close inspection of this figure, and of the discrete cross-correlation functions (Edelson & Krolik, 1988) in Fig. 15, reveals that the maxima in radial velocity occur 1 to 3 days before the maxima in FWHM and BIS respectively. Given that each peak has a width of about 7 days, this temporal offset between the RV and BIS weakens the correlation.

Similar temporal shifts between maxima in the RV, BIS, FWHM and chromospheric Ca II H&K emission flux have been reported for other intensively-observed stars such as GJ674 (Bonfils et al., 2007), GJ176 (Forveille et al., 2009), CoRoT-7 (Queloz et al., 2009) and HD41428 (Santos et al., 2014). While the authors of these papers generally attribute the phase lags in these cases to the RV signal being dominated by dark spots, the HARPS-N solar data indicate that the same phenomenon is present even when the RV signal is dominated by facular suppression of the convective blueshift (Meunier, Desort, & Lagrange, 2010; Haywood et al., 2016). It is also worth noting that HD41248 is an inactive G2V star, so the observation by Santos et al. (2014) that increases in RV precede increases in FWHM by 2-3 days is consistent with the solar behaviour, i.e., facular, not spot, dominated RV variations.

The reasons for this temporal offset between the RV and the line-profile parameters are subtle but important. The FWHM is expected to increase when a dark sunspot group crosses the centre of the solar disc. This phenomenon is well known from Doppler imaging of starspots (Vogt & Penrod, 1983): when the missing light from a dark spot crosses the stellar rotation profile, it reduces the depth of the centre of the normalised line profile while approximately preserving the area (subject to any centre-to-limb variation in the area of the local CCF). The width of the line must therefore increase to compensate. The lack of persistence seen in the FWHM in Figs. 9 and 10 supports an association of the FWHM with dark spots rather than the longer-lived surrounding magnetic regions.

In the absence of granular convection, we might expect the BIS to track the first derivative of the FWHM as the spectral signatures of dark spots cross the line profile from blue to red. Such an effect may be present at a low level, but in reality the BIS peaks only slightly later than the FWHM when an active region crosses disc centre. As noted above, the BIS peaks repeat for 3 or 4 solar rotations after the initial appearance of a large, isolated active region, even after the spots themselves have dispersed. This supports an association of the line asymmetry (quantified by the BIS) with the suppression of granular convection in the strongly-magnetic, bright facular regions that are the longer-lived counterparts of bipolar sunspot groups.

The close connection proposed by Meunier, Desort, & Lagrange (2010) between the solar radial velocity and magnetic suppression of granular convection has been confirmed observationally by Haywood et al. (2016) and Milbourne et al. (2019). The temporal offset between the RV and BIS signatures thus appears all the more surprising, until we consider how the measurements are actually made.

The cause of these temporal offsets was discussed by Dumusque, Boisse, & Santos (2014) and is illustrated in their Figure 6. The HARPS-N DRS determines radial velocities by fitting a Gaussian to the CCF and recording the location of the peak. The CCF profile represents the convolution of many different local line profiles with the limb-darkened stellar rotation profile. As discussed above, a dark spot crossing the disc creates a disturbance to the CCF that is symmetric in time about disc-centre crossing if the CCF profile itself is reasonably symmetric. If we now replace the dark spot with a highly-simplified active-region CCF which contributes similar total flux but whose line profile is displaced redward by suppression of convective blueshift, the net effect is to add a travelling perturbation whose form resembles the first derivative of the local line profile. Being antisymmetric, this disturbance influences the Gaussian fit, and hence the measured radial velocity, in a way that is not symmetric about disc-centre passage.

For a slowly-rotating star like the Sun, the rotation profile is narrower than the intrinsic linewidth. When the active region is on the approaching hemisphere, the disturbance weakens the blue wing of the CCF and moves the core redward. When the active region is at disc centre, the core of the line moves redward. When on the receding hemisphere, the disturbance boosts the red wing of the profile and weakens the blue side of the core. The displacement of the centre of the fitted Gaussian is therefore always redward, but is not necessarily symmetric about the time of disc-centre passage. The time lags between the maxima of the RV, the BIS and the FWHM thus depend crucially on the details of changes to the local line asymmetry that occur as the active region crosses the disc.

There is another effect that could introduce time delays between different activity indicators. Strong, mostly vertical (buoyant) magnetic fields inhibit, in particular, flows perpendicular to them. These however, are not visible at disc centre. Their area and limb darkened visiblity peaks $\sim 40$ degrees off disc centre. This leads to a natural phase offset between BIS (which is more affected by magnetic inhibition of velocities, peaking away from from disk center, and suppressed relative to local quiet Sun) and FWHM (which is more affected by thermal variations, peaking at disc centre). One implication is that the lag would then be, at least in part, a function of the rotation period.

Figure  14 also shows that at times when the Sun is quiet, the day-to-day scatter in the FWHM and BIS is comparable to the intra-day scatter. The same is not true of the radial velocity, which shows an additional day-to-day scatter of order 0.4 m s*-1*. The most likely source of this additional scatter is lack of repeatability in the daily calibrations used to establish the order locations and wavelength dispersion relations for the 2-D spectra.

5 Discussion and Conclusions

Three years of radial-velocity monitoring the Sun as a star with the HARPS-N instrument reveal variability on timescales ranging from minutes to days to years.

Within a single day, the radial-velocity signal is strongly correlated on timescales of minutes to hours. We find that 15-minute integrations attenuate the autocorrelation amplitude of the data by a factor 2 or so better than 5-minute sampling. The same autocorrelation analysis shows that samples taken 2 or more hours apart are effectively statistically independent. This supports the recommendations of the earlier study of the impact of stellar granulation by Dumusque et al. (2011), whose conclusions have determined the nightly observing strategy of surveys for terrestrial-mass planets around bright stars with HARPS and HARPS-N.

Finally, the solar data reveal a need for more sophisticated data-reduction strategies. High-precision instruments such as HARPS-N require comprehensive daily calibration sequences to monitor the long-term stability of the instrument. No calibration is perfect, however, and the solar data reveal radial-velocity zero-point errors with an RMS scatter of order 0.4 m s*-1*. These arise mainly from day-to-day zero-point uncertainty in the daily wavelength calibrations carried out by the version of the pipeline with which these data were reduced. A recent major upgrade to the HARPS/HARPS-N data-reduction system has introduced a better-constrained 2-dimensional fit across all échelle orders into these daily calibrations for stellar data. Improvements can therefore be expected in future data releases from this project, once the 58000 existing solar spectra and future observations have been reprocessed.

Long-term changes in the area of the CCF appear to arise from cycle-related changes in the magnetic network. There is no rotational modulation component to this signal, and the impact on the radial velocity is not yet fully clear. There appears to be no strong trend in RV over the three years of this study, but we caution that we have not yet had the opportunity to observe the full range of activity in a solar cycle.

The long-term evolution of the BIS, which measures the core-wing asymmetry of the CCF, has a baseline which mirrors the long-term evolution of the CCF area. Bipolar magnetic active regions passing across the solar disk cause perturbations lasting several days in both the radial velocity and the CCF profile parameters. The line-shape parameters of the CCF appear to respond to different components of the active regions. The FWHM shows modulation behaviour that tracks the sunspot number. The BIS shows more persistently repetitive modulation, suggesting that the line asymmetry responds primarily to disruption of the photospheric granulation pattern in the faculae, which are the longest-lived components of active regions.

There is no single proxy indicator that can recover the RV from an isolated measurement of the line-profile shape or other activity indicators without reference to time-domain information. Given that active regions in solar-type stars evolve on timescales comparable to the rotation period, it is therefore mandatory to sample stellar radial-velocity signals often enough that active-region passages can be fully characterised in all proxy indicators. The radial velocity tracks the BIS closely, but the BIS lags the radial velocity by about 3 days. As Figs. 14 and 15 show, these time lags are large enough that a simple linear decorrelation against BIS or FWHM cannot be used to correct the radial velocity for the effects of magnetic activity. A more sophisticated approach which uses time-domain information as well as line-shape information is needed. Although such a treatment is beyond the scope of this paper, we refer the reader to recent papers by Aigrain, Pont, & Zucker (2012), Rajpaul et al. (2015) and Jones et al. (2017). These multivariate Gaussian-process regression approaches use temporal derivative models of RV proxy indicators in a general way that has strong potential capability to model such time lags.

Acknowledgements

The HARPS-N project has been funded by the Prodex Program of the Swiss Space Office (SSO), the Harvard University Origins of Life Initiative (HUOLI), the Scottish Universities Physics Alliance (SUPA), the University of Geneva, the Smithsonian Astrophysical Observatory (SAO), and the Italian National Astrophysical Institute (INAF), the University of St Andrews, Queen’s University Belfast, and the University of Edinburgh. ACC acknowledges support from the Science & Technology Facilities Council (STFC) consolidated grant number ST/R000824/1 and UKSA grant ST/R003203/1, and the support of the visiting scientist program at Lowell Observatory, where part of this work was carried out. DP and RW acknowledge partial support by NASA award number NNX16AD42G. XD is grateful to the Branco Weiss Fellowship Society in Science for its financial support. SHS is grateful for support from NASA Heliophysics LWS grant NNX16AB79G. CAW acknowledges support by STFC grant ST/P000312/1. The research leading to these results received funding from the European Union Seventh Framework Programme (FP7/2007-2013) under grant agreement number 313014 (ETAEARTH). This work was performed in part under contract with the California Institute of Technology (Caltech)/Jet Propulsion Laboratory (JPL) funded by NASA through the Sagan Fellowship Program executed by the NASA Exoplanet Science Institute (RDH). DWL acknowledges partial support from the Kepler mission under NASA Cooperative Agreement NNX13AB58A with the Smithsonian Astrophysical Observatory. LM acknowledges the support by INAF/Frontiera through the "Progetti Premiali" funding scheme of the Italian Ministry of Education, University, and Research. Some of this work has been carried out within the framework of the NCCR PlanetS, supported by the Swiss National Science Foundation. FP further acknowledges support by the Swiss National Science Foundation SNSF through grant number 166227.

Appendix A Daily extinction and mixture model parameters

Table 2 lists Julian dates and daily extinction and mixture-model parameters for each day of observation. The Julian date is the median of the observation dates on the day to which each record applies. The values of the extinction coefficient $k_{60}$, the magnitude zero point $m_{60}(x=0)$ above the atmosphere, the background population mean magnitude offset $b_{\rm bg}$ and standard deviation $\sigma_{\rm bg}$, the transparency fluctuation amplitude $\sigma_{\rm jit}$ and the prior mixture fraction $Q$ are the medians of the final MCMC samples computed for each day. Only the first sixteen days of the campaign are listed here. The full table is available online-only, in machine-readable form as a comma-separated variable (CSV) file.

Appendix B Radial velocities and CCF profile parameters

Table 3 provides the full set of radial velocities and CCF profile parameters. It is provided online-only in machine-readable form as a CSV file. The columns are:

(1) Barycentric JD;

(2) JD;

(3) Barycentric RV in km s*-1*;

(4) Heliocentric RV in km s*-1*;

(5) Heliocentric RV in km s*-1*, corrected for differential extinction;

(6) RV error in km s*-1*;

(7) Observed CCF FWHM in km s*-1*;

(8) Sidereal CCF FWHM in km s*-1*, corrected for Earth orbital motion and solar obliquity;

(9) CCF contrast (percent);

(10) Observed CCF BIS in km s*-1*;

(11) Apparent hour angle in radians;

(12) Apparent Airmass;

(13) Pseudo-magnitude derived from SNR in order 60;

(14) Magnitude error derived from SNR in order 60;

(15) Solar declination $\delta$ (radian);

(16) Solar angular diameter (radian);

(17) Position angle of north solar rotation pole (radian);

(18) Mixture-model foreground membership probability $p(q_{i}=0)$.