2 Statistical methodology

  Denote the radial velocity (RV) at time t, measured at wavelength $\lambda$by $v_{t\lambda}$ and the true underlying RV by v_t. Further, let $\beta_0$ be the average RV.  A two-stage statistical model is then

$$\begin{eqnarray} v_{t\lambda_j}&=&v_t+\varepsilon_{t\lambda_j},\\ v_t&=&\beta_0+f(t\vert\gamma)+a(t),\end{eqnarray}$$ (1) (2)

($j=1,~\dots,n$), where $\varepsilon_{t\lambda_j}$ is a normally distributed error term with variance $\sigma^2_{v_t}$, due to measurement error and wavelength calibration, $f(t\vert\gamma)$describes the orbital motion where $\gamma$ groups the orbital parameters, and a(t) is a periodic fluctuation of v_t due to the pulsation of the star. A general expression for a(t) is  

$$a(t)=\sum_{k=1}^K\alpha_k\sin(\omega_kt+\phi_k).$$ (3)

Here, $\alpha_k$, $\omega_k$, and $\phi_k$ are the amplitude, the frequency, and the phase of the radial velocity due to the kth pulsation mode. It is assumed that the frequency of $f(t\vert\gamma)$ is much smaller than the pulsation frequencies $\omega_k$. This assumption ensures that it is reasonable to consider $f(t\vert\gamma)$ constant during a time period within which the sinusoidal terms complete their cycles. Clearly, Eqs. (1) and (2) can be combined into a single one:  

$$v_{t\lambda_j}=\beta_0+f(t\vert\gamma)+a(t)+\varepsilon_{t\lambda_j}.$$ (4)

When solely the estimation of v_t for fixed t is of interest, it is sufficient to obtain a number of replications at different wavelengths $\lambda_1,~\dots,\lambda_n$, whereafter they are averaged to yield $\overline{v}_t$. The variance is equal to  

$$\sigma^2_{v_t}=\frac{1}{n-1}\sum_{j=1}^n(v_{t\lambda_j}-v_t)^2.$$ (5)

Now, $\overline{v}_t$ is an unbiased estimator for v_t since the expectation of $\varepsilon_{t\lambda}$ is assumed to be zero. Its standard deviation is estimated by replacing v_t in (5) with $\overline{v}_t$, and taking the square root. In order to obtain the standard error we divide this expression further by n.

Moreover, $\overline{v}_t$ is unbiased for the orbital motion $\beta(t)=\beta_0+f(t\vert\gamma)$ since in addition a(t) is zero on average. In contrast, the variance estimators are different in both situations since $\sigma^2_{v_t}$ fails to acknowledge pulsational variability. Indeed, $\sigma^2_{v_t}$ is the variance of $v_{t\lambda}$ around v_t, while we are now interested in the spread of $v_{t\lambda}$ around $\beta(t)$.

Since the pulsational variability and the error term $\varepsilon_{t\lambda}$ are statistically independent, Eq. (4) yields

