\makeatletter
\@ifundefined{HCode}
{
\documentclass[CRPHYS, Unicode, screen, Thematic, HideArticleType]{cedram}
\newenvironment{noXML}{}{}
\newenvironment{Table}{\begin{table}}{\end{table}}
\def\thead{\noalign{\relax}\hline}
\def\endthead{\noalign{\relax}\hline}
\def\tabnote#1{\vskip4pt\parbox{.8\linewidth}{#1}}
\def\tsup#1{$^{{#1}}$}
\def\tsub#1{$_{{#1}}$}
\newcommand{\blambda}{\boldsymbol{\lambda}}
\usepackage{amssymb}
\RequirePackage{etoolbox}
\usepackage{upgreek}
\usepackage{amsmath}
\csdef{Seqnsplit}{\\}
\def\tbody{\noalign{\relax}\hline}
\def\jobid{crphys20221192}
%\graphicspath{{/tmp/\jobid_figs/web/}}
\graphicspath{{./figures/}}
\newcounter{runlevel}
\let\MakeYrStrItalic\relax
\def\refinput#1{}
\def\0{\phantom{0}}
\def\rmpi{\uppi}
\def\sfrac#1#2{{#1}/{#2}}
\def\stfrac#1#2{({#1})/{#2}}
\def\stofrac#1#2{(({#1})/{#2})}
\def\back#1{}
\def\xoverrightarrow#1{\overrightarrow{#1}}
\DOI{10.5802/crphys.144}
\datereceived{2022-12-06}
\dateaccepted{2023-02-07}
\ItHasTeXPublished
\makeatletter
\g@addto@macro{\UrlBreaks}{\UrlOrds}
\gappto{\UrlBreaks}{\UrlOrds}
}
{\documentclass[crphys]{article}
\def\CDRdoi{10.5802/crphys.144}
\def\sfrac#1#2{{#1}/{#2}}
\def\stfrac#1#2{({#1})/{#2}}
\def\stofrac#1#2{(({#1})/{#2})}
\def\href#1#2{\url{#1}}
\makeatletter
\def\selectlanguage#1{} 
\def\xoverrightarrow#1{\vec{#1}}
}
\makeatother
 
\usepackage{upgreek}
\makeatletter\def\cdr@extrapackages{upgreek}\makeatother

\dateposted{2023-04-25}
\begin{document}

\begin{noXML}
 
