An Adaptive Scheduling Tool to Optimize Measurements to Reach a Scientific Objective: Methodology and Application to Measurements of Stellar Orbits in the Galactic Center

The goal of most measurements performed is to extract some scientific information by fitting a model to the observations. We will assume the measurements errors to be normally distributed and we will denote one measurement by (m_i, t_i, σ_i) where m_i is the type of measurement, t_i is the time of the measurement, and σ_i is the measurement uncertainty. In addition, we will denote the model used to analyze the measurements as $M(\{{m}_{i},{t}_{i}\};{\boldsymbol{p}})$ where the vector ${\boldsymbol{p}}$ denotes all the N_param parameters that are included in the fit. These include the parameters that are scientifically relevant but can also include other parameters such as nuisance parameters or noise parameters. In some cases, one is interested in comparing two models ${M}_{1}(\{{m}_{i},{t}_{i}\};{{\boldsymbol{p}}}_{1})$ and ${M}_{2}(\{{m}_{i},{t}_{i}\};{{\boldsymbol{p}}}_{2})$ (which do not necessarily depend on the same number of parameters) in order to identify the one favored by the data.

In addition, we will denote the set of existing measurements by ${{ \mathcal O }}_{\mathrm{existing}}=\{({m}_{i},{t}_{i},{\sigma }_{i})\}$ and the set of possible (future) measurements by ${{ \mathcal O }}_{\mathrm{possible}}=\{({m}_{j},{t}_{j},{\sigma }_{j})\}$. Note that the set of possible future measurements is discrete and takes into account experimental/observational constraints (observational windows, maintenance time for experimental facilities, etc.).

The question we are addressing in this communication can be formulated as follows: given a set of existing measurements, how can the measurements from a set of possible future measurements be prioritized in order to increase the sensitivity with respect to a given scientific objective?

Three types of scientific objective will be considered: (i) a case where one is interested by maximizing the global quantity of information (the information entropy) related to the full posterior distribution, similarly to what has been developed by Ford (2008) and Loredo et al. (2012), (ii) a case where one focuses on one scientific parameter p_k and where one wants to minimize the uncertainty ${\delta }_{{p}_{k}}$ related to the estimation of that specific parameter, and (iii) a case where one wants to maximize the ratio of the Bayesian evidence between two models in order to maximize the chances of discriminating one model from the other.

The implementation of the adaptive scheduling tool presented in this paper is based on the Fisher matrix. This implicitly assumes that the posterior probability distribution function can be approximated by a multivariate normal distribution. In particular, this is the case when the errors follow a Gaussian distribution and the model can be approximated by its linearization around the maximum likelihood. Under this assumption, the Fisher matrix is a good estimate of the inverse of the covariance matrix for unbiased estimators. The Fisher matrix, a square matrix of dimension N_param, can be computed by (see, e.g., Vallisneri 2008; Heavens 2009 or the appendix from Albrecht et al. 2006)

$$\begin{eqnarray}&&{\boldsymbol{F}}={{\boldsymbol{P}}}^{T}\cdot {\boldsymbol{P}},\end{eqnarray} \tag{ 1 }$$

where ${\boldsymbol{P}}$ is the matrix of the partial derivatives of the model with respect to the parameters

$$\begin{eqnarray}&&{\left[{\boldsymbol{P}}\right]}_{{ik}}=\displaystyle \frac{1}{{\sigma }_{i}}\displaystyle \frac{\partial M(i;{\boldsymbol{p}})}{\partial {p}_{k}},\end{eqnarray} \tag{ 2 }$$

where σ_i corresponds to the uncertainty of the ith measurement. The matrix of the partial derivatives ${\boldsymbol{P}}$ has dimension ${N}_{\mathrm{measurements}}\times {N}_{\mathrm{param}}$.

An estimation of the uncertainties on the estimated parameters resulting from a fit of the model M using a set of measurements $\{\left({m}_{i},{t}_{i},{\sigma }_{i}\right)\}$ can be determined from the covariance matrix ${\boldsymbol{\Sigma }}$, which is the inverse of the Fisher matrix (see, e.g., Vallisneri 2008; Heavens 2009 or the appendix from Albrecht et al. 2006):

$$\begin{eqnarray}&&{\boldsymbol{\Sigma }}={{\boldsymbol{F}}}^{-1}.\end{eqnarray} \tag{ 3 }$$

In particular, the marginalized uncertainties of the estimated parameters are given by the diagonal of this matrix (Heavens 2009):

$$\begin{eqnarray}&&{\delta }_{{p}_{l}}=\sqrt{{{\rm{\Sigma }}}_{{ll}}}.\end{eqnarray} \tag{ 4 }$$

