© S. Lacour et al. 2019

1. Introduction

GRAVITY (Gravity Collaboration 2017) is an instrument used on the Very Large Telescope Interferometer (VLTI) situated at the Cerro Paranal Observatory. It can combine the light from four telescopes. These telescopes can either be the four Auxiliary Telescopes (ATs, with a primary mirror diameter of 1.8 m) or the four Unit Telescopes (UTs, with a diameter of 8 m). The specifications of the instrument were derived from the most demanding science case, which was to observe microarcsecond displacements of the light source causing the flares of the supermassive black hole Sgr A* (Genzel et al. 2010, and references herein). Such astrometric measurements are possible with 100 m baselines (Shao & Colavita 1992; Lacour et al. 2014) and were recently demonstrated on-sky by Gravity Collaboration (2018a,b,c).

In its quiescent state, Sgr A* can become fainter than K = 18 mag. Therefore, measuring its position reliably requires an integration time on the order of minutes. To enable such long integration times, it is important to correct the atmosphere effects in real time. The higher-order atmospheric wavefront distortions are compensated for by an adaptive optics (AO) system. However, the AO systems do not sense, and therefore cannot correct, the differential phase between the telescopes. This is the role of the fringe tracker: a phase-referencing target (IRS 16C in the case of Sgr A*) is used as a guide star. In real time, the optical path differences (OPD) between each pair of telescopes are computed, and are used to control the displacement of mirrors on piezoelectric systems. This is the counterpart of the AO system, but at interferometric scale.

To push the comparison a little farther: without fringe tracking, interferometry requires short integration times and deconvolution techniques. This was the time of speckle imaging (Labeyrie 1970; Weigelt 1977), when using bispectrum and closure phases was a good but not really sensitive technique. With fringe tracking, optical interferometry enters a new age: long detector integration times (DIT up to 300 s) give access to faint sources (Kmag of 19) and to the possibility of combining spectral resolution (up to 4000) with milliarcsecond spatial resolution. This is the historical equivalent to the emergence of adaptive optics: it enables a new level of science.

Fringe tracking is not new, however. Small observatories demonstrated the concept, with PTI and CHARA (Berger & Monnier 2006). On the Keck Interferometer (Colavita et al. 2013), comparable astrophysical objectives (Woillez et al. 2012) pushed a similar development for phase referencing (Colavita et al. 2010). Previous projects also existed at the VLTI: at first, the FINITO fringe tracker (Le Bouquin et al. 2008) was used in combination with the AMBER instrument (Petrov et al. 2007). More recently, ESO developed the PRIMA fringe tracker (Delplancke et al. 2006).

However, unlike AO, fringe tracking was not yet mature, and many complex problems had to be investigated for GRAVITY. A first problem was how to deal with limited degrees of freedom (the piston actuators) while many more optical path differences are measured (Menu et al. 2012). A second problem, which does not exist in AO, is that the phase signals are known only modulo 2π. A third problem, partially addressed by the AO community (Petit et al. 2008; Poyneer & Véran 2010), is how to set a correct state space control system that optimally uses a Kalman filter to cancel the vibrations (Choquet et al. 2014). A fourth difficulty is that both group-delay (GD) and phase-delay (PD) tracking need to be used in a control system to keep the best of both. The last hurdle of the project consisted of dealing with multiple closing baselines, some of them resolved. The GRAVITY fringe tracker is now the best and most modern instrument in the field of fringe tracking for optical interferometry. Below we describe the algorithms and mechanisms.

This paper builds upon earlier works from Cassaing et al. (2008), Houairi et al. (2008), Lozi et al. (2011), Menu et al. (2012), and Choquet et al. (2014) ¹. The Menu et al. (2012) paper theoretically describes modal control of the phase delay. The Choquet et al. (2014) paper simulates the expected performance of the Kalman controller. The present paper wraps up the series by presenting the final implementation on the VLTI. Section 2 is an overview of the technical implementation of the fringe tracker and is followed by Sect. 3, where the basis of the fringe sensing is defined by the observables. The control algorithm is presented in three sections: Sect. 4 defines the operational modes, Sect. 5 presents the group-delay controller, and Sect. 6 presents the phase-delay controller. Section 7 gives examples and statistics of on-sky observations. Last, Sect. 8 concludes the paper by a discussion of possible improvements to increase the sensitivity and accuracy of the fringe tracker.

2. Overview of the fringe-tracking system

2.1. Hardware

GRAVITY is equipped with two beam combiners (Perraut et al. 2018) that perform fringe tracking and scientific observations in the K band. GRAVITY has two main operational modes. In the on-axis mode, the light of one star is split 50:50, that is, equal portions of the flux go to the fringe tracking and to the scientific channel. In off-axis mode, the field is split into two with the help of a roof mirror: one of the two objects serves as fringe-tracking reference, while the scientific channel carries out longer integrations on the typically fainter science target (Pfuhl et al. 2010, 2012, 2014).

The two beam combiners are based on silica-on-silicium integrated optics (Malbet et al. 1999), optimized for 2 μm observations (Jocou et al. 2010). The beam combination scheme is pair-wise. Each telescope pair is combined using a static ABCD phase modulation. This means that for each of the six baselines, there are four outputs corresponding to a phase shift between the two beams of 0, π/2, π, and 3π/2 radians. The total 24 outputs of the beam combiner can be seen on the real time display (RTD) of the instrument. The relevant pixels are delimited by the green rectangles in Fig. 1. Within each rectangle, the flux is dispersed in six spectral channels. The 24 lines of rectangles correspond to the 24 outputs, while the two columns of rectangles correspond to the two linear polarizations.

Fig. 1. RTD of the SAPHIRA detector. The pixels in the green boxes are read by the fringe tracker and are used for tracking. They correspond to the six baselines, four ABCD outputs, two polarizations, and six wavelength channels. The names GV1 to GV4 correspond to the input beams. The values in yellow to the right correspond to the phase shift in degrees between the different ABCD outputs.

The detector is a HgCdTe avalanche photodiode array called SAPHIRA (Finger et al. 2016). It is running at 909 Hz, 303 Hz, or 97 Hz. It can also run at either low or high gain. The high gain corresponds to a gain of γ = 7 ADU per photodetection with a typical readout noise below σ_RON = 5 ADU (≈0.7 e⁻) per pixel. The low gain does not amplify the photodetections (γ = 0.5 ADU/e⁻) and is only used for very bright targets (K magnitudes below 5 on the UTs).

The flux is processed by a first local control unit (LCU) that yields values of the observables. The LCU is an Artesyn MVME6100 using an MPC7457 PowerPC®processor (Kiekebusch et al. 2010). The data are then transmitted to a second LCU by means of a distributed memory system called reflective memory network (RMN). This second LCU processes the observables to control four tip-tilt piston mirrors on piezoelectric actuators from Physik Instrument. Each actuator has its own position sensor, driven in closed loop. The cutoff frequency of the piezoelectric delay lines is then higher than 300 Hz, with a maximum optical path delay of 60 μm (Pfuhl et al. 2014).

Real-time monitoring of the fringe tracking, including live display, is done on a separate Linux workstation connected to the LCUs through the RMN. This workstation processes the data, computes the best-fit control parameters (including the Kalman parameters), and updates the control parameters of the second LCU (Abuter et al. 2016).

2.2. Software

The two LCUs use the VxWorks operating systems. The computation is done using the so-called tools for advanced control (TAC) framework, which uses the standard C language environment. The TAC processing is triggered synchronously with the SAPHIRA, following the predefined frequency of the detector.

The first LCU computes the necessary estimators for the controller:

• Four flux values (F_i); see Sect. 3.2.

• Six phase delays (Φ_i, j); see Sect. 3.3.

• Six phase-delay variance (Var(Φ_i, j)); see Sect. 3.4.

• Six group delays (Ψ_i, j); see Sect. 3.5.

• Four closure phases (Θ_i, j, k); see Sect. 3.6.

Inside this LCU, only a few parameters can be changed: the number of DITs over which each of these quantities can be averaged. The default values are presented in Sect. 3.

The second LCU is in charge of controlling the piezo-mirrors for adequate fringe tracking. Figure 2 is a block diagram of the controller algorithm:

• The group-delay control loop, based on an integrator controller, with a direct command to the piezo-actuator (in blue).

• A feed-forward predictor, based on the action to the actuators, to increase the gain of the closed-loop system (in green).

• The phase-delay Kalman controller (in red).

• Two peripheral blocks for searching the fringes and adding a π modulation to the phase delay.

Fig. 2. Block diagram of the GRAVITY fringe-tracking controller. The gray area corresponds to the GRAVITY open-loop transfer function, and the rest corresponds to the two controllers. The group-delay integrator controller is shown in blue, the phase-delay state controller in red, and the actuators predictive model is plotted in green. The group-delay controller is the main controller: it continues to track the fringes even if the instantaneous S/N is too low for phase-delay tracking. The phase-delay state controller is a closed-loop system that determines the atmospheric perturbations ${\widehat{\mathit{X}}}_V$. A proportional controller (K) corrects for the effect of the atmosphere. The last block of the group-delay controller, the quantization step, ensures that the group-delay control signal is always a multiple of 2π: the change in OPD caused by the group-delay controller is not seen by the phase-delay controller.

The dual architecture of the controllers is made to obtain both sensitivity and accuracy. In case of high signal-to-noise ratio (S/N), the Kalman filter can determine and predict the state of both atmospheric and vibrational perturbation for the best possible correction. In case of low S/N, the Kalman filter relies on its predictive model, which in the worst case, can be as simple as a constant value. In this case, the group-delay controller is still working efficiently and provides coherence instead of fringe tracking.

The third and last software element is on the Linux workstation. It is a python script that runs every 5 s on the last 40 s of data calculated on the first LCU. It computes the parameters that are then used by the second LCU, which does the real time control. It includes the parameters for the predictive control and the best gain for the Kalman filter.

3. Estimators of observables

3.1. Visibility extraction

The real-time processing load consists mostly of matrix multiplications. The pixel-to-visibility (P2VM) principle was first used for data reduction (Tatulli et al. 2007) on AMBER. AMBER used spatial modulation for fringe coding, but the formalism was subsequently adapted to work with ABCD beam combiners (Lacour et al. 2008). Using a similar notation, the relation between the incoming electric field E_i and the outgoing electric field S_k can be expressed as

$$\begin{array}{c}\hfill S_k={\displaystyle \sum _i}{T}_{k,i}{E}_i,\end{array}$$(1)

where T_k, i is the complex transfer function from the input i of the beam combiner to its output k. In the case of GRAVITY, i_max = 4 and k_max = 24. Averaged over the DIT, the flux is then equal to the average of its instantaneous intensity:

$$\begin{array}{c}\hfill ⟨|S_k{{|}^2⟩}_{\mathrm{DIT}}={⟨{\left|{\displaystyle \sum _n}{T}_{k,i}{E}_i\right|}^2⟩}_{\mathrm{DIT}},\end{array}$$(2)

or alternatively, after decomposition:

$$\begin{array}{c}\hfill ⟨|S_k{|}^2⟩=\mathfrak{R}[{\displaystyle \sum _n}|T_{k,i}{|}^2⟨|E_i{|}^2⟩+2{\displaystyle \sum _i}{\displaystyle \sum _{j>i}}⟨T_{k,i}{T}_{k,j}^{\ast }{E}_i{E}_j^{\ast }⟩].\end{array}$$(3)

The temporal average of the electric field leads to two types of coherence losses. One depends on the optical path inside the beam combiner, the other on the spatial brightness distribution of the astrophysical object. The first, C_k, i, j, is intrinsic to the device and has to be calibrated. The second, V_i, j, is the reason why we built the interferometer. The relation with the mean electric field intensity is then approximated by the following equation:

$$\begin{array}{c}\hfill ⟨T_{k,i}{T}_{k,j}^{\ast }{E}_i{E}_j^{\ast }⟩=|T_{k,i}{T}_{k,j}^{\ast }|C_{k,i,j}\sqrt{⟨|E_i{|}^2⟩⟨|E_j{|}^2⟩}{V}_{i,j}.\end{array}$$(4)

Hence, Eq. (3) can be written as a matrix product:

$$\begin{array}{c}\hfill \left(\begin{array}{c}⟨|S_1{|}^2⟩\\ ⋮\\ ⟨|S_{24}{|}^2⟩\end{array}\right)=\mathfrak{R}[\mathsf{V2PM}·\left(\begin{array}{c}⟨|E_1{|}^2⟩\\ ⟨|E_2{|}^2⟩\\ ⟨|E_3{|}^2⟩\\ ⟨|E_4{|}^2⟩\\ \sqrt{⟨|E_1{|}^2⟩⟨|E_2{|}^2⟩}{V}_{1,2}\\ \sqrt{⟨|E_1{|}^2⟩⟨|E_3{|}^2⟩}{V}_{1,3}\\ \sqrt{⟨|E_1{|}^2⟩⟨|E_4{|}^2⟩}{V}_{1,4}\\ \sqrt{⟨|E_2{|}^2⟩⟨|E_3{|}^2⟩}{V}_{2,3}\\ \sqrt{⟨|E_2{|}^2⟩⟨|E_4{|}^2⟩}{V}_{2,4}\\ \sqrt{⟨|E_3{|}^2⟩⟨|E_4{|}^2⟩}{V}_{3,4}\end{array}\right)],\end{array}$$(5)

where, using ^⊤ as the transpose operator:

$$\begin{array}{c}\hfill {\mathsf{V2PM}}^{\top }=\left(\begin{array}{ccc}|T_{1,1}{|}^2& \cdots & |T_{24,1}{|}^2\\ ⋮& \ddots & ⋮\\ |T_{1,4}{|}^2& \cdots & |T_{24,4}{|}^2\\ |T_{1,1}{T}_{1,2}^{\ast }|C_{1,1,2}& \cdots & |T_{24,1}{T}_{24,2}^{\ast }|C_{24,1,2}\\ ⋮& \ddots & ⋮\\ |T_{1,3}{T}_{1,4}^{\ast }|C_{1,3,4}& \cdots & |T_{24,3}{T}_{24,4}^{\ast }|C_{24,3,4}\end{array}\right).\end{array}$$(6)

Everything related to the transfer function of the instrument is inside the 10 columns and 24 rows of the 𝖵2𝖯𝖬 matrix. The 𝖵2𝖯𝖬 is calibrated during daytime on the internal source of the instrument. It is regularly computed for verification, but has been proven to be very stable over several months. It is part of the calibration files that are needed by the first LCU (Sect. 2.2).

The 𝖯2𝖵𝖬 is the pseudo-inverse matrix of V2PM. It can be obtained by splitting the 𝖵2𝖯𝖬 into real and imaginary parts, and taking the inverse using a singular-value decomposition (SVD). The 𝖯2𝖵𝖬 is used to retrieve the astrophysical information from the flux observed on the pixels. This information consists of the flux F and the coherent flux Γ. Both are computed from the pixel flux q_k, λ and 𝖯2𝖵𝖬 according to Eq. (7):

$$\begin{array}{c}\hfill \left(\begin{array}{c}{F}_{1,\lambda }\\ ⋮\\ F_{4,\lambda }\\ \mathfrak{R}[{\mathrm{\Gamma }}_{1,2,\lambda }]\\ ⋮\\ \mathfrak{R}[{\mathrm{\Gamma }}_{3,4,\lambda }]\\ \mathfrak{I}[{\mathrm{\Gamma }}_{1,2,\lambda }]\\ ⋮\\ \mathfrak{I}[{\mathrm{\Gamma }}_{3,4,\lambda }]\end{array}\right)={\mathsf{P2VM}}_{\lambda }·\left(\begin{array}{c}{q}_{1,\lambda }\\ ⋮\\ q_{24,\lambda }\end{array}\right).\end{array}$$(7)

In the above equation, we used a slightly different nomenclature, which is used hereafter. q_k, λ =< |S_k, λ|² > is the number of photons observed over one DIT of the fringe tracker on a given pixel. F_i, λ =< |E_i, λ|² > is the energy per DIT of the incoming beam i at wavelength λ. Last, the complex coherent flux Γ_i, j, λ corresponds to the visibility before normalization by the flux. It is obtained from the flux F_i, λ and the visibilities V_i, j, λ:

$$\begin{array}{c}\hfill {\mathrm{\Gamma }}_{i,j,\lambda }=\sqrt{F_{i,\lambda }{F}_{j,\lambda }}\mathfrak{R}[V_{i,j,\lambda }]+i\sqrt{F_{i,\lambda }{F}_{j,\lambda }}\mathfrak{I}[V_{i,j,\lambda }],\end{array}$$(8)

using $\sqrt{F_{i,\lambda }{F}_{j,\lambda }}=\sqrt{⟨|E_i{|}^2⟩⟨|E_j{|}^2⟩}$. All these values are computed in real time for each of the six wavelengths in the K band. In case of split polarization (when the Wollaston is inserted), the calculation is also done independently for both polarizations.

3.2. Flux estimator

The flux F_i is the value extracted after each DIT from Eq. (7), summed over the N_λ = 6 spectral channels:

$$\begin{array}{c}\hfill F_i={\displaystyle \sum _{\lambda =1}^{N_{\lambda }}}{F}_{i,\lambda },\end{array}$$(9)

where i corresponds to the input beam number.

3.3. Phase delay estimator

The phase delay Φ_i, j is derived from the complex coherent flux, but after a first step to correct for the phase curvature caused by the dispersion:

$$\begin{array}{c}\hfill {\mathrm{\Gamma }}_{i,j,\lambda }^{\prime }={\mathrm{\Gamma }}_{i,j,\lambda }exp\left(iD{[1-\frac{2.2\mu \mathrm{m}}{\lambda }]}^2\right),\end{array}$$(10)

which is only a first-order approximation of the dispersion. It is caused both by the atmosphere and by the fibered differential delay lines (FDDL)². Therefore, the corrective term D is a time-variable parameter that depends on the position of the star and the position of the FDDLs. The phase delay is then extracted by coherent addition of the six spectral channels:

$$\begin{array}{c}\hfill {\mathrm{\Phi }}_{i,j}=arg\left({\displaystyle \sum _{\lambda =1}^{N_{\lambda }}}{\mathrm{\Gamma }}_{i,j,\lambda }^{\prime }\right).\end{array}$$(11)

It is worth noting that the phase delay Φ_i, j is wrapped: it lies between −π and π. No unwrapping effort is made at this stage.

3.4. Phase variance estimator

Computing the S/N of the fringes is essential for the fringe tracker. It ensures that the controller does not track on noise. It is also needed for the state machine to know if it has the fringes locked or if it must start looking for the fringes elsewhere. The S/N is calculated from the variance of the phase delay: $S/N_{i,j}=1/\sqrt{\mathrm{Var}({\mathrm{\Phi }}_{i,j})}$. For each DIT, from the photon and background noise, the variance on each pixel is estimated by the relation:

$$\begin{array}{c}\hfill \mathrm{Var}(q_{i,\lambda })={\sigma }_{\mathrm{sky}}^2+\gamma (q_{i,\lambda }-q_{i,\lambda ,\mathrm{sky}}),\end{array}$$(12)

where σ_sky and q_{i, λ, sky} are the noise and flux, respectively, observed during sky observations. The term γ is the detector gain in ADU per detected photon. The covariance matrix of the real and imaginary part of the Γ terms can be obtained from the 𝖯2𝖵𝖬:

$$\begin{array}{c}\hfill {\mathrm{\Sigma }}_{\mathrm{\Gamma }}={\mathsf{P2VM}}_{\lambda }·\left(\begin{array}{c}\mathrm{Var}(q_{1,\lambda })\\ ⋮\\ \mathrm{Var}(q_{24,\lambda })\end{array}\right)·{\mathsf{P2VM}}_{\lambda }^{\top },\end{array}$$(13)

where ^⊤ is the transpose operator. To save processing time, only the diagonal values of the variance matrix are calculated. They correspond to the variance of the real and imaginary parts of Γ_i, j, λ. This assumes that the covariance between the real and imaginary parts is negligible (a good assumption for an ABCD with phase shift of 0, π/2, π,  and 3π/2). In the end, for simplicity, we estimated the variance of the phase by the following equation:

$$\begin{array}{c}\hfill \mathrm{Var}({\mathrm{\Phi }}_{i,j})=\frac{{\displaystyle {\displaystyle \sum _{\lambda }}{⟨\mathrm{Var}(\mathfrak{R}{\mathrm{\Gamma }}_{i,j,\lambda })+\mathrm{Var}(\mathfrak{I}{\mathrm{\Gamma }}_{i,j,\lambda })⟩}_{5\mathrm{DIT}}}}{{\displaystyle 2{\left|{\displaystyle \sum _{\lambda }}{⟨{\mathrm{\Gamma }}_{i,j,\lambda }^{\prime }⟩}_{5\mathrm{DIT}}\right|}^2}},\end{array}$$(14)

which is the variance of the amplitude of the coherent flux averaged over five DIT. The five-DIT average is a way to increase the precision of the calculation. However, this is done at the expense of accuracy: the coherent averaging of the complex coherent flux can add a negative bias to the phase variance estimator.

3.5. Group delay estimator

The group delay, Ψ_i, j, is also obtained from the complex coherent flux. Because it consists of a differential measure of the phase as a function of wavelength, this estimator is more noisy than the phase delay (Lawson et al. 2000). To increase its S/N, for each one of the N_λ = 6 spectral channels, the Γ_i, j, λ is first corrected for dispersion, cophased, and averaged over 40 DITs. The result is then used to derive the group delay by multiplying the phasor of consecutive spectral channels:

$$\begin{array}{cc}\hfill {\mathrm{\Gamma }}_{i,j,\lambda }^{″}& ={⟨{\mathrm{\Gamma }}_{i,j,\lambda }^{\prime }exp(-i{\mathrm{\Phi }}_{i,j})⟩}_{40\mathrm{DIT}}\hfill \end{array}$$(15)

$$\begin{array}{cc}\hfill {\mathrm{\Psi }}_{i,j}& =arg\left({\displaystyle \sum _{\lambda =1}^{N_{\lambda }-1}}{\mathrm{\Gamma }}_{i,j,\lambda +1}^{″}{\mathrm{\Gamma }}_{i,j,\lambda }^{″\ast }\right).\hfill \end{array}$$(16)

As for the phase delay, the group delay is estimated modulo 2π. In terms of optical path, however, Ψ_i, j corresponds to a value R times smaller than Φ_i, j, with R = 23, the spectral resolution of the GRAVITY fringe tracker. This is explicit in open-loop operation where the phase delay and group delay are compared for a response to top-hat piezo commands. Because both estimators wrap at 2π, this means that the estimator is valid over a long range equal to 23 times the wavelength. This is the main advantage of the group-delay estimator: to be able to find fringes far away from the central white-light fringe (the fringe of highest contrast). However, because it uses individual spectral channels, the group delay calculated on a single DIT would be extremely noisy. Thus the 40 DIT summation is a way to increase the S/N of the group delay, at the cost of losing response time. In the end, the group-delay estimator is a reliable, but slow, estimator of the optical path difference.

3.6. Closure-phase estimator

The closure phases, Θ_i, j, k, are calculated on all four triangles from the coherent flux. Before taking the argument, the bispectra are averaged over 300 DIT (corresponding to 330 ms at the fastest 909 Hz sampling rate). Closure phases are estimated from the phase delay:

$$\begin{array}{c}\hfill {\mathrm{\Theta }}_{i,j,k}^{\mathrm{PD}}=arg\left({⟨{\displaystyle \sum _{\lambda =1}^{N_{\lambda }}}{\mathrm{\Gamma }}_{i,j,\lambda }^{\prime }{\displaystyle \sum _{\lambda =1}^{N_{\lambda }}}{\mathrm{\Gamma }}_{j,k,\lambda }^{\prime }{\displaystyle \sum _{\lambda =1}^{N_{\lambda }}}{\mathrm{\Gamma }}_{i,k,\lambda }^{\prime \ast }⟩}_{300\mathrm{DIT}}\right),\end{array}$$(17)

but also from the group delay:

$$\begin{array}{c}\hfill {\mathrm{\Theta }}_{i,j,k}^{\mathrm{GD}}=arg\left({⟨{\displaystyle \sum _{\lambda =1}^{N_{\lambda }}}{\mathrm{\Gamma }}_{i,j,\lambda }^{″}{\displaystyle \sum _{\lambda =1}^{N_{\lambda }}}{\mathrm{\Gamma }}_{j,k,\lambda }^{″}{\displaystyle \sum _{\lambda =1}^{N_{\lambda }}}{\mathrm{\Gamma }}_{i,k,\lambda }^{″\ast }⟩}_{300\mathrm{DIT}}\right).\end{array}$$(18)

They are observables modulo 2π, no unwrapping is intended.

4. State spaces, projections, and state machine

4.1. OPD-state space

A main difficulty for the fringe tracker consists of dealing with the different dimensions of the vectors involved. The number of phase observables is six. The number of delay lines is four. Last, the number of degrees of freedom is three. In Menu et al. (2012), we proposed an ℝ³ modal-state space orthogonal to the piston space. However, as also mentioned in Menu et al. (2012), this modal control has an important drawback: it cannot work in a degraded mode where one or more telescopes are missing. Instead, the GRAVITY controller uses the OPD-state space. In some sense, the implemented fringe tracker is a downgraded version of the state controller proposed in Menu et al. (2012). However, it facilitates managing flux drop-out as well as working with a reduced number of baselines.

4.2. Reference vectors

For the system to work properly, the OPD-state space must be colinear to the piston space. However, because the astronomical object is not necessarily a point source, the closure phases are not necessarily zero. Therefore, the OPD component orthogonal to the piston space must be removed from the measurement. This is done by subtracting a reference position, or set point, which is computed from the closure-phase estimators ${\mathrm{\Theta }}_{i,j,k}^{\mathrm{PD}}$ and ${\mathrm{\Theta }}_{i,j,k}^{\mathrm{GD}}$. Then, the error terms, that is, the differences between the measured OPD and the set points, are colinear to the piston space.

However, the devil is in the details. There are four closure phases, and only three can be used. The noisiest closure phase is therefore discarded. The three other closure phases are applied on the three edges of the triangle forming the noisiest closure phase. The reference vector is therefore defined as:

$$\begin{array}{c}\hfill {\mathbf{Ref}}_{\mathrm{\Phi }}=\left(\begin{array}{c}{\mathrm{\Theta }}_{1,2,4}^{\mathrm{PD}}\\ {\mathrm{\Theta }}_{1,3,4}^{\mathrm{PD}}\\ 0\\ {\mathrm{\Theta }}_{2,3,4}^{\mathrm{PD}}\\ 0\\ 0\end{array}\right)\mathrm{or}\left(\begin{array}{c}{\mathrm{\Theta }}_{1,2,3}^{\mathrm{PD}}\\ 0\\ -{\mathrm{\Theta }}_{1,3,4}^{\mathrm{PD}}\\ 0\\ -{\mathrm{\Theta }}_{2,3,4}^{\mathrm{PD}}\\ 0\end{array}\right)\mathrm{or}\left(\begin{array}{c}0\\ -{\mathrm{\Theta }}_{1,2,3}^{\mathrm{PD}}\\ -{\mathrm{\Theta }}_{1,2,4}^{\mathrm{PD}}\\ 0\\ 0\\ {\mathrm{\Theta }}_{2,3,4}^{\mathrm{PD}}\end{array}\right)\mathrm{or}\left(\begin{array}{c}0\\ 0\\ 0\\ {\mathrm{\Theta }}_{1,2,3}^{\mathrm{PD}}\\ {\mathrm{\Theta }}_{1,2,4}^{\mathrm{PD}}\\ {\mathrm{\Theta }}_{1,3,4}^{\mathrm{PD}}\end{array}\right),\end{array}$$(19)

depending on which triangle has the lowest S/N, from left to right, the 123, 124, 134, or 234 triangle. The reference values for the group delay are calculated similarly:

$$\begin{array}{c}\hfill {\mathbf{Ref}}_{\mathbf{\Psi }}=\left(\begin{array}{c}{\mathrm{\Theta }}_{1,2,4}^{\mathrm{GD}}\\ {\mathrm{\Theta }}_{1,3,4}^{\mathrm{GD}}\\ 0\\ {\mathrm{\Theta }}_{2,3,4}^{\mathrm{GD}}\\ 0\\ 0\end{array}\right)\mathrm{or}\left(\begin{array}{c}{\mathrm{\Theta }}_{1,2,3}^{\mathrm{GD}}\\ 0\\ -{\mathrm{\Theta }}_{1,3,4}^{\mathrm{GD}}\\ 0\\ -{\mathrm{\Theta }}_{2,3,4}^{\mathrm{GD}}\\ 0\end{array}\right)\mathrm{or}\left(\begin{array}{c}0\\ -{\mathrm{\Theta }}_{1,2,3}^{\mathrm{GD}}\\ -{\mathrm{\Theta }}_{1,2,4}^{\mathrm{GD}}\\ 0\\ 0\\ {\mathrm{\Theta }}_{2,3,4}^{\mathrm{GD}}\end{array}\right)\mathrm{or}\left(\begin{array}{c}0\\ 0\\ 0\\ {\mathrm{\Theta }}_{1,2,3}^{\mathrm{GD}}\\ {\mathrm{\Theta }}_{1,2,4}^{\mathrm{GD}}\\ {\mathrm{\Theta }}_{1,3,4}^{\mathrm{GD}}\end{array}\right).\end{array}$$(20)

The closure-phase changes as a function of time, causing the reference position to adapt to any change in the phase closures. The closure phase is smoothed over a long enough time (300 DIT) to avoid adding additional noise. However, the choice of which triangle is the noisiest is made only once between each scientific frame to avoid sharp jumps in the reference vector over the integration time of the science detector.

This reference scheme works most of the time. However, problems arise in two specific instances. First, when the object is so highly resolved that the used closure phases contain a baseline with zero visibility. If that happens, an undefined reference value is applied to a perfectly sane baseline and the system can diverge. Second, if the fringe tracker is tracking on two unconnected baselines (e.g., between telescopes 1–2 and 3–4), the closure phases are undefined, and using their values would mean losing one of the two locked baselines. To resolve these two problems, the closure phases are modified as follows: if one of the baselines of any of ij, jk, or ik have an S/N below the value $S/N_{\mathrm{threshold}}^{\mathrm{GD}}$, then ${\mathrm{\Theta }}_{i,j,k}^{\mathrm{PD}}$ used in Eq. (19) is a fixed value, and ${\mathrm{\Theta }}_{i,j,k}^{\mathrm{GD}}=0$ is used in Eq. (20). The difference in treatment between the group and phase delay arises because the default group-delay tracking shall be zero, while the default phase-delay tracking can be any constant value.

4.3. Transfer matrices

After the reference values are subtracted from the OPD, we can freely project the data in piston-state space as well as back to the OPD-state space. Hereafter, we use the same nomenclature as in Menu et al. (2012). P corresponds to the four-dimension piston-state vector, while OPD corresponds to the six-dimension OPD-state vectors. The matrix M corresponds to the conversion between piston and optical path difference:

$$\begin{array}{c}\hfill \mathbf{\text{OPD}}=\mathsf{M}\mathbf{\text{P}},\end{array}$$(21)

where

$$\begin{array}{c}\hfill \mathsf{M}={\left(\begin{array}{cccccc}-1& -1& -1& 0& 0& 0\\ 1& 0& 0& -1& -1& 0\\ 0& 1& 0& 1& 0& -1\\ 0& 0& 1& 0& 1& 1\end{array}\right)}^{\top }.\end{array}$$(22)

The conversion OPD → P is ill-constrained, however: the rank of matrix 𝖬 is 3, not 4. This is because the global piston cannot be obtained from the differences in the optical path. Nevertheless, we can define a pseudo-inverse matrix:

$$\begin{array}{c}\hfill {\mathsf{M}}^{†}=\frac{1}{4}\left(\begin{array}{cccccc}-1& -1& -1& 0& 0& 0\\ 1& 0& 0& -1& -1& 0\\ 0& 1& 0& 1& 0& -1\\ 0& 0& 1& 0& 1& 1\end{array}\right).\end{array}$$(23)

where ^† denotes the pseudo-inverse operator.

4.4. Thresholds and S/N management

The system uses two distinct thresholds. A first threshold, $S/N_{\mathrm{threshold}}^{\mathrm{GD}}$, disables the baselines where the S/N is too low to make the baseline useful. This is the case when the fringes are not yet found, when the fringes are suddenly lost, or when the astronomical target is so highly resolved that the spatial coherence is close to zero. To detect these events, the $S/N_{\mathrm{threshold}}^{\mathrm{GD}}$ is compared to a moving average of the phase-delay variance estimator: ⟨Var(Φ_i, j)⟩_40 DIT. This group-delay threshold must be adapted to the target: it mustbe high enough to ensure that the fringe tracker does not track on side-lobes, but low enough to detect the fringes.

A second threshold, $S/N_{\mathrm{threshold}}^{\mathrm{PD}}$, is used solely for the phase-tracking controller. It disables the tracking on a given telescope in case of rapid S/N drop-off. The main target of this threshold is to be able to catch a drop in flux injection (caused by external tip-tilt), with the expectation that the S/N will increase soon hereafter. Typically, $S/N_{\mathrm{threshold}}^{\mathrm{PD}}=1.5$ and $S/N_{\mathrm{threshold}}^{\mathrm{PD}}<S/N_{\mathrm{threshold}}^{\mathrm{GD}}$.

Hence, we defined two matrices 𝕀 that convert OPD space error into a space of the same dimension colinear to the OPD space. They read

$$\begin{array}{cc}\hfill {\mathbb{I}}_{\mathrm{GD}}^6& =\mathsf{M}{({\mathsf{M}}^{\top }\mathsf{W}\mathsf{M})}^{{†}_{\mathrm{GD}}}{\mathsf{M}}^{\top }\mathsf{W}\hfill \end{array}$$(24)

$$\begin{array}{cc}\hfill {\mathbb{I}}_{\mathrm{PD}}^6& =\mathsf{M}{({\mathsf{M}}^{\top }\mathsf{W}\mathsf{M})}^{{†}_{\mathrm{PD}}}{\mathsf{M}}^{\top }\mathsf{W},\hfill \end{array}$$(25)

where the only difference is on the pseudo-inverse operator ^†. In both equations, the 6 × 6 weighting matrix 𝖶 distributes the weights among the different OPDs:

$$\begin{array}{c}\hfill \mathsf{W}=\mathrm{diag}\left(\begin{array}{cccccc}{w}_{1,2}& w_{1,3}& w_{1,4}& w_{2,3}& w_{2,4}& w_{3,4}\end{array}\right),\end{array}$$(26)

where