%\TopicEF{Exoplanets}{Exoplan\`etes}

\title{Astrometry of directly imaged exoplanets with optical
interferometry}

\alttitle{Astrom\'etrie des exoplan\`etes par interf\'erom\'etrie
optique astronomique}

\author{\firstname{Sylvestre} \lastname{Lacour}\CDRorcid{0000-0002-6948-0263}}
\address{LESIA, Observatoire de Paris, Universit\'e PSL, CNRS, 
Universit\'e Paris Cit\'e, Sorbonne Universit\'e, 5 place Jules Janssen, 
92195 Meudon, France}
\email{sylvestre.lacour@obspm.fr}

\def\thanksname{Note}
\thanks{Follows up on a conference-debate of the French Academy of
Sciences entitled ``Exoplanets: the new challenges'' held on 18 May
2021, visible via
\url{https://www.academie-sciences.fr/fr/Colloques-conferences-et-debats/exoplanetes.html}}

\thanks{\textbf{Note.}
Fait suite \`a une conf\'erence-d\'ebat de l'Acad\'emie des sciences
intitul\'ee \guillemotleft{}~Exoplan\`etes : les nouveaux
d\'efis~\guillemotright{} tenue le 18 mai 2021, visible via
\url{https://www.academie-sciences.fr/fr/Colloques-conferences-et-debats/exoplanetes.html}}

\keywords{\kwd{Astrometry}
\kwd{Exoplanets}
\kwd{Direct detection}
\kwd{Interferometry}
\kwd{High angular resolution}}

\altkeywords{\kwd{Astrom\'etrie}
\kwd{Exoplan\`etes}
\kwd{D\'etection directe}
\kwd{Interferom\'etrie}
\kwd{Haute r\'esolution angulaire}}

\begin{abstract}
Optical interferometry has always been seen as a promising but
difficult technique for astronomy. Until the year 2018, it was
excrucially limited in sensitivity to magnitudes below 10. However,
thanks to the advent of the GRAVITY instrument, fringe tracking and
dual field interferometry made it possible to observe very faint
target. This paper presents the different techniques used by optical
interferometry to perform astrometry. It shows why dual field
interferometry has become the key technique for faint companions.
Taking the exemple of the HD\,95086\,b exoplanet, we show how dual
field interferomtry can detect an exoplanet as faint as magnitude 19.5,
and how its astrometry was extracted. Use this example this paper
explains how and why an astrometric accuracy of 10 µas is
possible, and describes the remaining biases that can hinder this
measurement. Last, we conclude by presenting the orbital trajectory of
10 exoplanets in 4 systems, and conclude with the short and longer term
perspectives of the technique.
\end{abstract}

\begin{altabstract}
L'interf\'erom\'etrie aux longueurs d'ondes visibles a toujours 
\'et\'e consid\'er\'ee comme une technique prometteuse mais difficile
pour l'astronomie. Jusqu'\`a l'ann\'ee 2018, sa sensibilit\'e \'etait
limit\'ee aux objets c\'elestes les plus brillants, avec des magnitudes
inf\'erieures \`a 10. L'instrument GRAVITY, install\'e derri\`ere le VLTI
au Chili, a chang\'e cela --- gr\^ace notamment \`a sa capacit\'e de
suivre et corriger le d\'eplacement des franges et \`a son
interf\'erom\`etre double champ. Avec GRAVITY, l'on peut d\'esormais
observer des cibles tr\`es faibles. Cet article pr\'esente les
diff\'erentes fa\c cons de faire de l'astrom\'etrie avec un
interf\'erom\`etre optique. Il montre pourquoi l'interf\'erom\'etrie \`a 
double champ est devenue la technique cl\'e pour les compagnons
faibles. En prenant l'exemple de l'exoplan\`ete HD\,95086\,b, nous
montrons comment l'interf\'erom\'etrie \`a double champ peut d\'etecter
des exoplan\`etes aussi faibles que la magnitude 19,5. Nous expliquons
\'egalement comment nous obtenons des pr\'ecisions astrom\'etriques de
10 µas, et d\'ecrivons les biais restants qui peuvent entraver la
mesure. Enfin, nous terminons en pr\'esentant les orbites de 10
exoplan\`etes dans 4 syst\`emes, et concluons sur les perspectives \`a
plus long terme de cette technique.
\end{altabstract}

\maketitle

\vspace*{-1pc}

\end{noXML}

\let\rmmu\upmu

\section{Introduction}\label{sec1}
Our ability to resolve details in the sky is limited in large part by
the phenomenon of diffraction. When an astronomical object is observed
with a telescope, the Heisenberg principle implies that the direction
of the photon has a fundamental uncertainty because the photon passed
through the telescope aperture and its position has thereby been
constrained. More precisely, for a telescope of diameter $D$, the moment
of the photon, $\xoverrightarrow{p}$ can only be known to an accuracy 
of~$\sigma_{\xoverrightarrow p}$:
{\begin{equation}
\sigma_{\xoverrightarrow p} = \frac{\hbar}{D}.
\end{equation}}\unskip
Since $\xoverrightarrow p=\hbar/\lambda\xoverrightarrow \alpha$, the
uncertainty on  the direction of the photon $\xoverrightarrow \alpha$
is limited to  $\sigma_{\xoverrightarrow \alpha} = {\lambda}/{D}$.
In~practice, when observing at a  wavelength of $2.2~\upmu$m, with an 8
meter diameter telescope, the origin of a photon can not be determined
beyond a precision of 50~mas.

Several solutions have been proposed to overcome this limit. The
simplest is to accumulate many photons $N_p$ and to assume that they
all come from the same direction. This improves the astrometric
accuracy---the ability to locate a point source---by a factor
$\sqrt{N_{p}}$. Significantly more complex proposals use quantum
cloning or quantum squeezing, in order to improve the phase measurement
of a single photon beyond the diffraction 
limit~\cite{2014AandA...561A.118K}.

The technique discussed in the present article is interferometry, i.e.
the combination of light collected by several telescopes. By using the
Heisenberg principle on a photon that passes through a 
pair of telescopes, the precision on the moment of a photon is limited by
$\sigma_{\xoverrightarrow p}\cdot \xoverrightarrow B = {\hbar}$, 
where $\xoverrightarrow B$ is the vector joining the two telescopes.
Concretely, with 2 telescopes separated by 130 meters observing in the
infrared, the origin of a photon can be determined with an accuracy of
3.4~mas, along the direction of $\xoverrightarrow B$. Of course, as
with single dish telescopes, it is also possible to increase the
astrometric accuracy by $\sqrt{N_{p}}$ on a point source.

This article describes the interferometric technique as it is presently
implemented on the Very Large Telescope Interferometer (VLTI). The VLTI
is an observatory at Mount Paranal in Chile. It consists of four 8
meter diameter telescopes, plus several smaller auxiliary telescopes.
The light collected by these telescopes can be combined, i.a. with the
GRAVITY instrument that has been installed in 2018 as a VLTI new
generation  instrument~\cite{2017AandA...602A..94G}. GRAVITY includes an
original dual field combiner, which, as we will see below, enabled the
observation and astrometry of exoplanets through interferometry.

Section~\ref{sec:princ} presents the principle of an astrometric
measurement with an optical interferometer. Three techniques are
currently in use to extract astrometry: wide field observations, single
field observations and dual field observations. These are described in 
Section~\ref{sec:three}. The last technique is the only one which has
the sensitivity to directly detect exoplanets. It is described in more
details in  Section~\ref{sec:data}, with its limitations in  
Section~\ref{sec:bias}. Section~\ref{sec:results} presents direct
observations of exoplanets by our ExoGRAVITY team, and finalize this
paper with perspectives and possible improvements.

\section{Principle of astrometry with optical
interferometry}\label{sec:princ}
The purpose of astrometry is to determine the position of an
astronomical object. The Gaia observatory  measures absolute positions
of stars. It is a self-calibrated mission which yields absolute
astrometric measurements   with respect to sets of 
quasars~\cite{2016AandA...595A...1G}. With optical interferometry, all
measurements are done with respect to a physical length,
$\xoverrightarrow B$, which denotes the separation of the two telescopes
in an interferometer. $\xoverrightarrow B$ changes with time, as  Earth
rotates. It should therefore not be seen as a constant baseline, but as
a vector that rotates relative to astronomical targets.

More precisely, in an interferometer, the origin of a photon is
obtained from the difference of its optical paths through the two arms
of an interferometer. This so-called OPD is obtained by the phase
$\phi$ of the interference fringes:
{\begin{equation}
\mathrm{OPD~modulo}\ (\lambda)\ =\phi \frac{\lambda}{2\rmpi}.
\end{equation}}\unskip
The OPD varies as a function of the wavelength. For a given
measurement, at a given wavelength, the OPD is only known modulo the
wavelength. The absolute optical path through each arm of the
interferometer is not measured, but it includes the path of the light
from the star to the telescope, and from the telescope to the detector.

Assuming that the star is located at infinity, this difference in path
can be written as a scalar product between two vectors:
$\xoverrightarrow B$ the vector that joins the two telescopes, and
$\xoverrightarrow \alpha$ the position of the target in the sky:
{\begin{equation}
\mathrm{OPD}=\xoverrightarrow B \cdot \xoverrightarrow \alpha.
\end{equation}}\unskip
It is immediately observable that a phase measurement, thus the OPD,
only produces a scalar product measurement: the direction of the photon
can not be measured along the direction orthogonal to  $\xoverrightarrow
B$. If the baseline  $\xoverrightarrow B$ is  orthogonal to the
direction of the target, then we observe an OPD which is zero. The
information can however be recovered by using an interferometer with
more than 2 telescopes (with several, non-aligned  baselines) or by
using the rotation of Earth which modifies $\xoverrightarrow B$ in the
celestial reference frame.

In fact, an absolute measurement of the OPD is difficult. Therefore one
rather measures a relative of optical path difference ($\delta
\mathrm{OPD}$). This works by observing two objects in two directions
of the sky, $\xoverrightarrow \alpha$ and $\xoverrightarrow \beta$:
{\begin{equation}
\delta \mathrm{OPD}= 
\xoverrightarrow B\cdot
(\xoverrightarrow \beta-\xoverrightarrow \alpha).
\label{eq:dopd}
\end{equation}}\unskip

It might be useful to indicate some orders of magnitude: a modern
interferometer typically measures the phase of interference fringes
with an accuracy of  1\textdegree. In the near infrared, this
corresponds to an accuracy of 6~nm on the $\delta \mathrm{OPD}$. Using
an interferometer with a baseline of 130~m, this corresponds to an
angle of about $10~\upmu$as on sky, much below the angular resolution
of the interferometer  at 3.4~mas.

\section{Three techniques of astrometry in optical
interferometry}\label{sec:three}
Astrometry with an optical interferometer aims at measuring $\delta
\mathrm{OPD}$, as shown by Equation~(\ref{eq:dopd}). There are several ways
to measure $\delta \mathrm{OPD}$, which mostly differ by the hardware
of the optical interferometer. Wide angle astrometry was implemented by
the Navy Precision Optical Interferometer 
(NPOI)~\cite{1998ApJ...496..550A}. Single field astrometry was for
example used by the Sydney University Stellar Interferometer 
(SUSI)~\cite{2007MNRAS.377..415N}. Dual field interferometry is a
technique first implemented by the GRAVITY instrument.

\subsection{Wide angle astrometry}
The aim of an interferometer is to detect and record  fringes---any
photons which does not contribute to the fringe pattern, contribute
to the noise. At a given time, for a wavelength $\lambda$ and baseline
$\xoverrightarrow {B}$, an interferometric measurement yields 2
quantities: the phase and the amplitude of the fringes. Together, phase
and amplitude yield a complex value, which is called the coherent flux,
denoted $\mathbf{C}$ here.

In the case of a single object the coherent flux $\mathbf{C}$ writes:
{\begin{equation}
\mathbf{C}=
S_{\alpha}\,\exp
\left(2\rmpi j 
\frac{ \xoverrightarrow {B}\cdot \xoverrightarrow\alpha}{\lambda} 
+ j P_{\mathrm{atm}}\right)=
S_{\alpha}\,\exp
\left(2\rmpi j 
\frac{\mathrm{OPD}}{\lambda} + 
j P_{\mathrm{atm}}\right) .
\end{equation}}\unskip
$S_{\alpha}$ is the amplitude of the fringes, and of the coherent
flux. On an unresolved object, the coherent flux results from the
emission spectrum of the astronomical object, minus the atmospheric
absorption and the loss of coherence due to instrumental effects. The
atmospheric effect is represented by a chromatic piston
$P_{\mathrm{atm}}$. The position of the star in the sky,
$\xoverrightarrow\alpha$, introduces a linear slope on the phase of 
$\mathbf{C}$ as a function of $1/\lambda$. 

This slope of the phase of $\mathbf{C}$ can be theoretically used to
retrieve the OPD and therefore do an astrometric measurement.
Unfortunately, the uncertainty on $P_{\mathrm{atm}}$ makes that
impossible. The solution is for the interferometer to  point two stars,
one after the other, and to measure the phase of the fringes at two
moments, $t_1$ and $t_2$:
{\begin{equation}
\arg \mathbf{C}(t_1)-
\arg \mathbf{C}(t_2)=
2\rmpi  
\frac{\xoverrightarrow {B(t_1)} \cdot 
\xoverrightarrow \alpha - \xoverrightarrow {B(t_2)} \cdot 
\xoverrightarrow \beta }{\lambda}\ \%\ (2\rmpi).
\end{equation}}\unskip

The problem of this technique is that the atmospheric piston is not the
same at $t_1$ and $t_2$. A solution consists in averaging the piston
over time, but this can considerably impair the astrometric precision.
Another difficulty is Earth rotation, which implies that the baseline
$\xoverrightarrow B$ changes between the two measurements. This is not a
difficulty \textit{per se}, because the effect of Earth rotation can be
precisely estimated. However, the relative optical path difference
\mbox{$\xoverrightarrow {B(t_1)} \cdot \xoverrightarrow \alpha -
\xoverrightarrow {B(t_2)} \cdot \xoverrightarrow \beta$} can be huge, and
can introduce considerable biases because of errors in baseline
estimation, atmospheric dispersion, etc. This also makes it a challenge
to achieve high precision astrometry.

\subsection{Single field astrometry}
In the single field mode, or imaging mode, the interferometer detects simultaneously two objects, which are both in its field of view. In
this case, the interference fringes are superposed. The measured
quantity is a coherent flux $\mathbf{C}$ which is the sum of the 
coherent fluxes of the 2 objects:
{\begin{equation}
\mathbf{C}=
S_{\alpha}\,\exp
\left(2\rmpi j 
\frac{ \xoverrightarrow {B}\cdot 
\xoverrightarrow\alpha}{\lambda} + 
j P_{\mathrm{atm}}\right)+
S_{\beta}\,\exp
\left(2\rmpi j 
\frac{\xoverrightarrow {B}\cdot 
\xoverrightarrow\beta}{\lambda}+ 
j P_{\mathrm{atm}}\right),
\end{equation}}\unskip
or
{\begin{equation}
\mathbf{C}=
\left[S_{\alpha}\,\exp
\left(2\rmpi j 
\frac{ \xoverrightarrow {B}\cdot 
\xoverrightarrow\alpha}{\lambda} + 
j P_{\mathrm{atm}} 
\right) \right]
\left[1+
\frac{S_{\beta}}{S_{\alpha}}\,\exp
\left(2\rmpi j 
\frac{\xoverrightarrow {B}\cdot 
(\xoverrightarrow\beta - \xoverrightarrow\alpha) }
{\lambda}\right)\right],
\end{equation}}\unskip
where $S_{\alpha}$ and $S_{\beta}$ are the amplitudes of the two
coherent fluxes (from the star and planet for example).

The techniques to obtain $\delta \mathrm{OPD}= \xoverrightarrow
B\cdot(\xoverrightarrow \beta-\xoverrightarrow \alpha)$ are multiple. One can
look at the amplitude of the visibility $|\mathbf{C}|$ which changes as
a function of $\xoverrightarrow {B}/\lambda$. Another technique uses the
differential phase, i.e. the phase of  $\mathbf{C}$ as a function of
$\lambda$: one fits the modulation term
$[1+
\sfrac{S_{\beta}}{S_{\alpha}}\,\exp
(2\rmpi j
\stfrac{\xoverrightarrow {B}\cdot (\xoverrightarrow\beta -
\xoverrightarrow\alpha) }{\lambda})]$. Yet another technique
computes the closure phase, and thereby removes the phase term in
$[S_{\alpha}\,\exp
(2\rmpi j 
\stfrac{ \xoverrightarrow {B}\cdot
\xoverrightarrow\alpha}{\lambda}+ 
j P_{\mathrm{atm}}) ]$.

Single field astrometry is widely used and efficient for a number of
astronomical targets, but it is not practical for exoplanets. The main
reason is that it mixes the flux of the planet and the star. The
dynamic range is thus limited by the noise from the many stellar
photons. In addition, biases limit the accuracy of phase measurements
to half a degree. While an exoplanet, with a flux $10^{-4}$ fainter
than its star, produces a signal of $360/10^4=0.036\text{\textdegree}$.
Hence, single field astrometry doesn't reach the required precision on
phase measurements to detect planets in high contrast systems. 

\subsection{Dual field astrometry}
In  double field mode, the interferometer is in fact composed of two
interferometers: one observes the star, the other observes the planet.
The interferometer which observes the star usually acquires frames at
high frequency in order to track the fringes, and hence to cancel the
fast atmospheric perturbations. The second interferometer observes the
fainter companion and acquires frames with long integration times. This
minimizes the detector readout noise which could otherwise prevent the
observation of faint objects. In the following, parameters that relate
to the first, resp. second interferometer contain a subscript
$\mathit{FT}$, resp. $\mathrm{SCIENCE}$.

\begin{figure}
\includegraphics{fig01}
\caption{\label{fig:exmp}The principle of dual field interferometric
astrometry. Right panel: acquisition camera of GRAVITY, with the wide
field divided into two. The white arrow shows the position of the
planet. Left panel: The planet and the star are simultaneously observed
by 2 interferometers, one centered on the planet, the other on the
star. Both interferometers detect a fringe pattern. The relative
position of these two fringe patterns---also called relative optical
path difference---is related to the scalar product of the separation
between star ($\protect\xoverrightarrow{\alpha}$) and planet 
($\protect\xoverrightarrow{\beta}$) times the
baseline vector ($\protect\xoverrightarrow{B}$). 
Thus the position of the planet relative to
its star is precisely assessed from the relative location of the fringe
patterns.}
\end{figure}
 
Each interferometer has a narrow field of view, typically of the size
of the diffraction limit of a single telescope. With GRAVITY, the
interferometers use single-mode fibers with acceptance cones equivalent
to a field of view of 60~mas on the UTs. The right panel of 
Figure~\ref{fig:exmp} shows the field of view of the GRAVITY
acquisition camera, on which the 2 objects, the star and the companion,
can be seen. The acquisition field of view is of the order of 2
arcseconds with the UTs. Within this small field, the interferometer
selects an even smaller field of view, of the size of a point spread
function (${\approx}60$~mas on the UTs). Therefore, both interferometers
observe  different astronomical objects. The first interferometer,
points towards $ \xoverrightarrow\alpha$ and measures a coherent flux
equal to:
{\begin{equation}
\mathbf{C}_{\mathbf{FT}}=
S_{\alpha}\,\exp
\left(2\rmpi j 
\frac{ \xoverrightarrow {B}\cdot 
\xoverrightarrow\alpha}{\lambda}+ 
j P_{\mathrm{atm}}\right),
\end{equation}}\unskip
while the second interferometer, points towards $ \xoverrightarrow\beta$
and measures a different coherent flux:
{\begin{equation}
\mathbf{C}_{\mathbf{SCIENCE}}=
S_{\beta}\,\exp
\left(2\rmpi j 
\frac{ \xoverrightarrow {B}\cdot 
\xoverrightarrow\beta}{\lambda} +
j P_{\mathrm{atm}} +
j\mathrm{NCP} \right).
\end{equation}}\unskip
Precise metrology is necessary to measure the non-common optical path
between the two interferometers (the NCP term). The relative optical path
difference, $ \xoverrightarrow {B}\cdot (\xoverrightarrow\beta -
\xoverrightarrow\alpha)$, can be deduced from a combination of the
complex vectors: 
{\begin{equation}
\arg  
(\mathbf{C}_{\mathbf{SCIENCE}} \cdot 
\mathbf{C}_{\mathbf{FT}}^* \cdot 
\exp(-j\mathrm{NCP}  ) ) =  
2\rmpi  
\frac{\xoverrightarrow {B} \cdot
( \xoverrightarrow \beta - \xoverrightarrow \alpha) }
{\lambda}\ \%\ (2\rmpi). 
\end{equation}}\unskip
For exoplanets, however, the situation is a little bit more complicated
since the coherent flux on the science interferometer is often
contaminated by the stellar flux (see reduction technique in 
Section~\ref{sec:data}).

The direct observation of exoplanets uses the double field mode. This
mode is particularly interesting for several reasons. First, one
interferometer precisely measures the phase of the star, and uses this
to correct for atmospheric disturbances. This allows to freeze the
turbulence on the second interferometer which observes the exoplanet.
In that way, the second interferometer can make very long exposures on the
exoplanet, with integration times above 100~s. Second, the science
interferometer uses a single mode fiber which has the capacity to
filter out the stellar light. This drastically attenuates the photon
noise caused by the star on the science channel. These two effects have
allowed to observe faint exoplanets with K magnitudes as large as 19,
next to bright stars with K magnitudes down to 5.

\section{Data reduction and stellar light removal}\label{sec:data}
Despite our best efforts, on exoplanets at small angular separation,
most of the flux that enters the science interferometer is still
stellar light. This stellar light contributes to the detected signal in
two ways. The dominant contribution is incoherent light: photons with
phases that vary over milliseconds due to atmospheric perturbations.
Over the long integration times  (${>}10$~s), the fringes of such
atmospheric speckles are smeared out. This contribution is therefore
denoted ``incoherent flux''. The second contribution comes from
quasi-static speckles and is mostly due to imperfect optics: the phases
of the associated photons are constant over several seconds, and
produce coherent flux. This is a contribution from ``coherent flux''. 

The plot in  Figure~\ref{fig:contrast} traces the two contributions as
a function of the angular separation to the star. The three dashed
lines represent the expected stellar flux observed as function of the
performance of the adaptive optical (AO) correction, and as a function
of the distance to the star. This is the flux that we expect to record
in the science channel, the channel which is used to observe the
exoplanet. The doted points are actual measurements. These measurements
are observations of the exoplanetary system HD\,95086, taken during
technical time in the night of 14th February  2019.\footnote{ESO
program ID 60.A-9102(A).} The dotted blue points are the raw incoherent
flux values recorded on the science detector. They roughly follow the
prediction given by the red dashed line. The orange dots show the
amount of coherent flux, that is the stellar flux which produces
fringes, after leaking into the offset fiber of the science channel. 
Figure~\ref{fig:contrast} shows that the coherent flux from the
exoplanet is much smaller than the coherent flux from the star.

\begin{figure}
\includegraphics{fig02}
\caption{\label{fig:contrast}Incoherent (blue dots) and coherent
(orange dots) flux as a function of angular separation, normalized by
the total stellar flux, observed on HD\,95086. The dashed curves are
theoretical predictions of the incoherent flux as a function of the
performance of the adaptive optical correction. The red dot shows the
amount of coherent flux coming from the exoplanet HD\,95086\,b.}
\end{figure}
 
The exoplanet HD\,95086\,b~\cite{2013ApJ...772L..15R} was observed at a
separation of 624~mas. The contrast ratio between the planet and its
star is $1.6\times10^{-5}$: 
the planet is a factor 62,500 fainter than the star.
The flux leaked by the star into the science fiber is reduced by a
factor 500, and is thus still 125 times brighter than the exoplanet.
The coherent stellar flux is  12 times brighter than that of the
exoplanet. So the key is to disentangle the coherent energies from the
star and the exoplanet.

This is possible because the coherent energies from the star and the
exoplanet have different phases. The observing procedure actually
consists of two steps. In the first step, GRAVITY observes the star
in both interferometers:
{\begin{eqnarray}
\mathbf{C}_{\mathbf{FT}}(t_1)
&=&
S_{\mathrm{star}}\,\exp
\left(2\rmpi j 
\frac{ \xoverrightarrow {B}\cdot \xoverrightarrow\alpha}
{\lambda}\right) \cdot 
\exp(j P_{t_1} ) 
\Seqnsplit
\mathbf{C}_{\mathbf{SCIENCE}}(t_1)
&=&
S_{\mathrm{star}}\,\exp
\left(2\rmpi j 
\frac{ \xoverrightarrow {B}\cdot \xoverrightarrow\alpha}
{\lambda} \right) \cdot 
\exp(j P_{t_1} + 
j \mathrm{NCP}_{t_1}) .
\end{eqnarray}}\unskip
In a second step, the fringe tracker and science interferometer
respectively observe the star and the exoplanet:
{\begin{eqnarray}
\mathbf{C}_{\mathbf{FT}}(t_2)
&=&
S_{\mathrm{star}}\,\exp\left(2\rmpi j 
\frac{ \xoverrightarrow {B}\cdot \xoverrightarrow\alpha}{\lambda}\right) 
\cdot \exp(j
P_{t_2} )
\Seqnsplit
\mathbf{C}_{\mathbf{SCIENCE}}(t_2)
&=&
\left[S_{\mathrm{planet}}\,\exp
\left(2\rmpi j \frac{ \xoverrightarrow {B}\cdot \xoverrightarrow\beta}
{\lambda} \right) + 
\mathbf{U}_{\blambda} \cdot S_{\mathrm{star}}\,\exp
\left(2\rmpi j 
\frac{ \xoverrightarrow {B}\cdot \xoverrightarrow\alpha}{\lambda}
\right) \right] \cdot 
\exp(j P_{t_2}  +
j \mathrm{NCP}_{t_2}) .
\quad
\end{eqnarray}}\unskip
$\mathbf{U}_{\blambda}$ is the complex extinction caused by the off-axis
pointing of the single mode fiber. In the case of HD\,95086\,b in  
Figure~\ref{fig:contrast}, $|\mathbf{U}_{\blambda}|\approx1/5000$. But
still, $S_{\mathrm{planet}}\approx 1/12 \times |\mathbf{U}_{\blambda}|
\cdot S_{\mathrm{star}} $.

The goal is to derive $\xoverrightarrow{\alpha\beta} = \xoverrightarrow
\beta - \xoverrightarrow \alpha$. To do so, we assume the following:
\begin{itemize}
\item 
we know $S_{\mathrm{planet}}/S_{\mathrm{star}}$ to a constant, $\rho$,
which is the contrast ratio between the star and the planet.
\item 
$\mathbf{U}_{\blambda}$ is actually a polynomial as a function of lambda:
$\mathbf{U}_{\blambda}=\sum_{i=0}^{N} 
\mathbf{u}_{\mathbf{i},\mathbf{dit},\mathbf{baseline}}\lambda^i$. 
Where $\mathbf{u}_{\mathbf{i},\mathbf{dit},\mathbf{baseline}}$ is a
complex value which depends on the exposure time and baseline. The
typical order of the polynomial is either $N=4$ for large separations
(${>}200$~mas) or 6 for smaller separations.
\end{itemize}
So that we can estimate 
$\widetilde{\mathbf{C}}_{\mathbf{SCIENCE}}(t_2)$ from
the equation:
{\begin{eqnarray}
\widetilde{\mathbf{C}}_{\mathbf{SCIENCE}}(t_2)=
\frac{\mathbf{C}_{\mathbf{FT}}(t_2) \cdot 
\mathbf{C}_{\mathbf{SCIENCE}}(t_1) 
\exp(-j \mathrm{NCP}_{t_1})} 
{\mathbf{C}_{\mathbf{FT}}(t_1)
\exp(-j \mathrm{NCP}_{t_2}) } 
\left[ \rho\,\exp\left(2\rmpi j 
\frac{ \xoverrightarrow {B}\cdot 
\xoverrightarrow{\alpha\beta}}{\lambda} \right) +
\sum_{i=0}^{N}
(\mathbf{u}_{\mathbf{i},\mathbf{dit},\mathbf{baseline}})
\lambda^i \right].
\nonumber\\
\label{eq:CT}
\end{eqnarray}}\unskip
The least square solution can then be calculated, over the full
dataset, obtained for various wavelengths, integration times and
baselines:
{\begin{equation}
\chi^2=
\sum_{\lambda,\mathrm{dit},\mathrm{baseline}}
\left|\frac{\mathbf{C}_{\mathbf{SCIENCE}}-
\widetilde{\mathbf{C}}_{\mathbf{SCIENCE}}}
{\sigma_{\lambda,\mathrm{dit},\mathrm{baseline}}}
\right|^2 .
\label{eq:chi2}
\end{equation}}\unskip
$\mathbf{C}_{\mathbf{SCIENCE}}$ in Equation~(\ref{eq:CT}) is linear as
a function of $\rho$ and
$\mathbf{u}_{\mathbf{i},\mathbf{dit},\mathbf{baseline}}$. Hence a fast
weighted linear regression can be used to find the minimum $\chi^2$ for
a given $\xoverrightarrow{\alpha\beta}$. Therefore, the data reduction
consists in producing 2 dimensional $\chi^2$ maps as a function of the
right ascension and declination of  $\xoverrightarrow{\alpha\beta}$.  

\begin{figure}
\includegraphics{fig03}
\caption{\label{fig:chi2map}Search for the exoplanet in the HD95087\,b
data. The exoplanet is found by making a $\chi^2$ map, as a function of
$\protect\xoverrightarrow{\alpha\beta}$ and minimized over the terms $\rho$ and
$\mathbf{u}_{\mathbf{i},\mathbf{dit},\mathbf{baseline}}$. The 6 upper
maps correspond to the $\chi^2$ for a given pair of telescopes. The
lower maps correspond to the generalized $\chi^2$ map, with 2 different
ranges for $\protect\xoverrightarrow{\alpha\beta}$. The purple dots correspond
to the minima for the 8 individual exposures. The scatter of the purple
dots in the narrow $\protect\xoverrightarrow{\alpha\beta}$ is used for
computing the final error bars.}
\end{figure}

The $\chi^2$ maps obtained on the HD\,95086\,b dataset are shown in 
Figure~\ref{fig:chi2map}. The data acquired on the exoplanet is made of
3 exposures, with 16 DITs of 30~s, with the 2 polarisations. Each
polarisation was analysed separately, resulting into 6 data cubes. For
each data cubes, the position of the minimum $\chi^2$ is estimated
(purple dots in  Figure~\ref{fig:contrast}). It is interesting to see
the multiple minima in the $\chi^2$ maps, which can correspond to real
detections, but with an error of 10~mas in the detected position.
These minima reflect the fact that each baseline yields a degenerated
answer in the right ascension and declination space. The four upper
plots correspond to the $\chi^2$ maps, minimized over the linear
parameters $\rho$ and 
$\mathbf{u}_{\mathbf{i},\mathbf{dit},\mathbf{baseline}}$ for a single
baseline. We can see that for a given pair of telescopes, the
measurements are constrained mostly along one axis, modulo a $2\rmpi$
in phase. The degeneracy can be seen to be caused by the term
$\exp(2\rmpi j  \stfrac{ \xoverrightarrow {B}\cdot
\xoverrightarrow{\alpha\beta}}{\lambda})$. The uncertainty on the scalar
product, which causes the uncertainty perpendicularly to the baseline,
disappears as Earth rotates. The uncertainty on the $2\rmpi$ term along
the direction of the baseline can be reduced by increasing the width of
the spectral bandpass.

The estimation of the error bars are obtained from the scatter of the
$\chi^2$ minimum in a reduced $\chi^2$ map (purple dots in the lower
right panel of  Figure~\ref{fig:contrast}). The position of the
exoplanet relative to its star is this way estimated to be
$363.25\pm0.09$~mas in right ascension, and $-507,47\pm0.16$~mas in
declination. An extra term is calculated, which is the covariance
factor between the two directions. This value, also called Pearson
correlation coefficient, is $-0.31$ in our dataset. The anti-correlation
between right ascension and declination can be visualized by the 
$\chi^2$ valley in the lower right panel of  Figure~\ref{fig:contrast}.

\section{Astrometric biases}\label{sec:bias}

\subsection{Biases caused by the baselines}
It is crucial to have an accurate estimate of the interferometric
baseline $ \xoverrightarrow B$. Indeed, as can be seen in
Equation~(\ref{eq:dopd}) an error on $ \xoverrightarrow B$ translates
directly into a bias on the measurement of $\xoverrightarrow
\beta-\xoverrightarrow \alpha$. There are however at least 3 ways to
define the interferometric  baseline~\cite{2013ApJ...764..109W}, as
showed in  Figure~\ref{fig:baselines}:
\begin{description}
\item[$\xoverrightarrow {B_{\mathrm{WAB}}}$] 
the vector joining the pivot points of the telescopes. This baseline
does not change with the pointing of the telescopes. For this reason,
it is often called the wide angle baseline, and used for wide angle
mode observations. 
\item[$\xoverrightarrow {B_{\mathrm{IMG}}}$] 
the vector joining the telescopes' pupil centers. The definition of the
pupil as the place where all beams cross independantly of their origin,
makes this definition of the baseline independent of where the object
lies in the field of the interferometer. It is therefore used to
observe multiple objects in the field of a single interferometer, and
when doing imaging astrometry.
\item[$\xoverrightarrow {B_{\mathrm{NAB}}}$] 
the vector joining the positions of the metrological receivers. This
baselines is useful when metrology is used to measure non common path
values. It has the property of not changing whatever the path of the
metrological light inside the interferometer. This baseline is used in
dual field astrometry.
\end{description}
In a perfect interferometer the three vectors, corresponding to these
three definitions, have the same length and direction. If this is the
case, then the baseline is independent of the pointing direction of the
telescope, of the position of the object in the interferometric field,
and of the position of the optical metrology receivers. Unfortunately,
no interferometer is perfect, and the three vectors are never strictly
identical. An astrometric error budget has been established for GRAVITY
followed by a three-step approach:
\begin{enumerate}[(1)]
\item 
A campaign to measure the pivot points of all telescopes in the Earth
reference frame, using an external metrology system. The coordinates of
these pivot points in the celestial reference system are obtained using
the Earth Orientation Parameters given by the International Earth
Rotation Service (IERS).
\item 
To stabilize the telescope pupils, 4 beacons were installed on the
spiders of the 4 telescopes. These beacons send light through the
optical delay lines to the acquisition camera of the GRAVITY
instrument. The beacons are used to stabilize the pupil of the
instrument with a precision of 4~cm.
\item 
The measurement points of the metrology system consist of 4 InGaAS
photodetectors per telescope, placed on the spider arms. These
detectors sense light sent from the recombination lab and that passes
all optics of the interferometer up to M1, included. The exact
positions of these receivers were measured with respect to the pivot
points of the telescope with an accuracy below the millimeter level.
(Figure~\ref{fig:baselines}).
\end{enumerate}

\begin{figure}
\includegraphics{fig04}
\caption{\label{fig:baselines}The 3 baselines of optical interferometry
for astrometry: the wide angle baseline, the narrow angle baseline, and
the imaging baseline.}
\end{figure}

A more precise definition of these baselines, including mathematical
definitions, is presented  in~\cite{2013ApJ...764..109W}
and~\cite{2014AandA...567A..75L}. According to theoretical investigations,
and thanks to mitigation techniques, the bias introduced by
inaccuracies in baseline estimations should still allow astrometric
precisions of  $10~\upmu$as  in dual field observations.

\subsection{Other biases}
Inaccuracies in baseline estimations are the major source of bias.
However, a few other biases were identified during GRAVITY operations.
Most of them are related to non-common optical paths. They are of two
types: the NCP which are not seen via metrology, and the NCP that are
incorrectly seen through metrology.

Optical field aberrations are part of the first category. They
especially affect single field imaging observations, where the objects
are seen within the field of a single interferometer. These aberrations
introduce a field-dependent phase offset, up to several tens of
degrees. In single field observations, this NCP can be calibrated using
a phase  map~\cite{2021AandA...647A..59G}.  However, in dual field
interferometry, each target is maintained in the center of the field of
each interferometer, so this effect is not a problem.

Dispersion is also part of the first category of undetected NCP. Up to
the telescope entrance, the dispersion is caused by the atmosphere and
is mostly identical between telescopes, so it cancels out and can be
neglected. However, inside GRAVITY, if the differential delay lines
between the fringe tracker and the science are not in vacuum, they
introduce non-common-path dispersion. The GRAVITY differential delay
lines are fiber-based and introduce a lot of dispersion. The solution
used by GRAVITY is to probe the NCP with a laser which is at a
wavelength close to the astrophysical light, at  $1.9~\upmu$m. In
addition, the remaining dispersion between the K band and  $1.9~\upmu$m
is estimated by the pipeline and added to the measure from the
metrology.

High order aberrations are part of the second category of wrongly
sensed NCP. The photodetectors used for metrology sample a small part
of the pupil. High order aberrations are therefore sensed 
inaccurately, with errors up to 10\textdegree\ (60~nm). A solution
consists in averaging these aberrations over time. Other mitigation
techniques are currently developed by the GRAVITY collaboration in the
context of GRAVITY$+$. GRAVITY$+$ is an ongoing effort to enhance the sensitivity of GRAVITY~\cite{2019vltt.confE..30E}.

\section{Some results and perspectives}\label{sec:results}

\subsection{Astrometry of multi-planetary systems}
Following the first direct detection of an exoplanet by 
GRAVITY~\cite{2019AandA...623L..11G}, the ExoGRAVITY  team\footnote{ESO
program ID 1104.C-0651.} started a campaign to survey a set of directly
imaged exoplanets. At the beginning of the project, the only imaged
multi-planetary system was HR  8799~\cite{2010Natur.468.1080M}. PDS\,
70\,c followed in 2019 with a detection in H$\alpha$ by the  MUSE 
instrument~\cite{2019NatAs...3..749H}. 

Around the same time, radial velocity campaigns started to hint at new
planets that were young enough, and had long enough periods to be
imaged. Lagrange \etal\ obtained evidence of a second planet in the
$\beta$ Pictoris  system~\cite{2019NatAs...3.1135L}. With a semi major
axis of $\approx 2.7$ AUs, it is closer to the star than $\beta$
Pictoris b. GRAVITY observed it a few years later, using predictions on
its position from radial velocity 
data~\cite{2020AandA...642A..18L,2020AandA...642L...2N}.

The HD\,206893 system also happened to be very interesting. The
presence of a brown dwarf orbiting at 10 AU was already 
known~\cite{2017AandA...597L...2M}, but new HARPS data revealed the
presence of another orbiting body, which appeared to be an 
exoplanet~\cite{2019AandA...627L...9G}. By combining radial velocity
measurements, proper anomaly measurements from Gaia, and the GRAVITY
astrometry of HD 206893 B, the on-sky position of HD 206893\,c could be
predicted in  2021~\cite{2021AandA...652A..57K}. It was effectively
detected, and characterized as a $12.7\pm1$~Mjup exoplanet, a year
later by  GRAVITY~\cite{2022arXiv220804867H}.

The astrometric accuracy of these observations shows that a few
observations already constrain well the orbital parameters (Figure~\ref{fig:cohplanet}). For PDS70 b and c, two observations
separated by 6 months allowed to  derive the semi-major axis ($20\pm1$
for b and $33.2\pm2.5$ for c). The eccentricity could also be
established to $0.22\pm0.07$ and 
$0.05^{+0.05}_{-0.03}$~\cite{2021AJ....161..148W}, under the assumption
that the 2 exoplanets are co-planar. This assumption will be tested and
verified with a few additional observations over the coming years.

\begin{figure}
\includegraphics{fig05}
\caption{\label{fig:cohplanet}The four multi-planetary systems observed
by the ExoGRAVITY large  program~\cite{2020SPIE11446E..0OL}. The sizes
of the systems are quite different, with parallaxes of 8.8, 24.2, 51,
and 24.5~mas for PDS 70, HR 8799, Beta Pictoris, and HD 206893
respectively. It was possible to predict the orbital parameters
(period, eccentricity) quite accurately, with only a few GRAVITY
astrometric measurements, covering less than 3 years of orbital
motion. The orbital trajectories were adjusted using the Orbitize!
software~\cite{22}.}
\end{figure}

For the Beta Pictoris system, the 2 planets were shown to be likely in
a 7:1 mean motion resonance. The eccentricities of b and c were
constrained to respectively $0.103\pm0.003$ and $0.32\pm0.02$. As for
PDS\,70, the inner companion has a higher eccentricity. Simulations
show that the eccentricity of b can be explained by a transfer of
momentum from c. This, in turn, could mean that b migrated inward,
while preserving some of its eccentricity. It was also possible to
observe the imprint of the 3 years orbital motion of b, on the 
astrometry of c! This gave a fully independent measurement of the
dynamical mass of PDS\,70 c: $10.4\pm0.3$~Mjup compared to the prior
estimate of $8.89\pm0.75$~Mjup that includes radial velocity 
data~\cite{2021AandA...654L...2L}.

\subsection{Prospects of improvement}\label{sec:perspect}
Astrometry with optical interferometers is currently limited to 10 to 
$100~\upmu$as, depending on the brightness of the source. In the
context of the implementation of 
GRAVITY$+$~\cite{2019vltt.confE..30E}, several \mbox{ameliorations} are
currently implemented to increase the accuracy of the astrometry of
exoplanets. For faint planets the astrometry is limited by the
faintness of the signal (like HD95086\,b). Therefore, the goal of 
GRAVITY$+$ is to (i)~decrease the photon noise caused by the incoherent
stellar contamination and  (ii)~increase the planetary flux injected
into the fiber. This will be achieved thank to a new adaptive optics
correction which will reach high Strehl levels (green and orange dashed
line in  Figure~\ref{fig:contrast}). Another path taken is to improve
the fringe tracking, to increase the contrast of the fringes  and
therefore the signal. For brighter object, a lot of work is on-going
behind the scene to remove the last astrometric biases described in 
Section~\ref{sec:bias}.

Another solution for the future is to increase the angular resolution
of the optical interferometer. This can be done by observing exoplanets
at shorter wavelengths. The VLTI transmit the flux down to a wavelength
$\lambda= 1~\upmu$m, so observing at this wavelength would give a
factor 2 in accuracy compared to  $2.2~\upmu$m. But a thrilling new
project would be to build a new interferometer. This would allow to
increase significantly  $\xoverrightarrow {B}$. Such interferometer
could also be optimised for visible wavelength. It could allow to
resolve the exoplanets themselves, directly look for planetary
companions, rings \ldots 

\section*{Conflicts of interest}
The author has no conflict of interest to declare.

\section*{Acknowledgements}
Firstly, I would like to express my deep gratitude to the entire team
that contributed to the development of the GRAVITY instrument. Without
their tireless efforts and dedication, this project would not have been
possible. Secondly, I extend my sincere appreciation to the members of
the ExoGRAVITY team, whose contributions were instrumental in enabling
us to make these groundbreaking exoplanetary observations.
Additionally, I would like to thank Anthony Boccaletti and Aglae
Kellerer for their invaluable assistance in proofreading this paper.
Lastly, I am grateful for the support of the French ANR, which provided
funding for this project through the grant Project-ANR-21-CE31-0017.


\CDRGrant[ANR]{Project-ANR-21-CE31-0017}

\back{}

\bibliographystyle{crunsrt}
\bibliography{crphys20221192}
\refinput{crphys20221192-reference.tex}

\end{document}