It is important to note that this covariance matrix depends on (i) the measurement data set (the timing and the type of measurements and their related uncertainties {(m_i, t_i, σ_i)}) and (ii) the set of fitted parameters (i.e., ${\boldsymbol{p}}$).

In addition, Ford (2008) has shown that the expected difference in information entropy ${ \mathcal I }$ that is induced by an additional measurement k is given by (see Equation (25) from Ford 2008)

$$\begin{eqnarray}&&{{\rm{\Delta }}}_{k}{ \mathcal I }={\rm{ln}}\left(\displaystyle \frac{{\sigma }_{{\rm{pred}};k}}{{\sigma }_{k}}\right),\end{eqnarray} \tag{ 5 }$$

where σ_pred;k is the model uncertainty for the measurement k obtained while not including this measurement in the fit and σ_k is the expected measurement uncertainty. The notation Δ_kx denotes the difference in the quantity x induced by adding the measurement k to the existing data set. This relationship means that the optimal measurements in order to increase the global information entropy are those with large model uncertainties (or, quoting Loredo et al. 2012, "we learn the most by sampling where we know the least"), which is called maximum entropy sampling (Sebastiani & Wynn 2000).

Finally, under the assumption that the posterior is a multivariate normal distribution, the model uncertainty can directly be related to the inverse of the Fisher matrix through

$$\begin{eqnarray}&&\displaystyle \frac{{\sigma }_{\mathrm{pred};{k}}^{2}}{{\sigma }_{{k}}^{2}}={\tilde{{\boldsymbol{P}}}}_{{k}}^{T}\cdot {\boldsymbol{\Sigma }}\cdot {\tilde{{\boldsymbol{P}}}}_{{k}},\end{eqnarray} \tag{ 6 }$$

where ${\boldsymbol{\Sigma }}$ is the covariance matrix obtained when the measurement (m_k, t_k, σ_k) is not included in the fit and ${\tilde{{\boldsymbol{P}}}}_{k}$ is a column vector containing the partial derivatives of the measurement k with respect to the fitted parameters ${\boldsymbol{p}}$:

$$\begin{eqnarray}&&{\left[{\tilde{{\boldsymbol{P}}}}_{k}\right]}_{i}={\left[{\boldsymbol{P}}\right]}_{{ki}}=\displaystyle \frac{1}{{\sigma }_{k}}\displaystyle \frac{\partial M({m}_{k},{t}_{k};{\boldsymbol{p}})}{\partial {p}_{i}}.\end{eqnarray} \tag{ 7 }$$

One statistical tool to compare two models using a set of measurements is to use the Bayesian evidence, defined as the probability to obtain the data given a certain model (or in other words, the evidence is the normalization of the posterior) (see, e.g., Kass & Raftery 1995; Jaynes 2003; Liddle 2007; Gregory 2010; Gelman et al. 2013; Knuth et al. 2014). Under the assumption that the prior knowledge about the two models is uniform, the Bayesian evidence is also directly proportional to the probability of a certain model given the data set. Comparing the Bayesian evidences from fits using two different models is one way to compare these models and to discriminate them. Using the Laplace quadratic approximation, the Bayesian evidence for a model 1 is given by (Liddle 2007)

$$\begin{eqnarray}&&{{ \mathcal E }}_{1}={{ \mathcal P }}_{1,\max }\sqrt{\displaystyle \frac{{\left(2\pi \right)}^{{N}_{1}}}{\left|{{\boldsymbol{\Sigma }}}_{1}\right|}},\end{eqnarray} \tag{ 8 }$$

where ${{ \mathcal P }}_{1,\max }$ is the maximum of the posterior, N₁ is the number of parameters fitted in the model 1, and $\left|{{\boldsymbol{\Sigma }}}_{1}\right|$ is the determinant of the matrix ${{\boldsymbol{\Sigma }}}_{1}$. The difference of the log-evidence between two models is therefore given by

$$\begin{eqnarray}\begin{array}{rcl}\mathrm{ln}\displaystyle \frac{{{ \mathcal E }}_{1}}{{{ \mathcal E }}_{2}} & = & \mathrm{ln}\displaystyle \frac{{{\rm{\Pi }}}_{1,\max }}{{{\rm{\Pi }}}_{2,\max }}+\displaystyle \frac{{N}_{1}-{N}_{2}}{2}\mathrm{ln}2\pi +\mathrm{ln}\displaystyle \frac{{{ \mathcal L }}_{1,\max }}{{{ \mathcal L }}_{2,\max }}\\ & & +\ \displaystyle \frac{1}{2}\mathrm{ln}\left|{{\boldsymbol{F}}}_{1}\right|-\displaystyle \frac{1}{2}\mathrm{ln}\left|{{\boldsymbol{F}}}_{2}\right|,\end{array}\end{eqnarray} \tag{ 9 }$$