$$\begin{array}{c}\hfill w_{i,j}=\{\begin{array}{cc}0\hfill & \mathrm{if}1/{⟨\mathrm{Var}({\mathrm{\Phi }}_{i,j})⟩}_{40\mathrm{DIT}}<{|S/N_{\mathrm{threshold}}^{\mathrm{GD}}|}^2\hfill \\ 1/\mathrm{Var}({\mathrm{\Phi }}_{i,j})\hfill & \mathrm{otherwise}\hfill \end{array}.\end{array}$$(27)

This step is important to remove the risk of tracking on noise: if the variance of the phase reaches this threshold, the pseudo-inverse discards that baseline in its calculation. The pseudo-inversion is done using a SVD:

$$\begin{array}{c}\hfill {\mathsf{M}}^{\top }\mathsf{W}\mathsf{M}=\mathsf{U}\mathsf{S}\mathsf{V}.\end{array}$$(28)

The idea behind this decomposition is that 𝖴 and 𝖵 are two invertible orthonormal matrices (the left-singular and right-singular eigenvectors, respectively): 𝖨 = 𝖴𝖴^⊤ and I = 𝖵𝖵^⊤ , where 𝖨 is the identity matrix. 𝖲 is a square diagonal matrix where the values on the diagonal correspond to the square root of the eigenvalues: 𝖲 = diag(s₁, s₂, s₃, 0). For a four-telescope operation, three eigenvalues are non-zero. The number of non-zero eigenvalues decreases when the system cannot track all telescopes.

The pseudo-inverse of matrix S is calculated differently for the group-delay and phase-delay control loop. For the group-delay control loop, we have ${\mathsf{S}}^{{†}_{\text{GD}}}=\text{diag}(s_1^{{†}_{\text{GD}}},s_2^{{†}_{\text{GD}}},s_3^{{†}_{\text{GD}}},0)$, where

$$\begin{array}{c}\hfill s_i^{{†}_{\mathrm{GD}}}=\{\begin{array}{cc}1/s_i\hfill & \mathrm{if}{s}_i>0\hfill \\ 0\hfill & \mathrm{otherwise}\hfill \end{array}.\end{array}$$(29)

For the phase-delay control loop, we have ${\mathsf{S}}^{{†}_{\text{PD}}}=\text{diag} (s_1^{{†}_{\text{PD}}},s_2^{{†}_{\text{PD}}},s_3^{{†}_{\text{PD}}},0)$, where

$$\begin{array}{c}\hfill s_i^{{†}_{\mathrm{PD}}}=\{\begin{array}{cc}1/s_i\hfill & \mathrm{if}{s}_i>{|S/N_{\mathrm{threshold}}^{\mathrm{PD}}|}^2\hfill \\ s_i/{|S/N_{\mathrm{threshold}}^{\mathrm{PD}}|}^4\hfill & \mathrm{otherwise}\hfill \end{array}.\end{array}$$(30)

As a result, the two matrices 𝕀 are calculated for each DIT from a new SVD and the equations

$$\begin{array}{cc}\hfill {\mathbb{I}}_{\mathrm{GD}}^6& =\mathsf{M}{\mathsf{V}}^{\top }{\mathsf{S}}^{{†}_{\mathrm{GD}}}{\mathsf{U}}^{\top }{\mathsf{M}}^{\top }\mathsf{W}\hfill \end{array}$$(31)

$$\begin{array}{cc}\hfill {\mathbb{I}}_{\mathrm{PD}}^6& =\mathsf{M}{\mathsf{V}}^{\top }{\mathsf{S}}^{{†}_{\mathrm{PD}}}{\mathsf{U}}^{\top }{\mathsf{M}}^{\top }\mathsf{W},\hfill \end{array}$$(32)

where the difference between the two is that the eigenvalues in ${\mathbb{I}}_{\mathrm{PD}}^6$ are weighted down when they are below a value equal to $|S/N_{\mathrm{threshold}}^{\mathrm{PD}}{|}^2$.

4.5. State machine

The rank of matrix ${\mathbb{I}}_{\mathrm{GD}}^6$ (the number of non-zero eigenvalues) drives the decision-making of the state machine. For a four-telescope operation, a rank of 3 means that the position of the delay lines on all the telescopes is constrained. When only three telescopes are tracked, the rank is 2. The rank is 1 for two linked telescopes.

The state machine (Fig. 3) therefore has only three states: IDLE, SEARCHING, and TRACKING. When the operator starts the fringe tracker, it switches to the state SEARCHING. As soon as the rank of the ${\mathbb{I}}_{\mathrm{GD}}^6$ matrix is 3, the fringe tracker transitions to state TRACKING. When the rank of the matrix decreases and remains low for a period of 1 s or more, then the system automatically switches back to SEARCHING mode. In both states, the group-delay and phase-delay controllers are running. This means that whether in SEARCHING or TRACKING state, the system still tracks the fringes on the baselines with sufficient S/N.

Fig. 3. State machine of the fringe tracker. There are only three states: IDLE, SEARCHING, and TRACKING. The blue transitions are commands from the operator. The gray transitions are automatic decisions. Switching from searching to tracking only depends on the rank of the ${\mathbb{I}}_{\mathrm{GD}}^4$ matrix. If GRAVITY is running in a degraded mode, the TRACKING transition can happen for a rank of 2 or even 1.

5. Group delay tracking

5.1. Group-delay block diagram

Figure 2 represents the block diagram of the full control architecture. Within this design, the prevalent control loop is the group-delay loop (in blue). It has to be very reliable: the signal on the feedback signal Ψ is therefore enhanced by averaging over 40 DITs in Eq. (15). The control algorithm of the group delay is presented in Fig. 4. It generates a control signal u_GD whose unit is the radian of the phase delay. It is made of seven distinct blocks that are described below.

Fig. 4. Block diagram of the group-delay controller. The controller is centered on an integrator (block 5), but also includes several features to work in conjunction with the phase-delay controller.

The measured group delay is a vector defined as

$$\begin{array}{c}\hfill {\mathbf{\Psi }}_n={\left(\begin{array}{cccccc}{\mathrm{\Psi }}_{1,2}& {\mathrm{\Psi }}_{1,3}& {\mathrm{\Psi }}_{1,4}& {\mathrm{\Psi }}_{2,3}& {\mathrm{\Psi }}_{2,4}& {\mathrm{\Psi }}_{3,4}\end{array}\right)}_n.\end{array}$$(33)

The first block is a comparator that extracts the error between a reference vector (Ref_Ψ) and the measured group delay. The logic would be to have all six setpoints equal to zero to track on the white-light fringes. However, as explained in Sect. 4.2, this is not possible in the presence of non-zero group-delay closure phase Θ^GD. The use of a setpoint vector as defined by Eq. (20) therefore ensures that all baselines track the fringes with the highest contrast that do not contradict each other:

$$\begin{array}{c}\hfill {\mathit{\epsilon }}_{\mathbf{\Psi },n}={\mathbf{Ref}}_{\mathbf{\Psi }}-{\mathbf{\Psi }}_n.\end{array}$$(34)

The second block is there because the phase measurement is only known modulo 2π. Because ε_Ψ, n can have any value, this block adds or subtracts an integer number of 2π to ensure that the error is between −π and π.

The third block uses the ${\mathbb{I}}_{\mathrm{GD}}^6$ matrix to weight the errors between the different baselines. Following that matrix, the OPD error vector is now strictly, in a mathematical sense, colinear to the piston space. Moreover, the error on any baseline with no fringes is either estimated from other baselines or set to zero. After the third block, the error group-delay vector is now

$$\begin{array}{c}\hfill {\mathit{\epsilon }{}_{\mathbf{\Psi }}^{\mathbf{\prime }}}_{,n}={\mathbb{I}}_{\mathrm{GD}}^6([{\mathbf{Ref}}_{\mathrm{\Psi }}-{\mathbf{\Psi }}_n]\%2\pi ),\end{array}$$(35)

where the percent sign corresponds to the modulo function.

The fourth block is a threshold function that quenches the gain of the control loop if the absolute value of ε′_Ψ,n is lower than 2π/R. This value corresponds to an OPD of one wavelength, meaning the group-delay controller cannot converge on an accuracy better than 2.2 μm. This is necessary to let the phase-delay controller track within a fringe:

$$\begin{array}{c}\hfill {\mathit{\epsilon }{}_{\mathbf{\Psi }}^{\mathbf{″}}}_{,n}=\{\begin{array}{cc}{\mathit{\epsilon }{}_{\mathbf{\Psi }}^{\mathbf{\prime }}}_{,n}-\pi /R\hfill & \mathrm{if}{\mathit{\epsilon }{}_{\mathbf{\Psi }}^{\mathbf{\prime }}}_{,n}>\pi /R\hfill \\ {\mathit{\epsilon }{}_{\mathbf{\Psi }}^{\mathbf{\prime }}}_{,n}+\pi /R\hfill & \mathrm{if}{\mathit{\epsilon }{}_{\mathbf{\Psi }}^{\mathbf{\prime }}}_{,n}<-\pi /R\hfill \\ 0\hfill & \mathrm{otherwise}.\hfill \end{array}\end{array}$$(36)

The fifth block is the integrator. In the time domain, it writes

$$\begin{array}{c}\hfill {{\mathit{u}}^{\mathbf{\prime }}}_{\mathrm{GD},n}={{\mathit{u}}^{\mathbf{\prime }}}_{\mathrm{GD},n-1}+{\mathsf{G}}_{\mathrm{GD}}{\mathit{\epsilon }{}_{\mathbf{\Psi }}^{\mathbf{″}}}_{,n}.\end{array}$$(37)

The same controller gain 𝖦_GD is applied to all baselines. It minimizes the closed-loop response time and maximizes robustness. The open-loop transfer function is the same for each baseline. It is mostly a pure delay caused by the moving average of 40 DITs, as stated by Eq. (15) in Sect. 3.5. The −3 dB cutoff frequency therefore depends on the sampling rate of the fringe tracker. It is 13, 4.5, and 1.5 Hz for sampling rates of 909, 303, and 97 Hz, respectively.

The sixth block is the matrix M^† to transpose the control signal from OPD space to piston space.

The last block is a quantization function. Practically, it means that the group-delay controller, when it detected a group-delay error larger than a fringe, only makes 2π phase-delay jumps until it comes to the reference fringe, without disturbing the long-term phase measurement. The command issued from the group-delay controller is therefore for each piezo-actuator the closest value that is a multiple of 2π. Hence

$$\begin{array}{c}\hfill {\mathit{u}}_{\mathrm{GD},n}={\mathrm{round}}_{\{2\pi \}}\left({\mathsf{M}}^{†}{{\mathit{u}}^{\mathbf{\prime }}}_{\mathrm{GD},n}\right).\end{array}$$(38)

When M^† u′_GD, n is between −π and π, the control signal of the group-delay controller is constant, and the phase delay can work without interferences caused by the group-delay loop.

We note that the final control signal (as shown in Fig. 2) also includes the modulation function, the fringe-search function, and the phase-delay control signal:

$$\begin{array}{c}\hfill \mathit{u}={\mathit{u}}_{\mathrm{GD},n}+{\mathit{u}}_{\mathrm{modulation}}+{\mathit{u}}_{\mathrm{search},n}+{\mathit{u}}_{\mathrm{PD},n}.\end{array}$$(39)

5.2. Modulation function and 2π phase jumps

The rounding of the group-delay command ensures that the same phase is always tracked even though jumps between fringes are required when the white-light fringe is searched for. However, this means that the science detector always records the fringes at the same phase delay, and when the sky is not well subtracted, for example, this can bias the visibility. This can be explained in the case of a perfect ABCD beam combiner. Assuming the recorded flux on each one of the four pixels is q_A, q_B, q_C, and q_D, the 𝖯2𝖵𝖬 matrix calculates the raw complex visibility in this way:

$$\begin{array}{c}\hfill V=q_A-q_C+i(q_B-q_D).\end{array}$$(40)

When the sky removal on q_A is not perfect, we have an additional flux that biases the measurement: q̂_A = q_A + ε_A. To remove this error term, the solution is to record the fringes with π offsets:

$$\begin{array}{cc}\hfill \widehat{V}& =q_A+{\epsilon }_A-q_C+i(q_B-q_D)\hfill \end{array}$$(41)

$$\begin{array}{cc}\hfill {\widehat{V}}_{\pi }& =q_C+{\epsilon }_A-q_A+i(q_D-q_B),\hfill \end{array}$$(42)

giving

$$\begin{array}{c}\hfill \frac{\widehat{V}-{\widehat{V}}_{\pi }}{2}=\frac{2q_A-2q_C+i(2q_B-2q_D)}{2}=V.\end{array}$$(43)

This temporal π modulation is added to the group-delay command (Fig. 5). It is not seen by the group-delay control loop because of the quenching block (block 4 in Fig. 4). This command is synchronized with the reset of the science detector, and it is sequentially:

$$\begin{array}{c}\hfill {\mathit{u}}_{\mathrm{modulation}}=\left(\begin{array}{c}0\\ 0\\ 0\\ 0\end{array}\right)\mathrm{or}\left(\begin{array}{c}\pi /2\\ \pi /2\\ -\pi /2\\ -\pi /2\end{array}\right)\mathrm{or}\left(\begin{array}{c}\pi /2\\ -\pi /2\\ \pi /2\\ -\pi /2\end{array}\right)\mathrm{or}\left(\begin{array}{c}\pi /2\\ -\pi /2\\ -\pi /2\\ \pi /2\end{array}\right).\end{array}$$(44)

Fig. 5. Upper panel: signal from the modulation block, in radians (the Volts-to-radians gain is not accounted for here). At each science exposure (here 5 s), the umodulation values change to make π offsets between the baselines. The signal repeats every four science exposures so that each baseline is observed with as many +π/2 as −π/2 offsets. Lower panel: modified sawtooth signal to search for fringes when the fringe tracker is in SEARCHING state. Here the research is made on all telescopes, meaning that the rank of matrix ${\mathbb{I}}_{\mathrm{GD}}^6$ is 0.

Only with a minimum of 4 DIT with each of the u_modulation above can we have, on each baseline, as many exposures with 0 and π offsets. Therefore the number of science DITs within a GRAVITY exposure is recommended to be a multiple of 4.

5.3. Fringe search

The search is made through a modified sawtooth function of increasing amplitude (u_saw). This function is started when the fringe tracker enters SEARCHING state. It is disabled when the system transition to TRACKING state. This is the only difference between the two modes. The same sawtooth function is generated for each piezo-actuator, but with different scaling factors: −2.75, −1.75, 1.25, and 3.25 for the first, second, third, and fourth beam, respectively. The four commands are then multiplied by the kernel of the ${\mathbb{I}}_{\mathrm{GD}}^4$ matrix. This matrix is computed as

$$\begin{array}{c}\hfill {\mathbb{I}}_{\mathrm{GD}}^4=\mathsf{M}{\mathbb{I}}_{\mathrm{GD}}^6{\mathsf{M}}^{†}\end{array}$$(45)

giving the following signal:

$$\begin{array}{c}\hfill {\mathit{u}}_{\mathrm{search},n}=({\mathsf{I}}^4-{\mathbb{I}}_{\mathrm{GD}}^4){\left(\begin{array}{cccc}-2.75& -1.75& 1.25& 3.25\end{array}\right)}^{\top }{u}_{\mathrm{saw}},\end{array}$$(46)

where 𝖨⁴ is the 4 × 4 identity matrix. This ensures that the baselines with sufficient fringe signal are tracking and are not modulated by the sawtooth function. When no more fringes are found, the kernel is the identity, and each baseline is modulated. This gives the trajectory shown in Fig. 5.

An example of fringe research and acquisition is presented in Fig. 6. In this example, at t = 0, the fringes are found and tracked on two baselines: the baseline between telescopes 1 and 3, and the baseline between telescopes 2 and 4. The group delay for the two fringes is zero, and the rank of the ${\mathbb{I}}_{\mathrm{GD}}^4$ matrix is 2. The fringe tracker is therefore in SEARCHING state. At t = 0.8 s, the fringes are found on baseline 34. The rank of the ${\mathbb{I}}_{\mathrm{GD}}^4$ matrix increases, and the system switches to TRACKING state. The group delays of baselines 12, 14, 23, and 34 are brought to zero. Signals are later found on all six baselines and the system reaches a nominal tracking state at t = 0 s.

Fig. 6. Phase Φi, j (upper panel) and group Ψi, j (lower panel) delays. All values are phases modulo 2π. The dashed lines correspond to the piezo command. From bottom to top, the phases correspond to baselines i, j equal to 12, 13, 14, 23, 24, and 34. At the beginning of this recording, the system only has fringes on two baselines (green and purple). At t = 0.8 s, the system found fringes on the yellow baseline, and later on all the other baselines. Gray areas correspond to the baselines whose S/N is below the value $S/N_{\mathrm{threshold}}^{\mathrm{GD}}$, meaning that no tracking is performed. After t = 1.4 s, the system nominally tracks on all baselines.

6. Phase-delay tracking

6.1. Principle

The GRAVITY phase control loop uses both piston-space and OPD-space state vectors. This is the easiest way to properly handle both piezo-actuators and the atmosphere dynamics (Correia et al. 2008). The vibrations and atmospheric perturbations are represented by six OPD-space state vectors labeled together x_V. Each vector corresponds to a baseline. The piezo-actuators are characterized by four piston-space state vectors x_P, one for each delay line. The phase delay, Φ_n, is a vector of measured phases:

$$\begin{array}{c}\hfill {\mathbf{\Phi }}_n={\left(\begin{array}{cccccc}{\mathrm{\Phi }}_{1,2}& {\mathrm{\Phi }}_{1,3}& {\mathrm{\Phi }}_{1,4}& {\mathrm{\Phi }}_{2,3}& {\mathrm{\Phi }}_{2,4}& {\mathrm{\Phi }}_{3,4}\end{array}\right)}_n,\end{array}$$(47)

which results from a linear combination of the state vectors x_V, n and x_P, n at a time n.

The real-time algorithm of the fringe tracker follows a sequence:

1. It predicts the future state of the system from previous state using the equation of state (Sect. 6.2),

2. It uses Kalman filtering to update the state vectors (Sect. 6.5), and

3. It uses the system state to command the piezo-actuators to correct for vibrations and atmospheric effects. (Sect. 6.6).

The phase-delay controller is summarized by the block diagram presented in Fig. 7.

Fig. 7. Block diagram of the phase-delay controller. The two control signals uGD and uPD are summed before they are applied to the actuators. In this control scheme, the state vectors ${\widehat{\mathit{X}}}_V$ are regulated. The control signal uPD is then issued from a matrix multiplication of the state vectors: ${\mathit{u}}_{\mathrm{PD}}=\mathsf{K}{\widehat{\mathit{X}}}_V$.

6.2. Equations of state

The equations of state are

$$\begin{array}{cc}\hfill {\mathit{x}}_{V,n+1}& ={\mathsf{A}}_{\mathsf{V}}·{\mathit{x}}_{V,n}+{\mathsf{B}}_{\mathsf{V}}·{\mathit{v}}_n\hfill \end{array}$$(48)

$$\begin{array}{cc}\hfill {\mathit{x}}_{P,n+1}& ={\mathsf{A}}_{\mathsf{P}}·{\mathit{x}}_{P,n}+{\mathsf{B}}_{\mathsf{P}}·{\mathit{u}}_n,\hfill \end{array}$$(49)

where each matrix is three-dimensional. Because we assumed uncorrelated baselines, only the diagonals are populated:

$$\begin{array}{cc}\hfill {\mathsf{A}}_V& =\mathrm{diag}\left(\begin{array}{cccccc}{\mathsf{A}}_V^{1,2}& {\mathsf{A}}_V^{1,3}& {\mathsf{A}}_V^{1,4}& {\mathsf{A}}_V^{2,3}& {\mathsf{A}}_V^{2,4}& {\mathsf{A}}_V^{3,4}\end{array}\right)\hfill \end{array}$$(50)

$$\begin{array}{cc}\hfill {\mathsf{B}}_V& =\mathrm{diag}\left(\begin{array}{cccccc}{\mathsf{B}}_V^{1,2}& {\mathsf{B}}_V^{1,3}& {\mathsf{B}}_V^{1,4}& {\mathsf{B}}_V^{2,3}& {\mathsf{B}}_V^{2,4}& {\mathsf{B}}_V^{3,4}\end{array}\right)\hfill \end{array}$$(51)

$$\begin{array}{cc}\hfill {\mathsf{A}}_P& =\mathrm{diag}\left(\begin{array}{cccc}{\mathsf{A}}_P^1& {\mathsf{A}}_P^2& {\mathsf{A}}_P^3& {\mathsf{A}}_P^4\end{array}\right)\hfill \end{array}$$(52)

$$\begin{array}{cc}\hfill {\mathsf{B}}_P& =\mathrm{diag}\left(\begin{array}{cccc}{\mathsf{B}}_P^1& {\mathsf{B}}_P^2& {\mathsf{B}}_P^3& {\mathsf{B}}_P^4\end{array}\right).\hfill \end{array}$$(53)

Here the upper index corresponds to the telescope or baseline numbers. In the block diagram of Fig. 2, the equations of state are written in the frequency domain, but the transfer function writes

$$\begin{array}{cc}\hfill \mathcal{Z}\{{\mathit{x}}_{V,n}\}& =\frac{{\mathsf{B}}_{\mathsf{V}}{z}^{-1}}{1-{\mathsf{A}}_{\mathsf{V}}{z}^{-1}}\mathcal{Z}\{{\mathit{v}}_n\}\hfill \end{array}$$(54)

$$\begin{array}{cc}\hfill \mathcal{Z}\{{\mathit{x}}_{P,n}\}& =\frac{{\mathsf{B}}_{\mathsf{P}}{z}^{-1}}{1-{\mathsf{A}}_{\mathsf{P}}{z}^{-1}}\mathcal{Z}\{{\mathit{u}}_n\}.\hfill \end{array}$$(55)

The evolution of the system is driven on one hand by white noise v_n, and on the other hand by a user-controlled voltage applied to the piezo-actuators v_n. When v_n is white noise, x_V, n has a colored noise, as highlighted by Eq. (54). The problem is similar for AO systems, where it has already been mentioned and corrected for (Poyneer & Véran 2010).

In Menu et al. (2012) and Choquet et al. (2014), we have shown that using several autoregressive (AR) models of order 2 in parallel was effective to correct for both vibration frequencies and atmospheric turbulence. However, practically, two issues made that implementation difficult: i) the determination of the vibration peaks for low S/N data, and ii) the need to change the space model when vibrations appear or disappear. Both problems can be technically resolved, but to ensure maximum robustness, we preferred a fixed well-defined space model. The idea is that the state x_V do change with time, but the space-state model does not. We therefore used an autoregressive model of order 30. This means that the state-space model corresponds to the 30 last values of the phase delay:

$$\begin{array}{c}\hfill {\mathsf{A}}_V^{i,j}=\left(\begin{array}{cccccc}{v}_1^{i,j}& v_2^{i,j}& v_3^{i,j}& \cdots & v_{29}^{i,j}& v_{30}^{i,j}\\ 1& 0& 0& \cdots & 0& 0\\ 0& 1& 0& \cdots & 0& 0\\ 0& 0& 1& \cdots & 0& 0\\ \cdots & \cdots & \cdots & \cdots & \cdots & \cdots \\ 0& 0& 0& \cdots & 1& 0\end{array}\right)\end{array}$$(56)

and

$$\begin{array}{c}\hfill {\mathsf{B}}_V^{i,j}={\left(\begin{array}{cccccc}1& 0& 0& \cdots & 0& 0\end{array}\right)}^{\top }.\end{array}$$(57)

The value 30 allows for complexity in the vibrational pattern while characterizing the perturbations with a sufficiently low degree of freedom. It was chosen in light of Kulcsár et al. (2012): they showed that it is a good compromise with respect to other state-space models to correct atmosphere and vibration for the tip-tilt of AO systems. The transfer matrix $A_V^{i,j}$ is determined every 5 s on a workstation that is connected on the reflective memory network.

The transmission matrix $A_P^i$ is also an autoregressive model, but of order 5. It includes the response function of the full system (image integration time, processing time, inertia of the piezo mirror, etc.), from setting a control voltage u_n to measuring a phase delay Φ_n. It is measured monthly by injecting a top-hat signal into the system. The matrix for the AR5 model verifies

$$\begin{array}{c}\hfill {\mathsf{A}}_P^i=\left(\begin{array}{ccccc}{a}_1^i& a_2^i& a_3^i& a_4^i& a_5^i\\ 1& 0& 0& 0& 0\\ 0& 1& 0& 0& 0\\ 0& 0& 1& 0& 0\\ 0& 0& 0& 1& 0\end{array}\right)\end{array}$$(58)

and

$$\begin{array}{c}\hfill {\mathsf{B}}_P^i={\left(\begin{array}{ccccc}1& 0& 0& 0& 0\end{array}\right)}^{\top },\end{array}$$(59)

where the $a_0^i$, $a_1^i$, $a_2^i$, $a_3^i$, and $a_4^i$ values, obtained experimentally on the fringe tracker are presented without normalization in Table 1. The same values can be represented in a Bode plot, see Fig. 8. The cutoff frequencies, calculated at a phase of −90°, are around 60 Hz. This is considerably below the bandwidth of the piezo-actuator measured by Pfuhl et al. (2010), who showed a cutoff frequency above 220 Hz. This is caused by the pure delays inside the system: detector integration time, processing time, and data transfer between the different units (LCUs). The static gain of the piezo $G_{\mathrm{piezo}}^i=\sum _{k=1}^5{a}_k^i$ is also an important property of the piezo-actuators because they are used both in the group- and phase-delay controllers. Practically, to remove the static gain of the control loop, and because u is in unit of OPD radians instead of Volts, the $a_k^i$ values are normalized by the static gain of the $G_{\text{piezo}}^i$.

Fig. 8. Bode plot of the frequency response of the four pistons in open loop (including pure delay caused by integration and computation). The −90° cutoff frequencies are 61 Hz, 56 Hz, 54 Hz, and 59 Hz. The cutoff frequency is dominated by the pure delay caused by the control loop system.

6.3. Observation equation

The observable, Φ_n, depends on the state vectors of the piezo-actuators x_P and on the state vectors of the atmosphere and vibrations x_V. Both contribute to the phase delay. The equation for the measurement process of the phase delay is a linear combination of the two:

$$\begin{array}{c}\hfill {\mathbf{\Phi }}_n=\mathsf{M}·{\mathsf{C}}_P{\mathit{x}}_{\mathrm{P},n}-{\mathsf{C}}_V{\mathit{x}}_{V,n}+{\mathit{w}}_n+{\mathbf{Ref}}_{\mathbf{\Phi }},\end{array}$$(60)

where w_n is a white measurement noise and 𝖢_V and 𝖢_P are three-dimensional measurement matrices with two-dimensional matrices on the diagonal for the vibrations and the piezo-actuator, respectively:

$$\begin{array}{c}\hfill {\mathsf{C}}_V^{i,j}=\left(\begin{array}{cccccc}1& 0& 0& \cdots & 0& 0\end{array}\right)\end{array}$$(61)

$$\begin{array}{c}\hfill {\mathsf{C}}_P^i=\left(\begin{array}{ccccc}1& 0& 0& 0& 0\end{array}\right).\end{array}$$(62)

In the case of a resolved astronomical target, the phase-delay vector Φ_n also includes a term that corresponds to the spatial signature of the target on the phase. It appears in the form of a non-zero closure phase (Θ^PD ≠ 0) and is included in the phase delay by the term Ref_Φ, as defined in Eq. (19).

6.4. Parameter identification

The a_1≤k≤5 predictive terms for each of the piezo-actuators are determined offline during dedicated calibration laboratory measurements. This is not the case for the identification of the v_1≤k≤30 AR values. They are calculated on a distinct Linux workstation. The workstation senses in real time the phase delay that passes by the reflective network. Every 5 s, it collects the last 40 s of data (36 000 sample at 909 Hz) and computes the pseudo-open-loop phase delay Φ_ATM, n by removing the influence of the piezo command and closure phases. This phase delay corresponds to the phase delay produced by the atmosphere and the vibrations only:

$$\begin{array}{c}\hfill {\mathbf{\Phi }}_{\mathrm{ATM},n}={\mathbf{\Phi }}_n-\mathsf{M}·{\mathsf{C}}_P{\mathit{x}}_{P,n}-{\mathbf{Ref}}_{\mathbf{\Phi }},\end{array}$$(63)

from which it calculates the differences between consecutive values:

$$\begin{array}{c}\hfill \mathbf{\Delta }{\mathbf{\Phi }}_{\mathrm{ATM},n}={\mathbf{\Phi }}_{\mathrm{ATM},n}-{\mathbf{\Phi }}_{\mathrm{ATM},n-1}.\end{array}$$(64)

The ΔΦ_ATM, n differential phase values are then processed baseline after baseline to derive 29 p_1≤i≤29 AR parameters that best represent the data. This uses the Python toolbox “Time Series Analysis”. The generated parameters ensure stationarity (Jones 1980) and thus provide stability to the closed-loop algorithm. Last, the v_1≤i≤30 values are determined through

$$\begin{array}{c}\hfill v_i=p_i-p_{i-1},\end{array}$$(65)

after which the matrices 𝖠_V and 𝖠_P are sent by the RMN to the real-time LCU to adjust the parameters of the phase-delay controller.

6.5. Asymptotic Kalman filter

The vectors ${\widehat{\mathit{x}}}_{V,n|n-1}$ correspond to our best estimate of the state of the vibrations at a moment n from all the measurements available up to a time n − 1. It can be estimated from ${\widehat{\mathit{x}}}_{V,n-1}$ according to the equation of state in Eq. (48):

$$\begin{array}{c}\hfill {\widehat{\mathit{x}}}_{V,n|n-1}={\mathsf{A}}_V·{\widehat{\mathit{x}}}_{V,n-1}.\end{array}$$(66)

However, the goal of the Kalman filter is to update this estimate from the error between the new observable and this new estimate. This error can be derived from the measurement process in Eq. (60) and writes

$$\begin{array}{c}\hfill {\mathit{\epsilon }{}_{\mathbf{\Phi }}}_{,n}={\mathbf{Ref}}_{\mathbf{\Phi }}-{\mathbf{\Phi }}_n+\mathsf{M}·{\mathsf{C}}_{\mathrm{P}}{\mathit{x}}_{\mathrm{P},n}-{\mathsf{C}}_V·{\widehat{\mathit{x}}}_{V,n|n-1},\end{array}$$(67)

which becomes, after wrapping around 0 (between −π and π) and baseline weighting,

$$\begin{array}{c}\hfill {\mathit{\epsilon }{}_{\mathbf{\Phi }}^{\mathbf{\prime }}}_{,n}={\mathbb{I}}_{\mathrm{PD}}^6.(\left[{\mathit{\epsilon }}_{\mathbf{\Phi },n}\right]\%2\pi ),\end{array}$$(68)

where %2π corresponds to modulo 2π. The ${\mathbb{I}}_{\mathrm{PD}}^6$ matrix is here for two purposes. The first is to derive the errors on low S/N baselines through the baselines with higher S/N. The second purpose is to ensure that the low S/N data are either weighted down or discarded through the weighted pseudo-inversion in Eq. (32). The state estimator ${\widehat{\mathit{x}}}_{V,n}$ is finally updated through an integrator:

$$\begin{array}{c}\hfill {\widehat{\mathit{x}}}_{V,n}={\mathsf{A}}_V{\widehat{\mathit{x}}}_{V,n-1}+{\mathsf{G}}_{\mathrm{PD}}·{\mathit{\epsilon }{}_{\mathbf{\Phi }}^{\mathbf{\prime }}}_{,n}.\end{array}$$(69)

The gain 𝖦_PD is not a scalar but a three-dimensional matrix. The correct estimate of 𝖠PD is the basis of Kalman filtering. In GRAVITY, it is not identified at each DIT, but every 5 s on the sensing Kalman workstation. Therefore, it is obtained from the asymptotic Kalman equations. It is calculated from the two covariance matrices of the measurement noise Σ_w and the steady-state error Σ_∞ (Poyneer & Véran 2010; Menu et al. 2012):

$$\begin{array}{c}\hfill {\mathsf{G}}_{\mathrm{PD}}={\mathrm{\Sigma }}_{\infty }{\mathsf{C}}_V^{†}{({\mathsf{C}}_V{\mathrm{\Sigma }}_{\infty }{\mathsf{C}}_V^{†}+{\mathrm{\Sigma }}_w)}^{-1}.\end{array}$$(70)

The steady-state covariance matrix can be obtained from the algebraic Riccati equation, the vibrations input noise, and the vibration and atmospheric noise Σ_v:

$$\begin{array}{c}\hfill {\mathrm{\Sigma }}_{\infty }={\mathsf{A}}_V{\mathrm{\Sigma }}_{\infty }{\mathsf{A}}_V^{†}-{\mathsf{A}}_V{\mathrm{\Sigma }}_{\infty }{\mathsf{C}}_V^{†}{({\mathsf{C}}_V{\mathrm{\Sigma }}_{\infty }{\mathsf{C}}_V^{†}+{\mathrm{\Sigma }}_w)}^{-1}{\mathsf{C}}_V{\mathrm{\Sigma }}_{\infty }{\mathsf{A}}_V^{†}+{\mathrm{\Sigma }}_v.\end{array}$$(71)

The noise characteristic (Σ_v, Σ_w) makes the Kalman filter optimally adapted to the average noise on the system, but not to instantaneous noise variations. Adaptability is the role of the ${\mathbb{I}}_{\mathrm{PD}}^6$ matrix in Eq. (68).

6.6. Determination of the control signal

To use the Kalman filter as a controller, an optimal command is needed. The purpose of this section is to determine the control signal u_PD, n from the atmospheric state vectors: ${\mathit{u}}_{\mathrm{PD},n}=\mathsf{K}{\widehat{\mathit{x}}}_V$. The task of this signal is twofold. First, it must cancel the atmospheric perturbations on the measured phase delay (Φ_ATM, n). Second, it must ensure that the phase delay converges to the group-delay control signal u_GD, n + u_modulation + u_search, n.

The first task uses the predictive power of the equation of state:

$$\begin{array}{c}\hfill {\mathit{u}}_{\text{PD},n}={\mathsf{M}}^{†}{\mathsf{C}}_V{\mathsf{A}}_V^{n_{\mathrm{DIT}}}{\widehat{\mathit{x}}}_V,\end{array}$$(72)

where ${\mathsf{A}}_V^{n_{\mathrm{DIT}}}$ is the power of matrix 𝖠_V over n_DIT samples. The integer n_DIT is there to account for the delay between control signal and its effect on the phase delay. Ideally, n_DIT = 1, but because of the pure delay in the open-loop transfer function, it must be higher. This pure delay depends on the frequency of the fringe tracker. We currently use n_DIT = 3 for the 909 Hz frequency and n_DIT = 2 for the lowest frequencies.

The second task is achieved through the integrator in Eq. (69). This integrator updates the state vector of the atmospheric perturbations to minimize the error between estimation and observation. It causes the quantity ε_Φ, n to converge to zero: ε_{Φ, n → ∞} = 0. From Eq. (67), we can derive

$$\begin{array}{c}\hfill {\mathbf{\Phi }}_{n\to \infty }={\mathbf{Ref}}_{\mathbf{\Phi }}+\mathsf{M}·{\mathsf{C}}_P{\mathit{x}}_{P,n\to \infty }-{\mathsf{C}}_V·{\widehat{\mathit{x}}}_{V,n\to \infty }.\end{array}$$(73)

Under a steady-state assumption, 𝖢_Px_P,n→∞ = u_n→∞ because the piezo-gains are normalized. From Eq. (72) we also have, under the same assumption, ${\mathsf{C}}_V·{\widehat{\mathit{x}}}_{V,n\to \infty }=\mathsf{M}·{\mathit{u}}_{\mathrm{PD},n\to \infty }$ because the first row of 𝖠_V is normalized. Hence, we have what we desire, i.e.,

$$\begin{array}{cc}\hfill {\mathbf{\Phi }}_{n\to \infty }& ={\mathbf{Ref}}_{\mathbf{\Phi }}+\mathsf{M}·({\mathit{u}}_{\infty }-{\mathit{u}}_{\mathrm{PD},n\to \infty })\hfill \end{array}$$(74)

$$\begin{array}{cc}& ={\mathbf{Ref}}_{\mathbf{\Phi }}+\mathsf{M}·({\mathit{u}}_{\mathrm{GD},n\to \infty }+{\mathit{u}}_{\mathrm{modulation}}+{\mathit{u}}_{\mathrm{search},n\to \infty }).\hfill \end{array}$$(75)

6.7. Why not a simpler state controller?

The block diagram presented in Fig. 2 shows the complexity of the fringe tracker. A simpler version was used during the first commissioning in 2016. The phase and group-delay controllers were both different. First, the group-delay control signal was not used to directly command the piezo-actuators, but instead, it was used as a setpoint for the phase-delay controller. Second, the phase-delay controller was a proportional-integral controller. It was efficient and robust and a good commissioning tool. However, at low S/N, its performance was not adequate. The reasons were that

• the phase-delay estimator is noisy. The proportional controller (or worse, a derivative) directly injected the noise back into the control signal. This was especially problematic during flux dropouts when the S/N can be close to zero; and that

• the phase delay is modulo 2π. It needs to be unwrapped with respect to a prediction. Using the setpoint as the prediction resulted in many fringe jumps during poor atmospheric conditions.

We therefore decided to have the group delay directly command the piezo, skipping the phase-delay controller. In parallel to the group-delay controller, we included a Kalman filter on the phase-delay feedback and used the predictive model to unwrap the phase (the vector ε′_Φ,n). This gave the block diagram in Fig. 7.

7. On-sky observations

7.1. Operation

The simplicity of operating the fringe tracker lies in the simplicity of the state machine. With only three states, the operator interactions are limited to the transition between IDLE, SEARCH, and back to IDLE. When the fringes are detected, the system automatically switches to TRACKING and the instrument then starts recording data. However, there are two free parameters that could ask for the intervention of the operator: the $S/N_{\mathrm{threshold}}^{\mathrm{GD}}$ and $S/N_{\mathrm{threshold}}^{\mathrm{PD}}$ thresholds. An operational error could be, for example, to set a $S/N_{\mathrm{threshold}}^{\mathrm{GD}}$ too low and risk having the fringe-tracker tracking on the second lobe of a fringe packet.

When the thresholds are correctly set, the system is made to be fully autonomous, and able to deal with any glitch of the VLTI. For example, Fig. 9 shows the fringe tracker losing telescope AT2 and how it recovered. The figures are the same plots as those available to the support astronomers during real-time GRAVITY operation. The data were captured and plotted at t = 0 s. At t = −6.39 s, the AT2 looses its pointing target, and the injected flux in the fiber drops. Immediately, the S/N level drops below the purple line, and the FT discontinues tracking on all the three baselines that include AT2. At this point, the system is still in TRACKING state. At t = −5.39 s, after a time delay of 1 s, the system switches to SEARCHING state, and the VLTI delay lines start following the sawtooth function. At t = −3.26 s, the system recovers the fringes on the AT2-AT1 baseline and starts centering them. At this point, the rank of the ${\mathbb{I}}_{\mathrm{GD}}^4$ matrix is back to 3 and the fringe tracker switches back to TRACKING state. At t = −3.12 s, the fringes are detected on all six baselines, and the system again tracks nominally.

In summary, the complexity of dealing with multiple baselines is hidden behind the ${\mathbb{I}}_{\mathrm{GD}}^6$ and ${\mathbb{I}}_{\mathrm{PD}}^6$ matrices presented in Sect. 4.4. From the user’s point of view, the fringe tracker transitions from SEARCHING to TRACKING state, but the engine behind the scene does not change the way it operates.

Fig. 9. Observation of star GJ 65 A on 21 November 2018. This example shows the case of a glitch on the AT2 adaptive optics system that resulted in the loss followed by the recovery of the fringe by the fringe tracker. Upper two panels: fringe tracker phase Φi, j and group delay Ψi, j estimators. The color lines correspond to each of the six baselines. The π/2 phase jumps at −7 and −1.5 s are normal and correspond to the modulation synchronized with the 5 s DIT science camera (the high spectral resolution detector). The gray areas correspond to detection of low S/N fringes by the controller. Lower left panels: S/N for each of the baselines calculated as the inverse of the square root of the phase variance ($1/\sqrt{\mathrm{Var}({\mathrm{\Phi }}_{i,j})}$). The horizontal red lines correspond to the group-delay threshold ($S/N_{\mathrm{threshold}}^{\mathrm{GD}}$). If the black line dips below the red line, the FT stops tracking on that baseline. The horizontal purple lines correspond to the phase-delay threshold ($S/N_{\mathrm{threshold}}^{\mathrm{PD}}$). Lower center panel: flux Fi. Similarly, the red line is a threshold made by a moving average of all flux, which is used to detect the loss of a telescope. Lower right panels: commands to the piezo-actuators and VLTI delay line actuators. The piezo-actuators take care of the fast control signal, while the VLTI delay lines are used to offload and to search for fringes over large distances.

7.2. Sensitivity

The sensitivity is mostly a question of having enough photons on the detector to generate a feedback signal for the fringe tracker. In Fig. 10 we plot the flux versus magnitude of calibrators observed with GRAVITY during a period covering June 2017 to November 2018. The selection of the files with a tracking ratio above 80% leads to a total of 1117 exposures on 473 distinct calibrators. The orange dots correspond to the 814 AT observations. The others correspond to the MACAO (visible AO) and CIAO (infrared AO) UTs observations. The solid lines correspond to a total transmission of 1% from telescope to the fringe-tracker detector. The dashed lines correspond to the theoretical detector noise, scaled by $\sqrt{N_p}$, with N_p the total number of pixels divided by the number of telescopes.

Fig. 10. Transmission plot, i.e., photons detected on the GRAVITY FT receiver, per telescope and per second, as a function of the K-band magnitude (tracking ratio > 80%). The magnitudes of the targets are obtained from SIMBAD. Most of the observations are made on-axis, meaning that 50% of the flux is lost because of the beam splitter. However, the many K = 9.7 mag CIAO observations are taken off-axis, hence the higher flux. Targets observed with ATs can be as faint as 10 mag. Based on this plot, it is clear that because the flux observed with the UTs is more than 10 times higher, observations will be possible up to K = 13 mag.

On the ATs, the faintest target observed so far was TYC 5058-927-1, a star with K = 9.4 mag. It was observed in the night of 5 April 2017, when the atmospheric conditions were excellent: seeing down to 0.4″ and coherence time up to 12 ms. The observations were made in single-field mode, where only half the flux is sent to the fringe tracker. The fringe tracker efficiently tracked the fringes throughout the entire exposure time, at a frequency rate of 97 Hz, with OPD residuals between 250 and 350 nm.

Technically, if we were to extrapolate to the UTs the sensitivity observed on the ATs, we should be able to track stars of K magnitudes up to 12.5. However, the faintest star observed with the fringe tracker on the UTs so far is TYC 7504-160-1, a star with K = 11.1 mag. It was observed in the night of 2 July 2017 during good atmospheric condition: seeing between 0.4″ and 0.6″, coherence time between 4 and 7 ms. The frequency rate was 303 Hz, with OPD residuals between 350 and 400 nm. The relatively low UT sensitivity is still not understood. One explanation could be that the 97 Hz integration time cannot be used: the fringe contrast is attenuated by the vibrations of the UT structure. Another explanation could come from moderate AO performances. With a low Strehl ratio, difficulties in injecting the light in the fibers decrease the number of available photons, but also create flux drops that decrease the visibility and complicate the fringe tracking. One last possibility is a selection effect caused by the lower availability of the UTs.

7.3. Signal-to-noise ratio

The limit for fringe tracking could in theory be an S/N of 1 per coherence time. Below that number, the phase varies faster than our ability of measuring it: our knowledge of the phase decrease with time, hindering the convergence of the fringe tracker. However, practically, we observed that an S/N of 1.5 per DIT is required to keep the fringes in the coherence envelope. The signal is proportional to the number of coherent photons received. The noise is created by quantum noise on one hand and the background noise on the other hand. At 909, 303, and 97 Hz, we measured during sky observations a standard deviation of σ_sky = 4, 5, and 8 ADU per pixel, respectively. These values come from the quadratic sum of the noises caused by the scattered metrology light (7 ADU at 97 Hz), the sky and environmental background (6 ADU at 97 Hz), and the read-out noise (4 ADU).

The detector variance observed during sky observations is used to compute the real-time S/N shown in the lower left panel of Fig. 9. In Fig. 11, for a dataset covering the seven months between April and November, we have plotted this S/N as a function of the measured flux. The three solid lines correspond to the theoretical S/N calculated from σ_sky and photon noise. These values lie below these theoretical lines because of a loss of visibility contrast. This can be caused either by non-equilibrium of the flux between the different telescopes, or by OPD variations within a DIT of the fringe tracker detector.

Fig. 11. Signal-to-noise ratio as a function of detected photon per telescope per DIT. The S/N is the computed by the real-time computer, defined by $1/\sqrt{\mathrm{Var}({\mathrm{\Phi }}_{i,j})}$ , as stated in Eq. (14). The three colors correspond to the three different frame rates of the fringe tracker. The solid lines correspond to theoretical values assuming 100% fringe contrast. The vertical lines are theoretical flux assuming observations on UTs with 1% throughput. The horizontal line corresponds to an S/N of 1.5.

The vertical lines correspond to the flux of a star of K magnitude of 8, 12, and 16 observed with the UTs under the assumption of 1% throughput. The drop at low flux of the theoretical S/N curves correspond to the effect of the observed sky noise. This noise is mostly detector noise at high frequency, and a combination of background and metrology noise at low frequency. The dotted lines are theoretical computations of the noise assuming only photon noise. Under this assumption, 1% throughput, and 100% visibilities, the UT sensitivity could technically reach magnitudes up to K = 16 mag.

7.4. OPD residuals and S/N

The fringe tracker performance degrades when the limiting sensitivity is approached. This is partly because the phase-delay control-loop gain decreases at low S/N because of the weighted inversion of matrix 𝖲^†_PD in Eq. (32). This is also partly caused by the difficulty of estimating the correct fringe-tracking state from a noisy measurement. This decrease in performance is shown in Fig. 12 for the same dataset as presented in Fig. 11. Concretely, below a S/N of 4, the residuals are above 300 nm. These residual include the variance caused by the measurement noise, which is equal to λ/(2π S/N): ≈90 nm for a S/N of 4. The fringe tracker residual does not correct the phase to that level, however: the residual remains a factor 2 or 3 above that value. At an even lower S/N (<3), the Kalman filter has problems to identify the state of the atmosphere and sometimes cannot update them fast enough. Fringe tracking is then only possible when the atmosphere varies slowly, when the atmospheric conditions are best.

Fig. 12. OPD residuals as a function of the fringe tracker S/N. Above an S/N of 10, the best fringe-tracking residuals are limited to a constant value around 80 nm, caused by the intrinsic bandpass of the fringe tracker hardware (dashed line). At lower S/N, the fringe tracker performance decreases because it is difficult to predict and track the evolution of the perturbations. The phase error caused by a readout noise on the phase goes as λ/(2π S/N) (lower dotted line). The actual performances of the fringe tracker shows that the degradation is more likely a factor 2 or 3 above that limit.

7.5. OPD residuals and τ₀

At high S/N (≤10), the accuracy of the fringe tracking depends on several environmental parameters: vibrations, wind speed, seeing conditions and coherence time (τ₀). Between wind seeing and coherence time, Lacour et al. (2018) showed that the strongest correlation is observed with τ₀. Using the same dataset as in Fig. 10 and using only the calibrators observed at 909 Hz, we plot in Fig. 13 the OPD residuals as a function of the coherence time as observed by the ESO site monitor at 500 nm. On average, the ATs perform better, with a median residual of 150 nm. Under the best conditions, the OPD residuals can be as low as 75 nm. The UTs residuals are higher, with a median value of 250 nm, and a minimum at 180 nm.

Fig. 13. OPD residuals as a function of coherence time (at 500 nm) for all calibrators and GC data taken between July 2017 and November 2018 in the 909 Hz mode. The dashed vertical lines correspond to the first quintile (3 ms) and last quintile (7 ms) of the coherence time distribution as observed by the ESO astronomical site monitor (ASM) at visible wavelength. The thick lines correspond to the first and last quintile of the OPD residual distribution (purple for UTs and red for ATs). The median fringe-tracking residuals are 150 nm on the ATs and 250 nm on the UTs.

Depending on the seeing conditions, the fringe tracker shows different limitations. In the worst atmospheric conditions, which make up 20%, (τ₀ <  3 ms), the UT and AT performances are limited by the coherence time of the atmosphere. Under these conditions, 20% of all observations have OPD residuals above 380 nm. The consequence is then that the jitter of the phase significantly affects the visibility in the science channel with long integration times. We can show that this contrast loss can be directly estimated from the variance of the phase according to the relation

$$\begin{array}{c}\hfill V_{\mathrm{residuals}}=⟨exp(i\mathrm{\Phi })⟩=exp(-{\sigma }_{\mathrm{\Phi }}^2/2).\end{array}$$(76)

With 380 nm rms jitters (σ_Φ ≈ 1.1 rad), the contrast of the fringes of the science detector will only be at V_residuals ≈ 55% of its maximum.

During the best 20% of atmospheric conditions (τ₀ >  7 ms), however, the fringe tracker can reach very low residuals. On the ATs, 80% of all observations are below 140 nm (V_residuals >  93%). However, on the UTs, the environmental vibrations cause the residuals to remain between 220 and 320 nm. The explanation is the higher level of vibrations observed on the UTs.

7.6. Power spectral density

In Figs. 14 and 15 we show the power spectral density of the OPD for two astronomical objects: IRS 16C observed with UTs, and GJ 65 observed with ATs. The observation conditions are listed in Table 2. The IRS 16C galactic center target was observed with the CIAO off-axis AOs. The conditions were good to excellent, and the residual OPD was around 220 nm over the entire duration of the observations. The GJ 65 binary was also observed in good seeing conditions, without AO, and OPD residuals of 140 nm.

Fig. 14. Power spectrum density and cumulative sum of the power spectrum density of phase residuals (2.2 μm/2π × Φn) toward target IRS 16C. The red curves are the pseudo-open-loop values, i.e., in the hypothesis of no fringe tracker (2.2 μm/2π × ΦATM, n). The blue curves correspond to power spectral densities of 1 f−2 μm2/Hz and 0.001 μm2/Hz. The six lower plots are the reverse cumulative sum of the power spectrum. At zero frequency, the sum reaches 220 nm on average over the six baselines.

Fig. 15. Same as in Fig. 14, except for the binary star GJ 65, which was observed with the ATs. The fringe-tracking residuals are 140 nm, lower than what is observed with the UTs. The main difference is the lower level of vibrations in the 10−100 Hz range.

To compute the power spectral density, the phase delay Φ_n is unwrapped and the π/4 modulation function removed. The pseudo-open-loop phase delay Φ_ATM, n is then computed from Eq. (63) using the piezo transfer function estimated during calibration. In Figs. 14 and 15, the six upper plots correspond to the square root of the power spectral density. The black curves correspond to the measured OPD estimated from the phase delay Φ, and the red curves show the power spectrum of the reconstructed atmosphere Φ_ATM, n. The atmospheric power spectrum (which also includes the vibrations) shows, as expected, the prevalence of the low frequencies in the atmospheric perturbations. Below 10 Hz, the power spectral density fits the relation

$$\begin{array}{c}\hfill {\mathrm{PSD}}_{{\mathbf{\Phi }}_{\mathrm{ATM},n}}(f<10\mathrm{Hz})=1f^{-2}\mu {\mathrm{m}}^2{\mathrm{Hz}}^{-1},\end{array}$$(77)

very well, which is not far, but different, from the power of −5/3 of a Kolmogorov type atmosphere. After fringe-tracking correction, the residual OPDs are attenuated below the spectral density:

$$\begin{array}{c}\hfill {\mathrm{PSD}}_{{\mathbf{\Phi }}_{\mathrm{n}}}(f<10\mathrm{Hz})\le 0.001\mu {\mathrm{m}}^2{\mathrm{Hz}}^{-1}.\end{array}$$(78)

The fringe tracker is therefore clearly a high-pass filter. The cutoff frequency of the fringe tracker correction is not well defined because it depends on the efficiency of the Kalman filter and the accuracy with which the equations of state reflect the system. On the UTs, the 3 dB cutoff frequency is on the order of 10 Hz. On the ATs, maybe because the vibrations are less frequent and the predictive model more accurate, the fringe tracker performance has a higher cutoff frequency at 30 Hz. This 30 Hz shows that the system is optimized because it is close to the cutoff frequency of the open-loop system (≈60 Hz in Fig. 8).

8. Discussion

8.1. Have we reached the ultimate sensitivity?

It is often assumed that the sensitivity of an optical interferometer must decrease as a function of N, the number of telescopes. This is not necessarily true. The limiting sensitivity is when each of the degrees of freedom of the fringe tracker reaches a variance of one during one coherence time of the atmosphere, and each of the degrees of freedom corresponds to a non-zero eigenvalue of matrix 𝕀 (Sect. 4.4). Therefore the threshold on the phase-delay control loop ($S/N_{\mathrm{threshold}}^{\mathrm{PD}}$) is applied on the eigenvalues in Eq. (30). For N ≫ 1, the signal grows like the flux (∝N), the degrees of freedom like N, and the number of pixels like the number of baselines, N². For an increasing number of telescopes, the detector noise therefore becomes ever stronger. However, in the case of a system without background and detector noise, the S/N on the eigenvalues no longer depends on N (Woillez et al. 2017). The ultimate sensitivity of the fringe tracker can then be established at one photon per coherence time and per telescope, regardless of the number of telescopes. For the VLTI, assuming 8 m telescopes, a coherence time of 10 ms, a throughput of 1%, and using the full K band, the ultimate sensitivity is Kmag = 17.5. Even assuming the need for an S/N of 1.5 per baseline and per coherence time, we should be able to reach a Kmag of 16 (Fig. 11).

What can be done? There are several paths forward to reach this magnitude. The first is to decrease the background and detector noise. The sky brightness in the K band at Paranal is about 12.8 mag/square arcsec³. Therefore, the fraction of sky background light entering the 60 mas single-mode fiber is almost negligible (Kmag_sky ≤ 19). A more important light emitter is the VLTI thermal background because of the ≈75% absorption of the optical surfaces. All of these surfaces are typically about T = 283 K to 288 K and contribute to the majority of the background light. Regarding the detector noise, even if the SAPHIRA detector has a readout noise below 1 e⁻, the number of pixels used by the fringe tracker per telescope is 36 (72 in split polarizations). For faint objects, this noise can therefore dominate. Solutions to use fewer pixels have been proposed (Petrov et al. 2014, 2016) and could lead to a better sensitivity. The development of infrared detectors could also be a promising path forward.

The second path is to increase the coherent throughput. On average, only 1% of the total flux reaches the fringe tracker detector. The 20 mirrors between M1 and the GRAVITY cryostat absorb up to 75% of the light. In addition, 50% of the light is scattered inside the IO beam combiner, 50% is lost by using the beamsplitter on the on-axis mode, and 50% because of the amplification of the detector (the so-called multiplicative noise). The last part (≈60% of the remaining flux) is lost when the light is injected into the single-mode fiber. Throughput improvements could therefore come from using fewer optical elements and possibly switching more of the mirrors involved from aluminium to gold coatings. Focus can also be placed on a more efficient light coupling into the single-mode fiber. This would mean a better AO system for a better Strehl ratio. The advantages of a good Strehl ratio are manyfold. It increases the mean throughput. It also maximizes the fringe contrast by giving a better instantaneous flux equilibrium. Last, it avoids flux drop-out, which is detrimental for good fringe tracking and model prediction.

Third, but not least, special care will be taken in monitoring and removal of the vibrations. The vibrations cause two main problems. Because they are usually at high frequency, they are difficult to predict and affect the fringe-tracking performance (in contrast to the atmospheric perturbations, which are easier to correct because they are at a lower frequency). The main problem, however, is that they limit how slow the fringe tracker can run because the decrease in fringe contrast hurts more than the additional integration time.

8.2. Have we reached the ultimate accuracy?

A false assumption is that sensitivity can be gained by trading accuracy for sensitivity because coherencing (i.e., keeping the fringes within the coherent length) requires fewer photons per coherence time than fringe tracking. However, by simultaneously using the group delay and the phase delay, it is possible to have the best of both worlds: the group-delay controller does the coherencing, while the phase-delay controller works in parallel as far as the S/N permits (Fig. 12).

Despite the high sensitivity, we routinely track high S/N fringes within 100 nm residual rms with GRAVITY. However, when the coherence time is short (τ₀ <  3 ms at 500 nm), the fringe tracker performance degrades (Fig. 13). This is caused by the open-loop latency of the fringe tracker (of about 4 ms). Observing during these conditions could clearly benefit from a fringe tracker with a shorter response time. Could we still improve the control loop during good atmospheric conditions, however?

The answer is yes. It lies in the proper management of the S/N by the Kalman filter. For convenience, we used an asymptotic estimation of the Kalman gain from the Riccati equation. A better Kalman filter would propagate errors as well as the state by also applying the equation of state to the covariance matrix:

$$\begin{array}{c}\hfill {\widehat{\mathrm{\Sigma }}}_{\mathit{x},n|n-1}={\mathsf{A}}_V·{\widehat{\mathrm{\Sigma }}}_{\mathit{x},n-1}·{\mathsf{A}}_V^{\top },\end{array}$$(79)

and derive the optimum gain each DIT. This could be achieved with additional computing power.

An additional amelioration could come from modal control. This was proposed in Menu et al. (2012) and could theoretically be implemented. However, we have currently not been able to find any practical implementation that would make it robust for a realistic environment.