$$\begin{eqnarray} \sigma^2_{\beta(t)}&=&E[(v_{t\lambda}-\beta(t))^2]\\ &=&\mbox{v... ...&\mbox{var}[a(t)]+\frac{1}{n-1}\sum_{j=1}^n(v_{t\lambda_j}-v_t)^2.\end{eqnarray}$$ (6) (7) (8)

When interest lies in the accuracy of the estimator $\widehat{\beta(t)}= \overline{v}_t$ one should compute

$$\begin{eqnarray} \mbox{var} ({\widehat{\beta(t)}})&=&E[(\overline{v}_{t}-\beta(t... ... \right)\nonumber\\ &=&\mbox{var}[a(t)]+\frac{1}{n}\sigma^2_{v_t}.\end{eqnarray}$$ (9)

These variances will be used as weights for the determination of the unknown orbital parameters.

Available data can consist of both full pulsation cycles as well as measurements at a single time t (but at n different wavelengths $\lambda_j$). In the second case, the second term in (9) is straightforward to evaluate but the first one is not. Indeed, even when we can assume that $\beta(t)$ can be considered constant during a pulsation cycle, a(t) is generally non-zero. Unfortunately, measurements at a single t provide no information about the discrepancy between the two. In order to overcome this problem, we will propose a method to determine the first variance component from external information.

In the case of a monoperiodic pulsation, K=1 in (3). When more than one sine term is present, it is often reasonable to assume that the function a(t) can be approximated by a single sine function during one night, but that different approximations would be necessary for different nights. In other words, we will assume that the sinusoidal amplitude $\alpha$ is constant during one night, but fluctuates in a complicated fashion over longer periods of time. To model this, we assume that the amplitudes are drawn randomly from a population of amplitudes with mean $\mu$ and variance $\tau^2$.This implies that, if several values for $\alpha$ have been obtained, $\alpha_1,~\dots,\alpha_m$ say, one can estimate the average $\mu$,$\overline{\mu}$ say, as well as the variance

$${\tau}^2=\frac{1}{m-1}\sum_{\ell=1}^m(\alpha_\ell-\mu)^2.$$

These results can be used to obtain the variance of a(t). Now, a(t) is a sinusoidal function of t, however with variable amplitude $\alpha$. To account for this extra source of variation, we use a fundamental result in mathematical statistics (Bickel & Doksum 1977, p. 36):  

$$\mbox{var}(a)=\mbox{var}[E(a\vert\alpha)]+E[\mbox{var}(a\vert\alpha)],$$ (10)

where $E(a\vert\alpha)$ is the expectation of a(t), given a value of $\alpha$ and similarly $\mbox{var}(a\vert\alpha)$ is the conditional variance of a random variable. The unconditional counterparts are denoted by E(.) and $\mbox{var}(.)$. Now, assuming time t is uniformly distributed within a pulsation cycle $[0,2\pi/\omega]$, the conditional mean of a(t), given an amplitude $\alpha$is

$$E(a\vert\alpha)= \frac{\omega}{2\pi}\int_0^{2\pi/\omega}\alpha\sin[\omega(t+\phi)]{\rm d}t=0$$

whence the first term in (10) cancels. Secondly, the conditional variance of a(t), given an amplitude $\alpha$ is

$$\mbox{Var}(a\vert\alpha)=\frac{\omega}{2\pi}\int_0^{2\pi/\omega} \alpha^2\sin^2[\omega(t+\phi)]{\rm d}t= \frac{\alpha^2}{2}$$

of which the mean is

$$\begin{eqnarray} E[\mbox{var}(a\vert\alpha)]&=&E\left(\frac{\alpha^2}{2}\right)\... ...box{Var}(\alpha)\nonumber\\ &=&\frac{1}{2}\mu^2+\frac{1}{2}\tau^2.\end{eqnarray}$$ (11)

Substituting (11) into (10) leads to the following variance formula for the estimated $\beta$: 

$$\mbox{var}({\widehat{\beta(t)}})= \frac{1}{2}\mu^2+\frac{1}{2}\tau^2 +\frac{1}{n}\sigma^2_{v_t},$$ (12)

which is estimated by plugging in estimated values for $\mu$, $\tau^2$, and $\sigma^2_{v_t}$. It is instructive to contrast this quantity with the variance of a single measurement $v_{t\lambda}$ about the average RV $\beta(t)$: 

$$\mbox{var}(v_{t\lambda})= \frac{1}{2}\mu^2+\frac{1}{2}\tau^2 +\sigma^2_{v_t}.$$ (13)

When the number n of measurements $v_{t\lambda_k}$ ($k=1,~\dots,n$) increases, the determination of v_t is done with increasing precision and the third term in (12) approaches zero. In contrast, the intrinsic source of uncertainty, due to the sinusoidal fluctuation, cannot be circumvented.

Note that, for a star which is monoperiodic, K=1 and hence the same amplitude will be found during each observational period. Hence, $\tau^2\equiv 0$. In practice, one might still find a small but non-zero value for $\tau^2$ since the amplitudes will be determined with measurement error.