where ${{ \mathcal L }}_{{i},\max }$ and π_i,max are the likelihood and prior of the model i for the parameters that maximize the posterior. The first two terms (the ratio of the priors and the number of fitted parameters) do not depend directly on the scheduling of the measurements. On the other hand, the third part of this equation depends on the actual measurements themselves. This dependence makes it hard to find the sequence of measurements optimal to discriminate between two models (since the actual measurements are not known). One solution is to consider a Gaussian likelihood and to simulate data using one of the two models. Let us assume that model 1 is the model we are trying to confirm from the data and let us use this model to simulate the measurements so that the last equation becomes

$$\begin{eqnarray}\begin{array}{rcl}\mathrm{ln}\displaystyle \frac{{{ \mathcal E }}_{1}}{{{ \mathcal E }}_{2}} & = & \mathrm{ln}\displaystyle \frac{{{\rm{\Pi }}}_{1,\max }}{{{\rm{\Pi }}}_{2,\max }}+\displaystyle \frac{{N}_{1}-{N}_{2}}{2}\mathrm{ln}2\pi +\displaystyle \frac{1}{2}\mathrm{ln}\left|{{\boldsymbol{F}}}_{1}\right|-\displaystyle \frac{1}{2}\mathrm{ln}\left|{{\boldsymbol{F}}}_{2}\right|\ \\ & & +\ \displaystyle \sum _{i}\displaystyle \frac{{\left({M}_{1}({m}_{i},{t}_{i},{{\boldsymbol{p}}}_{1})-{M}_{2}({m}_{i},{t}_{i},{{\boldsymbol{p}}}_{2})\right)}^{2}}{2{\sigma }_{i}^{2}}.\end{array}\end{eqnarray} \tag{ 10 }$$

In conclusion, in this section we have shown how the Fisher matrix can be used to characterize the different objectives that one may want to optimize on:

• 1.  
  The global information entropy: it requires one to optimize over Equation (6), which depends on the inverse of the Fisher matrix.
• 2.  
  The uncertainty of one specific parameter of interest (${\delta }_{{p}_{l}}$): it requires one to optimize over one component of the diagonal of the inverse of the Fisher matrix.
• 3.  
  A criterion to optimize the chances of discriminating a model M₁ from a model M₂: one suggestion is to optimize over an estimate of the difference in the log of the Bayesian evidence, given by Equation (10). This quantity depends on the evaluation of the models themselves and on the determinant of their Fisher matrices.

In the next two subsections, we will discuss how to estimate efficiently the Fisher matrix, its determinant, and its inverse.

As can be noticed from Equation (2), one needs to compute the partial derivatives of the model with respect to the parameters in order to compute the Fisher matrix. Three methods can be used to compute these derivatives:

• 1.  
  Analytically compute the derivatives of the model. This method has the advantage of being quick to evaluate, but it requires long calculations and code. When the model relies on the integration of a set of ordinary differential equations of motion, this method also requires integration of the equations of variations (as an example, see Section 4 from Lainey et al. 2004).
• 2.  
  Use automatic differentiation. The main advantage of this method is that one does not need to derive explicitly the expression for the partial derivatives.
• 3.  
  Use finite differences. This method has the advantage of easy implementation (especially when the model requires the integration of differential equations) but has the disadvantage that it can be numerically unstable.

Each of these methods has pros and cons and they need to be chosen on a case-by-case basis considering each specific situation.

One of the main steps in the algorithm developed in this paper is to compute the Fisher matrix, its determinant, and its inverse. Indeed, the goal of the adaptive scheduling tool is to optimize either on one diagonal component of ${\boldsymbol{\Sigma }}$ or on its determinant $\left|{\boldsymbol{\Sigma }}\right|$. Therefore, we will need to evaluate these quantities for a huge number of different combinations of measurements, and it is important to compute them using a method that is computationally efficient. Recomputing the full Fisher matrix and inverting it for all the sets of measurements is computationally intensive. An alternative approach is to update the Fisher matrix as well as its determinant and its inverse.

In other words, if we consider the Fisher matrix ${{\boldsymbol{F}}}^{(n)}$ obtained for a set of measurements ${{ \mathcal O }}_{n}$, we can efficiently compute the Fisher matrix ${{\boldsymbol{F}}}^{(n+1)}$ which is obtained from the same set of measurements with an additional measurement (m_k, t_k, σ_k). The procedure is straightforward and is given by

$$\begin{eqnarray}&&{{\boldsymbol{F}}}^{(n+1)}={{\boldsymbol{F}}}^{(n)}+{\tilde{{\boldsymbol{P}}}}_{k}^{T}\cdot {\tilde{{\boldsymbol{P}}}}_{k}.\end{eqnarray} \tag{ 11 }$$

Following the same principle, the determinant of the Fisher matrix can easily be updated:

$$\begin{eqnarray}&&\left|{{\boldsymbol{F}}}^{(n+1)}\right|=\left|{{\boldsymbol{F}}}^{(n)}\right|\left(1+{\tilde{{\boldsymbol{P}}}}_{k}^{T}\cdot {{\boldsymbol{\Sigma }}}^{(n)}\cdot {\tilde{{\boldsymbol{P}}}}_{k}\right),\end{eqnarray} \tag{ 12 }$$

where ${{\boldsymbol{\Sigma }}}^{(n)}={\left[{{\boldsymbol{F}}}^{(n)}\right]}^{-1}$.

Finally, the procedure to update the covariance matrix (the inverse of the Fisher matrix) is given by

$$\begin{eqnarray}&&{{\boldsymbol{\Sigma }}}^{(n+1)}={{\boldsymbol{\Sigma }}}^{(n)}-{{\boldsymbol{U}}}_{k}^{(n)},\end{eqnarray} \tag{ 13 }$$

where ${{\boldsymbol{U}}}_{k}^{(n)}$ is the update matrix, which depends on the measurement (m_k, t_k, σ_k) and is given by

$$\begin{eqnarray}&&{{\boldsymbol{U}}}_{k}^{(n)}=\displaystyle \frac{\left({{\boldsymbol{\Sigma }}}^{(n)}\cdot {\tilde{{\boldsymbol{P}}}}_{k}\right)\cdot {\left({{\boldsymbol{\Sigma }}}^{(n)}\cdot {\tilde{{\boldsymbol{P}}}}_{k}\right)}^{T}}{1+{\tilde{{\boldsymbol{P}}}}_{k}^{T}\cdot {{\boldsymbol{\Sigma }}}^{(n)}\cdot {\tilde{{\boldsymbol{P}}}}_{k}},\end{eqnarray} \tag{ 14 }$$

where the vector ${\tilde{{\boldsymbol{P}}}}_{k}$ is given by Equation (7).

To summarize, adding one measurement (m_k, t_k, σ_k) to a given data set ${{ \mathcal O }}_{n}$ will lead to a change in the criterion that we are optimizing on, given by the following expressions.

• 1.  
  Optimization over the global information entropy (maximization):

  $$\begin{eqnarray}&&{\varepsilon }_{k}={{\rm{\Delta }}}_{k}{ \mathcal I }=\mathrm{ln}\left({\tilde{{\boldsymbol{P}}}}_{k}^{T}\cdot {{\boldsymbol{\Sigma }}}^{(n)}\cdot {\tilde{{\boldsymbol{P}}}}_{k}\right),\end{eqnarray} \tag{ 15a }$$

  where the column vector ${\tilde{{\boldsymbol{P}}}}_{k}$ is given by Equation (7).
• 2.  
  Optimization over one specific parameter of interest (minimization of the related uncertainty):

  $$\begin{eqnarray}&&{\varepsilon }_{k}=-{{\rm{\Delta }}}_{k}{\delta }_{{p}_{l}}^{2}={\left[{{\boldsymbol{U}}}_{k}^{(n)}\right]}_{{ll}},\end{eqnarray} \tag{ 15b }$$

  using Equation (14).
• 3.  
  Optimization for model comparison (maximization of the ratio of the Bayesian evidence):

  $$\begin{eqnarray}\begin{array}{rcl}{\varepsilon }_{k}={{\rm{\Delta }}}_{k}\mathrm{ln}\displaystyle \frac{{{ \mathcal E }}_{1}}{{{ \mathcal E }}_{2}} & = & \displaystyle \frac{{\left({M}_{1}({m}_{k},{t}_{k},{{\boldsymbol{p}}}_{1})-{M}_{2}({m}_{k},{t}_{k},{{\boldsymbol{p}}}_{2})\right)}^{2}}{2{\sigma }_{k}^{2}}\\ & & +\ \displaystyle \frac{1}{2}\mathrm{ln}\displaystyle \frac{1+{\tilde{{\boldsymbol{P}}}}_{1,k}^{T}\cdot {{\boldsymbol{\Sigma }}}_{1}^{(n)}\cdot {\tilde{{\boldsymbol{P}}}}_{1,k}}{1+{\tilde{{\boldsymbol{P}}}}_{2,k}^{T}\cdot {{\boldsymbol{\Sigma }}}_{2}^{(n)}\cdot {\tilde{{\boldsymbol{P}}}}_{2,k}},\end{array}\end{eqnarray} \tag{ 15c }$$

  where subscripts 1 and 2 refer to the two models to be compared.

To start the algorithm, we need a set of existing measurements ${{ \mathcal O }}_{\mathrm{existing}}=\left\{({m}_{i},{t}_{i},{\sigma }_{i})\right\}$ and a set of possible future measurements ${{ \mathcal O }}_{\mathrm{possible}}=\left\{({m}_{j},{t}_{j},{\sigma }_{j})\right\}$. We also need to know the scientific objective we are optimizing for: (i) maximize the global information entropy or (ii) minimize the uncertainty related to a given parameter p_j or (iii) maximize the odd ratio between two models.

Our implementation of the adaptive scheduling tool is iterative. There are two different ways to control the number of iterations. The first way is to chose arbitrarily the desired number of measurements (N) and to find a scheduling for these N measurements that is optimized with respect to a given criterion. The second way is to stop the iterations when we reach a given objective (for example when the uncertainty on an estimated parameter is below a given threshold). In this second case, the number of measurements N is free as well.

The algorithm consists in the following steps.

• 1.  
  Fit the model $M({m}_{i},{t}_{i};{\boldsymbol{p}})$ to the existing measurements ${{ \mathcal O }}_{\mathrm{existing}}$ in order to determine the optimal parameters ${{\boldsymbol{p}}}_{\mathrm{opt}}$ (in the case of a comparison of models, both models need to be fitted). For this step, the actual data need to be known, not just their uncertainties.
• 2.  
  Compute the partial derivatives of the model (of both models in the case of a comparison):

  $$\begin{eqnarray*}&&\displaystyle \frac{1}{{\sigma }_{i}}{\left.\displaystyle \frac{\partial M({m}_{i},{t}_{i};{\boldsymbol{p}})}{\partial {p}_{k}}\right|}_{{{\boldsymbol{p}}}_{\mathrm{opt}}}\end{eqnarray*}$$

  for all the existing and possible future measurements. This evaluation can be computationally intensive (depending on the model) but is required only once.
• 3.  
  Compute the starting covariance matrix ${{\boldsymbol{\Sigma }}}^{(0)}$ from the set of existing measurements ${{ \mathcal O }}_{\mathrm{existing}}$ by inverting the Fisher matrix, see Equation (3).
• 4.  
  Iteratively determine the set of N optimal measurements by using the following procedure (starting with n = 0 and with an empty list of optimized measurements ${{ \mathcal O }}_{\mathrm{opt}}=\{\}$):

  • (a)  
    loop on all the possible future measurements $({m}_{k},{t}_{k},{\sigma }_{k})\in {{ \mathcal O }}_{\mathrm{possible}}$ and for each of them compute the quantity ε_k from Equations (15) depending on the scientific objective considered.
  • (b)  
    find the measurement (m_k, t_k, σ_k) that maximizes ε_k, add it to the list of optimal measurements ${{ \mathcal O }}_{\mathrm{opt}}$, and remove it from the list of future possible measurements ${{ \mathcal O }}_{\mathrm{possible}}$.
  • (c)  
    update the covariance matrix by adding the measurement (m_k, t_k, σ_k) using Equation (13), increase n and iterate by going back to step (a) as long as the number of measurements in the optimal set of future observations is smaller than N, or as long as we have not reached the given threshold.

This algorithm is very efficient and allows one to prioritize quickly measurements based on a criterion related to a scientific objective.

The set of existing measurements used in this analysis consists in astrometric measurements performed at the Keck observatory between 1995 and 2017 and spectroscopic measurements performed at the Keck observatory and at the VLT between 2000 and 2017. These measurements are reported in the literature (Gillessen et al. 2017; Do et al. 2019), and here we consider them to be independent. The model $M(\{{m}_{k},{t}_{k}\};{\boldsymbol{p}})$ used for this analysis relies on the integration of the first post-Newtonian equations of motion and includes the Römer time delay and the first-order relativistic redshift. This model is fully described in the appendix of Do et al. (2019). There are 15 parameters that characterize this model: the six orbital parameters for the star S0-2, the SMBH mass parameter GM, the distance to the GC R₀, the 2D position and 3D velocity of the SMBH, a parameter that encodes a deviation from the relativistic redshift, and a parameter that encodes a deviation from the first post-Newtonian order in the equations of motion (or in other words that encodes a deviation from the relativistic advance of the periastron).

In this section, we give a full example of telescope time scheduling using the tool presented in Section 2, including considerations related to systematics and to weather and instrumental losses. We focus on the measurement of the relativistic redshift of S0-2 and therefore fix the parameter that encodes a deviation from the first post-Newtonian equations of motion to its relativistic value (i.e., to unity). Therefore, the model contains 14 parameters. The objective is to find a measurement sequence in 2018 that would optimize the chances of a successful detection of the relativistic redshift given the set of existing measurements (up to 2017). The set of future possible measurements considered comprises two distinct sets: the first one consists of astrometric measurements for every day between 2018 March 1 and September 15 and the second one consists of radial velocity measurements for every day in the same time period. The uncertainty for the future measurements has been chosen as 20 km s⁻¹ for the radial velocity and 0.9 mas for the astrometry.

The sequence of 10 measurements optimized in order to detect the relativistic redshift (we are optimizing on the redshift parameter) using the algorithm presented in Section 2 is presented in Table 1 and Figure 1. It is important to keep in mind that this result depends on the existing measurements: different existing measurements could have led to another sequence. This sequence also depends on the expected accuracy of future possible measurements and on the parameters that are included in the fit.

Zoom In Zoom Out Reset image size

Two interesting facts need to be noticed about the optimized sequence of measurements obtained in this analysis. First, 9 out of the 10 measurements are spectroscopic observations. This makes sense considering that the redshift impacts directly on the radial velocity but not the astrometry. A more surprising fact is that the adaptive scheduling tool does not favor measurements at the maximum of the redshift signal. Rather, it favors measurements that are at the radial velocity turning points (i.e., the maximum and minimum of the radial velocity curve, see Figure 1). This is due to correlations between the redshift parameter and the other model parameters, in particular with T₀ (the time of closest approach for S0-2) and with the SMBH mass parameter GM. These correlations are maximal exactly when the redshift signal is maximal, making that epoch not optimal in order to measure the redshift signal. In other words, the redshift signal can easily be absorbed by a small change in T₀ or GM. To demonstrate the effect of correlations, we run a test case and search for the optimized sequence of measurements if we fit for all the model parameters except for T₀ and GM. The orange measurements from Figure 1 are the radial velocity measurements optimized for that case using the algorithm from Section 2. When T₀ and GM are not included in the fit, these measurements are all located close to the maximum of the redshift signal, as expected.

It is interesting to quantify the impact of an optimized scheduling compared to a "naive" one. For this, let us consider the gain in the signal-to-noise ratio (S/N) of the relativistic redshift due to one RV measurement in 2018. The redshift S/N using measurements existing up to 2017 is 1.16. A naive scheduling consisting in performing one additional measurement at the time when the redshift signal is maximal induces an increase in S/N of 4.5% (the S/N being increased up to 1.205), while an optimized scheduling consisting in performing the measurement at the time when the RV is maximal increases the S/N by 335% (the S/N is now 4.5). This illustrates clearly the benefit of an optimized scheduling and how such a scheduling tool can help to increase the outcome of a project.

At this stage, the important point learned from our adaptive scheduling tool is that, from a statistical point of view, nearly all the power to measure the relativistic redshift comes from spectroscopic measurements around the two RV turning points. This result relies purely on statistical considerations, and practicalities and systematics need to be considered as well. First of all, other longer-term scientific objectives will be pursued with S0-2's measurements. It is important to realize that measurements around closest approach, while not of prime importance for the redshift, are important for measuring other scientific parameters such as, e.g., R₀ (see Section 4.4) and the relativistic advance of the periastron (see Section 4.3). Therefore, we included measurements at closest approach, an important event as well. In addition, two strategies have been used to control systematics for the relativistic redshift measurement: (i) we decided to use different telescopes (Keck and Gemini) to measure the radial velocities in order to assess hypothetical biases related to one instrument and (ii) we decided to also take astrometric measurements at the same epochs as radial velocities and vice versa (see the discussion in Section 3). Finally, the results from Table 1 have been obtained assuming a perfect knowledge of S0-2's orbital parameters. In practice, there is an uncertainty associated with these parameters. In particular, at the end of 2017, the uncertainty related to the 2018 time of closest approach T₀ was of the order of five days. In order to control for a possible bias in the estimate of T₀, we included measurements very early in 2018.

In addition, it is important to consider instrumental losses and weather losses. Using different telescopes (Keck and Gemini) for the capture of the important measurements is one way to mitigate the risk of missing important measurements because of instrumental losses (in addition to allowing us to check for systematics). In addition, we increased the number of requested measurements in 2018 in order to take into account possible measurements losses. A careful analysis of the Keck Observatory statistics shows that roughly 50% of GC measurement nights are lost because of weather or instrumental problems. This has been included in our scheduling by making a large number (10⁵) of simulations where each measurement has a probability of 50% of being lost and by assessing the distribution of the relativistic redshift S/N at the end of 2018.

To summarize, the scheduling of the 2018 GC measurement campaign for the GCOI team was based on: (i) measurements of radial velocities very early in 2018 to control our estimate of T₀; (ii) radial velocities at the two radial velocity turning points, the most important measurements to detect the redshift, using at least two different instruments to control for systematics; (iii) measurements at closest approach important for other scientific objectives; (iv) dual spectroscopic–astrometric measurements at the important epochs to control for systematics; (v) the total number of measurements being determined by requiring that the relativistic redshift S/N at the end of 2018 has a probability of at least 99.9% of being above 5, considering that each night has a 50% chance of being lost, and (vi) the exact scheduling being determined by considering a weather pattern lasting three days. The scheduling resulting from this analysis is presented at the top of Figure 2 and consists of 36 spectroscopic measurements (for Keck and Gemini) and 10 astrometric measurements. For the sake of comparison, the actual measurements performed in 2018 are presented in the lower panel of Figure 2. They consist of 19 radial velocity measurements and five astrometric measurements. The main sources of losses for 2018 are related to bad weather and to eruptions from the volcano in Hawaii, one of them being accompanied by an earthquake that damaged the instruments for a couple of days. The distribution of the relativistic redshift S/N at the end of 2018 using the scheduling described in this section and assuming that each measurement has a 50% chance of being lost is presented in the right panel of Figure 3. The probability of having an S/N below 5 is 0.03%. The S/N presented in this figure is a statistical S/N. No systematics are considered in this analysis and systematics will diminish the S/N calculated here. The analysis of the real data has shown that including systematics such as offsets between instruments, correlations in the astrometric positional uncertainty, or additional systematics uncertainties will produce a systematic uncertainty of the same order of magnitude as the statistical uncertainty (see Do et al. 2019 for the analysis of the measurements and a detailed study of systematics). The right panel of Figure 3 shows the evolution of the S/N with time in 2018. The blue curve corresponds to the evolution expected from the original scheduling while the orange curve corresponds to the evolution of the real measurements. Although nearly half of the measurements in 2018 have been lost, the final S/N obtained in the analysis remains far above 5. This is due to the careful planning and to the number of measurements scheduled but also to the fact that some of the measurements were actually better than anticipated.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

This first example shows the utility of the adaptive scheduling tool presented in Section 2 and how it can lead to surprising results. As a result of this analysis, the GCOI team decided to put more effort into measuring the radial velocity of S0-2 around the two turning points in order to maximize the chances of a successful detection of the relativistic redshift (Do et al. 2019). We also developed the details behind the GCOI's decision to include systematics and possible weather/instrumental losses in the scheduling of the 2018 measurement campaign.

The relativistic redshift is only the first relativistic effect measurable using stellar orbits in the GC. The next effect expected to be measured is the well-known advance of periastron, which is due to the first post-Newtonian relativistic correction to the equations of motion. In this section, we are interested in prioritizing measurements in order to measure this effect. For this reason, we fix the redshift parameter to its relativistic value (i.e., to unity) but we let the parameter that encodes a deviation from the first post-Newtonian equation of motion be free. As a result, the model again contains 14 parameters. The set of "existing measurements" consists of astrometric and radial velocity measurements of S0-2 up to 2018. The set of future possible measurements consists of astrometric and spectroscopic measurements for each year between March and September (the yearly window for which GC measurements are feasible from Earth) on a daily basis. The uncertainty considered for the future measurements is 0.2 mas for the astrometry and 20 km s⁻¹ for the radial velocity.

As a first step, we run the adaptive scheduling tool considering possible measurements between 2018 and 2058. The distribution of the measurements from an optimized sequence of 30 observations is presented in Figure 4. While of limited interest, this figure shows three interesting features. First, it shows that, as anticipated, the detection of the advance of the periastron relies mainly on the astrometry. Second, it shows that the optimal strategy to measure the advance of the periastron given a certain measurement window is to take measurements as late as possible. This is easily explained considering that this effect is a secular effect that increases with time. Finally, we can notice that measurements around periastron (2050 corresponds to a time of closest approach) are of prime importance. This last point will be developed further below. Changing the possible measurement window does not change these three results qualitatively.

Zoom In Zoom Out Reset image size

Figure 4 is of somewhat limited interest because "late" measurements are always favored by the adaptive scheduling tool. A more interesting result consists in identifying the regions within one orbit that are particularly sensitive to the advance of the periastron. One way to derive such a result is to run the adaptive scheduling tool multiple times by increasing the measurement window one year at a time and to save all the optimized sequences of measurements after removing the measurements that correspond to the last year of the observational window (to discard the effect of "late" measurements due to the secular effect). The result from such an analysis is presented in Figure 5 and is particularly interesting. The most sensitive measurements are astrometric measurements around closest approach. It is interesting to note that the closest approach itself is not favored by the adaptive scheduling tool (see the bottom left panel from Figure 5) because of correlations with other fitted parameters. Rather, it is important to have measurements bracketing the closest approach. The second important set of measurements are those around apastron. This can easily be interpreted since this is the location where a precession of the orbit will have the largest impact. The next set of important measurements are those of the R.A. turning point (astrometric measurements around the maximum of the R.A.). In addition, while radial velocity measurements are less crucial for the detection of the relativistic advance of periastron, the important spectroscopic measurements are all located around the closest approach and at the two turning points in the radial velocity curve.

Zoom In Zoom Out Reset image size

In order to quantify the impact of an optimized scheduling obtained using the algorithm from Section 2, we compare the gain in the relativistic (advance of periastron) S/N after one orbital period for two different scheduling strategies. More precisely, we use the current existing measurements and compare two schedulings of 32 astrometric measurements for the next 16 yr. The first scheduling consists in a "naive" sampling for which two measurements are taken every year, and the second scheduling results from our optimization such that the measurements are distributed in time following the results presented in the top left panel of Figure 5. The S/N resulting from a regular sampling is 4.1, while that resulting from the optimized scheduling is 5.9, a relative increase of 45%. This shows that, although the scheduling resulting from the algorithm presented in Section 2 is not guaranteed to be the global optimum (see the discussion in Section 3), it is nevertheless superior to standard non-optimized scheduling and it can help scientists to prioritize their measurements for a given scientific objective.

As a conclusion, the relativistic advance of the periastron is a secular effect such that a long time baseline increases the chance of its detection. On top of that, measurements at the astrometric turning points (closest approach, apastron, and R.A. turning point) are important. In the mid-term, this means that 2019 will be an important year since it corresponds to the R.A. turning point, and the years between 2022 and 2026 will be important for measuring precisely the motion of S0-2 around the apastron.

In this section, we will focus on the distance to our GC (R₀), and prioritize measurements to improve our knowledge of R₀ using the algorithm from Section 2. For this reason, we fix both parameters of general relativity to their relativistic values (i.e., to unity) so that the model contains 13 free parameters. The set of "existing measurements" is exactly the same as that used in the previous section, i.e., astrometric and radial velocity measurements of S0-2 up to 2018. The set of future possible measurements consists of astrometric and spectroscopic measurements for each year between March and September (the yearly window for which GC measurements are feasible from Earth) on a daily basis. The uncertainty considered for the future measurements is 0.2 mas for the astrometry and 20 km s⁻¹ for the radial velocity.

First, a similar behavior to that for the advance of the periastron is observed: late measurements are always favored by the adaptive scheduling tool and a long time baseline will increase our knowledge of R₀. Therefore, we use the same strategy as the one described in the previous section: we run the adaptive scheduling tool by increasing the time baseline one year at a time, record the optimized sequences of measurements but discard the observations related to the last year of the observational window. The result from such an analysis is presented in Figure 6. The estimation of R₀ relies on both astrometry and radial velocity. The important RV measurements are located at closest approach and bracket it. The important astrometric measurements are located at the R.A. turning point (maximum of the R.A. curve) and before closest approach. Of particular interest, this analysis shows that the 2019 astrometric measurements will be of prime importance for improving our knowledge about R₀.

Zoom In Zoom Out Reset image size

In order to assess the impact of the optimized scheduling, we compare the S/N of R₀ between two schedulings extending over one orbital period: one regular sampling and one optimized scheduling using the algorithm from Section 2. More precisely, we use the current existing measurements and compare two schedulings of 32 astrometric measurements and 16 RVs for the next 16 yr. The first scheduling consists of a "naive" scheduling for which two astrometric measurements and one radial velocity measurement are performed every year, and the second scheduling results from our optimization such that the measurements are distributed in time following the results presented in the top left panel of Figure 6. The S/N resulting from a regular sampling is 287 while that resulting from the optimized scheduling is 380, a relative increase of 30%. Similarly to what has been discussed in Sections 4.2 and 4.3, this shows the superiority of the scheduling obtained from the algorithm presented in Section 2 over standard non-optimized and uniform scheduling and it can help scientists to prioritize measurements for a given scientific objective.