© The Authors 2024

1 Introduction

Polarization is an essential part of the information contained in the electromagnetic radiation of astronomical sources. The use of polarimetric observations enables a better understanding of the source of radiation as well as its environment. Polarimetric observations are nowadays used in a very broad range of science cases, and an increasing number of instruments are equipped with a polarimetric mode (see, e.g., Witzel et al. 2011; Dorn et al. 2014; Norris et al. 2015; van Holstein et al. 2020). With the enormous success of the Very Large Telescope Interferometer (VLTI) beam combiner instrument GRAVITY in recent years in various science fields (see, e.g., GRAVITY Collaboration 2018a,b,c, 2019a,b), interest in polarimetric observations with GRAVITY and the VLTI has also grown. The fundamental capabilities of GRAVITY to make polarimetric observations have already been demonstrated via observations of polarization of flares from the supermassive black hole SgrA* (GRAVITY Collaboration 2018b, 2020b). With the help of the polarization data, it was possible to constrain the magnetic fields around SgrA*. Similarly, the Event Horizon Telescope (EHT) collaboration has studied the magnetic fields around the black hole M87* with their recently released polarimetric image (Event Horizon Telescope Collaboration 2021). However, polarimetry enables more than the study of magnetic fields; many other research areas profit from the availability of polarization measurements. For example, disks around young stellar objects can be studied with the help of polarimetry (Hunziker et al. 2021), and measurements of the dust properties of evolved stars benefit from polarization measurements (Ireland et al. 2005; Norris et al. 2012b; Haubois et al. 2019). A more complete overview of polarized observations is available in Elias et al. (2008) and Trippe (2014). To combine polarimetric measurements with the unique angular resolution of GRAVITY, we want to characterize the polarization properties of GRAVITY and the VLTI.

Ideally, a telescope and its instrument would not alter the polarization of incoming light. In reality, however, the optical train of a telescope influences the polarization signal. This can mean that the instrument produces a polarization signal, so-called instrumental polarization (IP), or alters the incoming polarization by introducing crosstalk, which mixes the incoming polarization states. To compensate for these effects, the telescope and its instrument must be carefully calibrated for their effect on the measured polarization signal. In this paper, we show the results of a series of measurements carried out to calibrate the polarimetric properties of the VLTI. This includes characterizing the amount of crosstalk between different polarization states and the IP introduced by the VLTI. In the case of an interferometer, this is more difficult than for a single-telescope instrument, as there are a significantly higher number of reflections and one also has to account for both of the rotations of the telescope, in elevation and in azimuth. For this reason, polarimetric observations with optical interferometers are not common yet, but the foundations were laid in the early 2000s (Elias 2001, 2004). Soon after, first steps were taken to study variable stars and circumstellar environments (Ireland et al. 2005; Rousselet-Perraut et al. 2006), followed by work in aperture masking (Norris et al. 2012a, 2015). Similar to GRAVITY at the VLTI, the MIRC-X at CHARA has conducted its first polarimetric observations (Setterholm et al. 2020).

While the modeling is more complicated for an interferometer than for a single telescope, there is no fundamental difference in the calibration model for absolute polarization. We were able to use similar calibration models to those used in solar physics (see, e.g., Beck et al. 2005; Harrington et al. 2019) and for the NACO and SPHERE instruments at the Very Large Telescope (VLT; Witzel et al. 2011; van Holstein et al. 2020). We then used our test data to adapt the model and constrain the polarimetric properties of the VLTI. With this approach, we have constructed a full calibration model to correct polarimetric observations.

Apart from the absolute effect the VLTI has on the polarization measurement, there is an additional effect that has to be considered for interferometers. If the light paths of the individual telescopes have different polarimetric properties, this can introduce differential birefringence between the different telescopes. The VLTI was built with great care to ensure that the different light paths and reflections within are as similar as possible, but of course they cannot be identical, as there are imperfections in the trains, as well as individual upgrades such as the adaptive secondary mirror at UT4. It is important to understand these differential effects since differential birefringence leads to a loss of fringe contrast and therefore limits the sensitivity of an interferometer (Beckers 1990; Rousselet-Perraut et al. 1996). Furthermore, differential effects can also introduce errors to the visibility phase, limiting the astrometric accuracy for polarized targets. This was already explored for a part of the VLTI by Lazareff et al. (2014a) and is continued with this work.

While most of the light path and the reflections are similar for the Unit Telescopes (UTs) and the Auxiliary Telescopes (ATs), we focus solely on UT observations in this work. For the ATs, the de-rotation of the field is done in the telescopes themselves, which adds more complexity to the polarization measurement. Furthermore, the ATs are not fixed in place and can be repositioned. This could affect the polarization, mainly if the telescopes are located on different sides of the delay line. Considering this and the scientific importance of the UTs, we decided to limit this study to the UTs.

The work presented here is split into two main parts. In the first part (Sects.2 –8), we develop a calibration model for the VLTI and GRAVITY to calibrate polarimetric observations with GRAVITY. For this part, we assume that all telescopes are identical, and we use the Stokes formalism. This is the formalism typically used for modeling instrumental effects, and it minimizes the necessary degrees of freedom in the model. This formalism is introduced in Sect. 2. We discuss the instrumental effects of the VLTI on the polarization with the telescope model and the calibration measurements in Sects. 3 and 4. In Sect. 5, we add GRAVITY to the model to complete it, and apply the model to on-sky data in Sects. 6 –8.

In the second part of the work (Sect. 9), we investigate the differential effects between the telescopes and how they affect observations. For this part, we analyzed the data of each telescope individually to measure the differential effects between them. We also switched to the Jones formalism, which requires more parameters to describe the polarization but can describe the propagation of the light phase through the telescopes; this is needed to understand the effects of differential behavior on the interferometric signal. We therefore introduce the Jones formalism and the necessary concepts in Sect. 9.1. Similar to the first part, we then fit the model again, but this time with the Jones formalism and for each telescope individually, to be able to constrain differential effects. The results of both parts are briefly summarized in Sect. 10.

2 Conventions – Stokes formalism

There are two different conventions for describing polarization (Collett 1992; Tinbergen 2005). One is the Stokes formalism (with Stokes vectors and Mueller matrices), and the other is the Jones formalism (with Jones vectors and Jones matrices). The Stokes formalism is often used to describe instrumental effects on polarization, as the components of the Stokes vector directly relate to the measurable intensities. It can also describe partial polarization and has simple formulas to measure and calculate the fundamental properties of polarized light. The Stokes values are also easily measured using a half- and a quarter-wave plate. One disadvantage is that the Stokes formalism does not include phase information, which we need to describe interferometric quantities. We therefore start with the Stokes formalism and later switch to the Jones formalism when the phase information is needed. This is the case when we look at differential effects between two telescopes in Sect. 9.

In the Stokes formalism, the light and its polarization are described by a Stokes vector. An electric field that is described by $$E(z,t)=\left(\begin{array}{c}{E}_x\\ E_y\end{array}\right)e^{i(kz-\omega t)}=\left(\begin{array}{c}{A}_x\cdot e^{i{\varphi }_x}\\ A_y\cdot e^{i{\varphi }_y}\end{array}\right)e^{i(kz-\omega t)},$$(1)

and the Stokes vector is defined as $$s=\left[\begin{array}{c}I\\ q\\ u\\ v\end{array}\right]=\left[\begin{array}{c}\langle E_x^2+E_y^2\rangle \\ \langle E_x^2-E_y^2\rangle \\ \langle 2E_x{E}_y\cos \delta \rangle \\ \langle 2E_x{E}_y\sin \delta \rangle \end{array}\right]=\left[\begin{array}{c}{A}_x^2+A_y^2\\ A_x^2-A_y^2\\ 2A_x{A}_y\cos \delta \\ 2A_x{A}_y\sin \delta \end{array}\right],$$(2)

with δ = ϕ_x − ϕ_y. For practical purposes, it is easier to define the Stokes vector with the measured flux at different angles: $$s=\left[\begin{array}{c}I\\ q\\ u\\ v\end{array}\right]=\left[\begin{array}{c}{F}_{00}+F_{90}\\ F_{00}-F_{90}\\ F_{45}-F_{135}\\ F_{\text{RH}}-F_{\text{LH}}\end{array}\right],$$(3)

where F₀₀ is the flux after a linear polarization filter at 0° and F_RH and _LH are the flux measurement for right- and left-handed circular polarization. As the absolute intensity is not important for the polarization properties, in the following we only consider normalized Stokes vectors: $$S=\frac{s}{I}=\left[\begin{array}{c}1\\ Q\\ U\\ V\end{array}\right]$$(4)

The first parameter is the intensity in the non-normalized Stokes vector and is one in the normalized Stokes vector. The second and third parameters, Q and U, represent linear polarization. Positive Q shows linear polarization in the vertical direction and negative Q in the horizontal direction. U is 45° rotated in the counterclockwise direction with respect to Q, looking toward the source. The V describes circular polarization, with a positive V being right-handed and a negative V being left-handed¹. The normalized Stokes vectors Q, U, and V range from −1 to 1.

A 4 × 4 real matrix, the Mueller matrix, describes the polarization change for any optical system. For an input state S_in the output state S_out is calculated as follows: $$S_{\text{out }}=M\cdot S_{\text{in }}.$$(5)

The general components of the Mueller matrix are described as follows: $$M=\left(\begin{array}{cccc}I\to I& Q\to I& U\to I& V\to I\\ I\to Q& Q\to Q& U\to Q& V\to Q\\ I\to U& Q\to U& U\to U& V\to U\\ I\to V& Q\to V& U\to V& V\to V\end{array}\right).$$(6)

Regarding instrumental effects, we looked at two main contributions in the Mueller matrix. The first one is induced polarization, often referred to as IP. This is described by the first column of the Mueller matrix: I → Q/U/V and corresponds to the polarization signal produced by the instrument. The second effect is the crosstalk between the states Q/U/V ↔ U/V/Q. The crosstalk introduces a mixing of the polarization states by the instrument. The first row of the Mueller matrix is often considered not important in astrophysical implications, as the elements Q/U/V → I do not play a significant role as, for most cases, the source polarization is small (Q, U, V ≪ 1). The quantities on the diagonal are the element I → I, which is one when we work with normalized Stokes vectors and the terms for polarimetric efficiency Q/U/V ↔Q/U/V, which describes how well polarimetric states are maintained.

One of the advantages of the Stokes parametrization is that it is very easy to calculate the essential polarization properties. From the Stokes parameters, one can calculate the degree of polarization (DOP), the degree of linear polarization (DOLP), and the polarization angle (Θ_pol) as follows: $$\text{DOP}=\sqrt{Q^2+U^2+V^2},$$(7) $$\text{DOLP}=\sqrt{Q^2+U^2},$$(8) $${\text{Θ}}_{\text{pol}}=\frac{1}{2}\arctan \left(\frac{U}{Q}\right)+n\cdot \frac{\pi }{2},$$(9)

where n is 1 for Q < 0 and otherwise 0.

Fig. 1 Schematic layout of the VLTI platform. The positions of the four UTs are shown as big circles, and the possible stations for ATs as small circles. The delay lines and the VLTI Lab are indicated in the center of the platform.

3 VLTI model

To build up a calibration model for the VLTI, we first modeled the light path with all its mirrors. This was done only once and not for all UTs individually, as the light paths of the four UTs are almost identical. The only differences are the distances between some mirrors and the directions of the first reflection in the delay line. This reflection changes depending on whether the telescope beam is on the left or right of the VLTI lab, so it is different between UT 1 & 2 and UT 3 & 4 (see Fig. 1). As this changes the incident plane by 180°, it does not affect the propagation of polarization. The overall model is based on what was previously developed by Lazareff et al. (2014a; for more details, see also Lazareff et al. 2014b).

3.1 VLTI light path

The light path for two telescopes is shown in Figs. 2 and 3. After the primary and secondary mirrors, the light is sent to the Nasmyth platform by M3. It then travels to the center below the telescopes, where it is guided into the Coudé room. In the Coudé room, it travels through the star separator, where some of the light is directed into the adaptive optics system. As the adaptive optics system is not sensitive to polarization and the light is not fed back, it is irrelevant for studying the polarization properties. From there, all light beams are sent to the delay lines. Here, the positions of the mirrors change slightly for each telescope, but as the distances only differ in the direction of light propagation and all the reflections are identical, these differences do not affect the polarization. From the delay lines, the light enters the VLTI lab, reaching the beam compressor, which adapts the beam size to fit the beam size required by GRAVITY. After the beam compressor, the light continues to the VLTI switchyard, where it can be sent to the individual instruments. For GRAVITY, there is one more reflection to feed the light into the instrument. In the instrument, the field is de-rotated by a K-mirror and then passes a half-wave plate (HWP) before it is fed into the fiber coupler. More details on the exact components of GRAVITY are given in Sect. 5.

Fig. 2 Schematic model of the VLTI light path. The light path is shown for two UTs, including GRAVITY, in red and blue. The rounded arrows indicate the possible rotations in the VLTI – for changes in azimuth (Az), elevation (El), and parallactic angle (Pa) – and in GRAVITY, at the field de-rotator and the half-wave plate (HWP). The straight arrows indicate the movement of the delay lines. The location of the star separator is only shown for the blue beam but is in the same place for both telescopes. The fiber optics and elements within the GRAVITY beam combiner are only indicated here and further explained in Sect. 5.

3.2 Modeling

An electromagnetic wave incident on a mirror can be decomposed into a component parallel (p-component) and a component perpendicular (s-component) to the plane of incidence. Reflections on a metallic mirror can introduce a linear polarization if the reflectivity of the two components is different or a circular polarization when there is a different phase shift for the two components. The Mueller matrix that describes such a reflection is given by (see, e.g., Collett 1992) $$M=\frac{1}{2}\left(\begin{array}{cccc}{r}_s^2+r_p^2& r_s^2-r_p^2& 0& 0\\ r_s^2-r_p^2& r_s^2+r_p^2& 0& 0\\ 0& 0& 2r_s{r}_p\cos (\delta )& 2r_s{r}_p\sin (\delta )\\ 0& 0& -2r_s{r}_p\sin (\delta )& 2r_s{r}_p\cos (\delta )\end{array}\right),$$(10)

where r is the reflection coefficient of each component and δ the relative retardation: δ = ϕ_s − ϕ_p. The r and δ can be directly calculated from the Fresnel formula: $$\sin {\Theta }_i=n\sin {\Theta }_t,$$(11)

where Θ_i and Θ_t are the angles of incident and transmitted light and n is the material- and wavelength-dependent refractive index. While this is the original Fresnel formula, the refractive index for metals is a complex number, and therefore, the reflection angle Θ_t is complex and is not a regular angle anymore. With the incident angle and the complex Θ_t, one can now calculate the reflectance: $$R_s=-\frac{\sin \left({\Theta }_i-{\Theta }_t\right)}{\sin \left({\Theta }_i+{\Theta }_t\right)}=r_s\exp \left(i{\varphi }_s\right),$$(12) $$R_p=\frac{\tan \left({\Theta }_i-{\Theta }_t\right)}{\tan \left({\Theta }_i+{\Theta }_t\right)}=r_p\exp \left(i{\varphi }_p\right).$$(13)

While the transmitted part is not relevant for metal surfaces, one can still use n and Θ_i to calculate the reflection coefficients and formulate the Mueller matrix for the reflection of a mirror.

As the form of the Mueller matrix given in Eq. (10) is rather nonintuitive, one can modify it by introducing the diattenuation, D: $$D=\frac{r_s^2-r_p^2}{r_s^2+r_p^2}.$$(14)

With this, the Mueller matrix of a metal surface takes the form (Chenault & Chipman 1993; Keller 2002; van Holstein et al. 2020) $$M=\frac{r_s^2+r_p^2}{2}\left(\begin{array}{cccc}1& D& 0& 0\\ D& 1& 0& 0\\ 0& 0& \sqrt{1-D^2}\cos (\delta )& \sqrt{1-D^2}\sin (\delta )\\ 0& 0& -\sqrt{1-D^2}\sin (\delta )& \sqrt{1-D^2}\cos (\delta )\end{array}\right).$$(15)

This form has the advantage that one can disentangle the different effects of a single mirror: the diattenuation shows the amount of IP. It has values between −1 and 1 and is 0 for no IP. The relative retardation, or retardance, δ introduces crosstalk for values below 180°. The prefactor to the Mueller matrix in Eq. (15) is only important for the total transmission and cancels out when working with normalized Stokes vectors. For the later calculations, we set it to one.

To get the incidence angles of the mirrors, we used the positions of the individual mirrors as given by Michel (2000). We show the notation for the mirrors, which we use in the following, in Fig. 3. From the positions, we could calculate the light path and the incident angle at each mirror. Together with the material of the mirrors, this is enough to set up the VLTI model. However, we implemented the following simplifications:

• We did not model M1 and M2 as the incidences are near normal, and we could neglect their contribution.

• The star separator (see Fig. 2) is not implemented as a special element in our model. As it is not rotating, we instead approximated it with the initial positions of M10 and M11 from the VLTI setup before the implementation of the star separator, as given by Michel (2000), see also our Fig. 3).

• The reflection in the delay line is done with a cat-eye retroreflector. In our model, we simplify this to three mirrors.

• The beam compressor is modeled as just one equivalent mirror, as all the incidences are very close to normal.

The materials for the individual mirrors and the incidence angles are listed in Table 1. The mirror M9 is a dichroic mirror, which reflects the infrared light and passes the optical light to the optical adaptive optics system. As this is not a simple metal surface, the polarimetric quantities for this mirror were measured by Lazareff et al. (2014a). The mirrors in the cat-eye are modeled as gold mirrors, but as mentioned before, this is just an approximation as the cat-eye retroreflector should not have a significant influence as all reflections are near normal. One thing to consider is that the silver mirrors in the train have a protective coating, which will lead to a different effect on polarized light.

The refractive indices of the mirrors are taken from the initial model from Lazareff et al. (2014a), as given in an online database². For the three used mirror materials, the refractive indices at a wavelength of 2250 nm are

• gold: n = 0.99 + 13.81i

• silver: n = 0.77 + 13.41i

• aluminum: n = 2.75 + 22.28i

GRAVITY operates in the K-band between 2000 and 2500 nm, which is why we used the refractive indices at the center of this band, at 2250 nm. The differences in refractive indices over the ±250 nm are on the order of ±(0.1 + 2i) and should not have a significant effect on the calibration model we are building up. With the entire model, we later analyzed the wavelength dependence more, which is discussed in Appendix B.

For the gold and aluminum, we could directly use these values, but for the silver mirrors, we had to consider that they have a protective coating. For these mirrors, we assumed a protective layer of Al₂O₃ with a thickness of 210 nm. These values fit the measurement we had available for one of the protected mirrors. The diattenuation and the phase shift can still describe the polari-metric properties of the protected mirrors, but the calculation is more complicated than outlined before. We used the method given by Jellison (1999) and Goldstein (2003). The method and a complete treatment of the protected mirrors are provided in Appendix A.

With the phase shift and diattenuation for each mirror, we could calculate the Mueller matrix for each mirror. To combine several Mueller matrices, one can take their product to get the combined Mueller matrix: $$M=M_n{M}_{n-1}\cdots M_2{M}_1.$$(16)

Our model starts at M3 and goes all the way down into the VLTI lab until the GRAVITY feeding optics.

Fig. 3 Simplified version of the VLTI light path from Fig. 2 to show the modeling and the experimental setup. The black rectangles show where the laser is launched and where the polarimeter is mounted. The names of the mirrors used in the text are given. The color of the mirror number shows the grouping that was used for the fitting. Gray mirrors are not fitted in our calibration model. Each colored group is located in one common plane: M4−M8 are in one vertical plane, and M10−M18 are in one horizontal plane.

3.3 Coordinate system and field rotation

To finalize the polarization model, one has to consider the 3D nature of the light path and that there are several fixed and varying field rotations in the path of the VLTI. The initial coordinate system is chosen so that Q aligns with north. From this start, the model is constructed by using the mirror positions from Michel (2000) and defining the light path as the vector from one mirror to the next. For each mirror, we calculated the angle of incidence and used this to define the direction of the s- and p-component before and after the mirror. For each pair of consecutive mirrors, we then compared the direction of the s- and p-components of the outgoing light with the same components of the incoming light of the next mirror. If there is a field rotation between the two mirrors, the comparison shows exactly this rotation, which we added to the model. In this way, we followed a plane through the VLTI that is orthogonal to the direction of light propagation. With this approach, we get the direction of s- and p- components and all field rotations. This model of following the propagating light while allowing for rotations in between does fully describe the light path, and no additional assumptions have to be made (for examples of similar approaches, see, e.g., Capitani et al. 1989; Beck et al. 2005; Balthasar et al. 2011; Harrington et al. 2019).

As the analysis is done for a telescope at the reference position (azimuth at 0° and elevation at 90°), we had to add the field rotation for azimuth and elevation by hand. The resulting field rotations are the following (see also Fig. 3):

• Between M3 and M4, there is a rotation due to the telescope movement in elevation (El). The rotation is z = 90° − El, as an elevation of 0^° corresponds to a zenith angle of 90^°.

• Between M8 and M9, there is a rotation depending on the telescope’s position in azimuth (Az). The rotation is given by ϕ = −(Az + 18.98°) + 6.02°. The 18.98° comes from the fact that the VLTI baselines are rotated by −18.98° compared to the east-west direction (see Fig. 1). As the zero position of the UTs is toward the south, this introduces an offset in the azimuth position. The 6.02° comes from a field rotation in the star separator mirrors.

• Between M9 and M10, there is a 90° rotation as the plane of reflection changes from reflections perpendicular to the ground (in the telescopes) to reflections in a horizontal plane (in the delay lines and the VLTI lab).

• One additional rotation, which is not coming out of the model, is the paralactic angle. This has to be taken into account as the Stokes parameters are defined toward the source and not in our chosen reference system (toward the north).

All these rotations are identical for the light paths of the four UTs. This leads to a total rotation of the field in the light path by $$\begin{array}{l}{\text{Θ}}_{\text{VLTI}}=\left({90}^{°}-\text{El}\right)-\left(\text{Az}+{18.98}^{°}\right)+{6.02}^{°}+{90}^{°}+\text{Pa}\\ =\text{Pa}-\text{El}-\text{Az}+{167.04}^{°}.\end{array}$$(17)

This is also stated in Gitton & Wilhelm (2003 , with the only difference being the sign of the azimuth angle. This is because the angle is defined as east of south in Gitton & Wilhelm (2003), while we used the convention of the ESO ISS system, which is east of north (Perraut & Berger 2010).

A rotated optical element would usually be implemented by multiplying the Mueller matrix of the element with a rotation matrix R. The Mueller matrix of that element is given by M_Θ = R(−Θ)⋅M⋅R(Θ), where R is (Collett 1992) $$R(\text{Θ})=\left(\begin{array}{cccc}1& 0& 0& 0\\ 0& \cos 2\Theta & \sin 2\Theta & 0\\ 0& -\sin 2\Theta & \cos 2\Theta & 0\\ 0& 0& 0& 1\end{array}\right).$$(18)

This was done to ensure that the input coordinate system was preserved. In our case, the field rotations are part of the optical system, and there is no advantage in preserving the input coordinate system. We therefore rotated the coordinate system at each field rotation. The response of an optical element and the field rotation is described by S = R(Θ)⋅(M⋅S′) = (R(Θ)⋅M)⋅S′ This way, the field rotation as given in Eq. (17) is automatically included in the final calibration Mueller matrix, with the resulting reference system then being defined as Q vertical in the lab. This approach gives the same result as using the usual convention that conserves the coordinate system and then rotating the reference system by the full angle given in Eq. (17).

One must take into account that each metallic mirror introduces a 180^° phase shift, which is equivalent to a change of coordinate system for a Stokes vector (Keller 2002). Due to this change in the coordinate system, rotations after an odd number of mirrors go into the total field rotation in the opposite direction. The full Mueller matrix of the VLTI, including all necessary rotations, is then given by $$\begin{array}{l}{M}_{\text{VLTI}}=M_{rmM18}\cdots M_{rmM10}\hfill \\ \cdot R\left({90}^{°}\right)\cdot M_{rmM9}\cdot R\left(-\left(\text{Az}+{18.98}^{°}\right)+{6.02}^{°}\right)\cdot M_{rmM8}\cdots M_{rmM4}\hfill \\ \cdot R\left({90}^{°}-\text{El}\right)\cdot M_{rmM3}\cdot R( \text{Pa}).\hfill \end{array}$$(19)

As mentioned before, M₁ and M₂ are omitted as they can be neglected for the polarization analysis.

3.4 Analyzing the input model

We use calibration data to fit the input model in the next section, but we can use the theoretical model to understand some of the principal properties of the VLTI mirror train and how it will affect polarization. In Fig. 4, the diattenuation and the phase shift of all the mirrors are shown. In an ideal case, the diattenuation would be zero (no IP), and the phase shift would be 180° (no crosstalk) for each mirror. One can immediately see that the mirrors with large reflection angles are the worst offenders. In terms of IP, the culprit is the dichroic mirror M9. As it is not a metal mirror and the values are measured individually, it does not lie on the lines of the three used metals. In terms of phase shift, the worst mirrors are the silver mirrors. This is due to the protective coating, which introduces a significantly stronger phase shift than an uncoated silver mirror would experience (see Appendix A). The dominant mirrors here are the mirrors M4, M12, and M16. M4 is at the Nasmyth platform of the telescope, and M12 and M16 are in the VLTI delay line. This introduces another effect, as there is the telescope rotation of the azimuth angle between them. This rotation worsens the effect of crosstalk, as we have mirrors with strong crosstalk terms with a field rotation in between them. This rotation adds a strong correlation with the azimuth position of the telescope to the polarization effects. The second rotation on the telescope is the elevation rotation, but the only mirror on a different side of the rotation than the others is M3. As shown in Fig. 4, M3 has, despite its 45° reflectance angle, comparably good polarization properties. We therefore expect fewer effects from the change in elevation than from the change in azimuth. A further effect we can infer from Fig. 4 is that while the values for the retardation are significantly different from the ideal value of 180^°, the values for the diattenuation are closer to zero. We therefore expect, in general, more crosstalk than IP from this model.

To verify these conclusions, we modeled the response of an arbitrary input Stokes vector. We used the vector S = (1, 0.3, 0.3, 0)^T. This corresponds to an input state with a DOP of 42.5% and an equal amount of linear polarization. The polarization angle is 22.5^°. This is the maximum source polarization one expects in astronomical sources in the near-infrared, but the exact vector is randomly chosen as an example. To see how this input vector propagates through the VLTI, we calculated the VLTI Mueller matrix for all telescope states and show the result for the output states in Fig. 5. In a second test, we repeated the calculation for unpolarized light. The polarization degree for this case is shown in Fig. 6. From these plots, we can verify a couple of conclusions:

• While there is a strong correlation with the Azimuth angle, there is only a small dependence on the telescope’s elevation.

• As the polarization quantities are only defined over a 180^° range, there is a repetition after 180^° of azimuth rotation.

• The DOP only varies by approximately ±2% from the input state, showing that IP is, in comparison to crosstalk, not the dominant effect for highly polarized targets.

• We see large variations in DOLP (±12%) and polarization angle (±13^°). This shows the large amounts of crosstalk in the VLTI path, shifting power between the Stokes parameters Q, U, and V.

• While the crosstalk dominates for highly polarized targets, we still expect some IP. This is, depending on the telescope position, in the range of 0–4%.

From this analysis, we conclude that our model produces output states that behave as expected and move on to calibrating it with test data.

Fig. 4 Diattenuation (top) and relative phase shift (bottom) of the individual VLTI mirrors in our model, shown as a function of the reflection angle for a wavelength of 2250 nm. The worst offenders are labeled in both plots. The colored lines show the continuous values for the three used materials: gold, silver, and aluminum. M9 does not lie on a line, as it is a dichroic mirror.

4 Calibration measurement

To verify and calibrate our model, we took calibration data at the VLTI. As a light source, we used a high-power thulium laser from IPG Photonics with a laser wavelength of 1908 nm. As the IP is expected to change with the wavelength, we used a laser at 1908 nm to be as close as possible to the science wavelength of GRAVITY (between 2000 and 2500 nm). The differences over 500 nm should be minor, as the refractive indices do not change significantly. In Appendix B, we analyze what differences our model predicts for this change in wavelength and confirm that the differences are minor compared to the uncertainties. Later, we could use calibration observations on sky to verify that it is not a limiting factor. The polarization measurements were done with a PAX polarimeter from Thorlabs. This is a rotatingwave plate-based polarimeter that Thorlabs customized to work at near-infrared wavelengths. To have the full light path of the VLTI, we launched the laser in the VLTI lab from the reference plates just in front of the GRAVITY feeding optics. With a linear polarizer, a half-wave, and a quarter-wave plate, we could modify the polarization of the laser and set it to arbitrary input states. The measurement head of the polarimeter was mounted onto a spider arm of one UT. This allowed us to measure the full light path at different telescope positions. (For a sketch of the experimental setup, see Fig. 3).

To get the complete polarization information, the goal was to measure the Mueller matrix of the light path at different telescope positions. As a Mueller matrix has 16 free parameters, we needed at least four input states to determine the entire matrix. As shown in Layden et al. (2012), Sabatke et al. (2000), and Reddy et al. (2014), it is best to equally space the input states over the possible parameter space, which minimizes the error from the matrix inversion.

To keep our input states as simple as possible, we used six different input states with equal distribution over the Poincaré sphere. The input states were four linearly polarized states with a distance of 45 degrees from each other and two fully circularly polarized states (left and right). The linear states were intentionally chosen not to coincide with the geometric axis of the light path (i.e., not 0^°, 45^°, …). The six input states we used are the following: 100% linear polarized at 75^°, 100% linear polarized at 30^°, 100% linear polarized at –15^°, 100% linear polarized at –60^°, 97% circular polarized, left-handed, and 97% circular polarized, right-handed.

The two circular states are only 97% circular polarized, as it was difficult to get a state with zero ellipticity in our test setup. For better reproducibility, we settled on those states. The measurement was done for all four UTs with an average of 12 telescope positions in altitude and azimuth. To test how well we can measure the polarization and how reproducible the input states are, we did separate test measurements. In these measurements, we placed the polarimeter in the VLTI lab directly behind the optics to modify the polarization of the laser. By repeatedly going through all different input states in the same manner as for the telescope measurements, we could estimate the reproducibility of these states. We conclude that the uncertainty on the polarization angle is on the order of 0.5° and 0.2% for the DOP. Those uncertainties are added to the error derived from the temporal scatter of each measurement.

Fig. 5 Properties of the output polarization for the input state S = (1, 0.3, 0.3, 0)T and all possible telescope positions. The plots show, from left to right, the DOP, the DOLP, and the polarization angle. The input state has a DOP and DOLP of 42.5% and a polarization angle of 22.5°. These values are shown as triangles in the color bars.

Fig. 6 Expected DOP for unpolarized light at all telescope positions.

4.1 First results

The first test we did was to check if the DOP is maintained or if we have significant depolarization in the light path of the VLTI. Overall, we measure a DOP of (98.1 ± 0.4)% and, therefore, a polarization loss of around 2%. Such a small amount of depolarization is expected from scattering on dust in the optical train and again indicates that there is not a large amount of IP.

More interesting is the effect of birefringence on the polarization angle and the DOLP. For the polarization angle, we looked at the four linear input states, 100% linearly polarized at 75°, –60°, –15°, and 30°. This measurement must be corrected for the field rotation, according to Eq. (17), with an additional correction for the fact that the polarimeter measurement head was mounted on the telescope spider, which is at an angle of 5.5^° from the central axis.

The measured polarization angles as a function of azimuth position are shown in Fig. 7. These data were taken at an elevation of 50^° and 70^° . As the change in elevation is not the dominant factor, the data are shown in the plot together. This is only true for these figures. Later, each telescope state will be fitted with its correct telescope elevation. While the measured values lie around the input values, there is some modulation of around 15^° . This is the crosstalk introduced by the mirror train, which clearly depends on the telescope’s position. We see the same effects and the same order of magnitude here, as we have seen in Sect. 3.4 and Fig. 5. With the data, we also show the prediction of the polarization angle by the polarization model from Sect. 3. The model is shown for an elevation of 70^°, as well as 50^° , illustrating again that the change in azimuth introduces a more dominant effect than the elevation change. We also see clearly that there is a 180^° ambiguity with the azimuth angle, which is, as discussed earlier, expected because the polarization properties are only defined in a range of 180^° .

The figure generally illustrates that the data roughly follows the model, and the amount of crosstalk we measured is nicely predicted by the model. The data from all telescopes are shown in one plot to illustrate that the telescopes behave very similarly. The differences between the UTs will be discussed in Sect. 9.

The third effect we could investigate is how much crosstalk there is between the linear and circular polarization states, that is to say, how elliptic the input states become. The result is shown in terms of the DOLP in Fig. 8. Here, we see that the linear states, which should have 100% linear polarization, have much lower values, going down to below 60%, again depending on the telescope position. The inverse effect is clearly shown for the circular states, which reach very high values in the DOLP with a maximum of 90%.

From the calibration data, one can conclude that the UTs behave very similarly and do not show significant depolarization. However, we clearly see polarization effects, which would modify a polarization angle measurement by up to 15^°. There is also substantial crosstalk between linear and circular states, which could decrease the measured linear polarization degree by up to 40%. Both effects are dependent on the telescope’s position. At this point, we have not done any fitting yet, but we can already say that the model predicts the data very well.

Fig. 7 Measured polarization angles for the different telescope positions, only corrected by the geometric rotation of the field. The four panels show the data for each linear input state, with the different colors showing the data from the four UTs. The data from 50 and 70° elevations are shown together. The gray horizontal line shows the input values for the polarization angles, the solid black line the value expected from the model for 70° elevation, and the dotted black line the model for 50° elevation.

Fig. 8 Measured DOLP for the different input states. Left and middle: linear input states, where the input DOLP is 100%. Right: circular states, where the input is fully circularly polarized, so 0% linear polarization. The data points are shown in different colors for the different states. The data from 50 and 70° elevations are shown together. The lines in the same color show the model prediction for each state (solid line for 70° elevation and dotted for 50° elevation). The dashed gray line shows the input values in all plots.

4.2 Fitting the calibration model

To reach a full calibration model, we improved the purely analytic VLTI model by fitting it to the obtained calibration data. The model includes 18 mirrors with two input values for the refractive indices and several rotations in the train. This has proven to be almost impossible to fit to our sparse data. To overcome this, we grouped all mirrors that have no rotation between them: M3, M4 to M8, M9, and VLTI lab and delay lines (M10 to M18). The groups of mirrors are also indicated in Fig. 3 by different colors in the mirror notation. The groups rotate then with a change of elevation between M3 and M4, a change in azimuth between M8 and M9, and a constant field rotation after M9: $$\begin{array}{l}{M}_{\text{VLTI}}=M_{\text{Lab}}\cdot R\left({90}^{°}\right)\cdot M_{rmM9}\\ \cdot R\left(-\left(\text{Az}+{18.98}^{°}\right)+{6.02}^{°}\right)\cdot M_{rmM4-8}.\\ \cdot R\left({90}^{°}-\text{El}\right)\cdot M_{rmM3}\cdot R( \text{Pa}).\end{array}$$(20)

The advantage of this approach is that the form of the Mueller matrix for a group of reflections stays the same as for a single reflection (see Eq. (10)). As the values in this matrix do not correspond to the values from a single Fresnel equation anymore, we could modify the Mueller matrix to the matrix in Eq. (15). This leaves us with two quantities for each mirror group, the diattenuation and the phase shift. The values for M9 were measured by Lazareff et al. (2014a), which leaves us with only six values to fit. Furthermore, the fitted matrices do not include rotations, which makes it possible to apply the model for each telescope position and in both propagation directions.

With the fitted values, we therefore have a polarization model that calculates a Mueller matrix for the whole VLTI light path and depends on the telescope position. For an ideal mirror, we assumed D = 0 and δ = 180° The values derived from the fit are listed in Table 2. The uncertainties are derived from the scatter of fitting each telescope individually and from bootstrapping the dataset. The results for the IP have a significant error bar relative to the values, but the values are all on the order of 10^-2 – 10^-3, which again shows that there is only a little IP. The values for the phase shift differ more strongly from δ = 180^, with the mirrors in the lab and the delay line (M10-M18) and the mirrors in the telescope (M4-M8) contributing equally. This again confirms the findings from the model in Sect. 3.4. The comparable values for 6 for both parts are expected, as we show in Fig. A.1, M4, M12, M16, and less strongly M6 should have the most significant impact on the retardance. The retardance of several mirrors adds up, and given that the worst offenders are situated in both parts of the train, one would not expect one part to be significantly better than the other. From our model, we would have expected that the lab and delay line part would be worse as it contains M12 and M16, but the fitting results show a comparable phase shift for the group M4-M8. This group contains exclusively protected silver coatings, and as mentioned earlier, they might have different protective coatings and are more challenging to model. We assumed this is the case for the strong retardance of this group.

With the fitted values, the calibration model is a simple function of telescope position. We obtain a Mueller matrix for each telescope position, which describes the IP of the VLTI, by executing Eq. (20). The sky polarization can be calculated from the measured Stokes vector and the Mueller matrix of the VLTI by applying Eq. (5).

With the fitted calibration model, we calibrated the test dataset, shown in Fig. 9. We obtained the Mueller matrix of each telescope position from our calibration model and multiplied it to the data. For the calibrated data, the polarization angle of the input states is recovered well, and the DOLP is 100% for the linear input states and very low for the circular states, which matches the input states. If one compares the calibrated data with the original in Figs. 7 and 8, this is a very clear improvement. The calibrated polarization states match the input states with an average error of 0.5 deg in the polarization angle and 0.4% in the DOLP. This excellent agreement of the fitted model also validates the simplifications made to the mirror train in the model (see Sect. 3.2).

Fig. 9 Test data calibrated with the fitted polarization model. The data are the same as in Figs. 7 and 8 but now calibrated. In both plots, the data are shown as dots, and the input states as dashed lines. The different colors show the different input states. The left plot shows the polarization angle, and the right plot shows the DOLP. Since all states are fully polarized, the recovered degree of circular polarization is (100 - DOLP).

5 Instrumental polarization of GRAVITY

So far, the results have been independent of the interferometric instrument and generally valid for VLTI observations in the near-infrared. However, to calibrate polarized observations, the instrument has to be taken into account as well. Here, we discuss the IP of the GRAVITY beam combiner (for a complete overview of GRAVITY, see GRAVITY Collaboration 2017).

5.1 GRAVITY light path

In GRAVITY, the light first passes the fiber coupler (Pfuhl et al. 2014). Part of the fiber coupler is a K-mirror to de-rotate the field and a HWP. The K-mirror, as well as the HWP, rotate in a fixed way during the observation: the K-mirror is used as a de-rotator and moves according to the field rotation described in Eq. (17). With the de-rotation of the field, it also de-rotates the sky polarization. Normally, we would therefore not need to derotate the polarization signal. However, the HWP in GRAVITY, which follows the K-mirror in the light path, rotates opposite the K-mirror and reintroduces the field rotation in the polarization. The reason for this is that it allows the metrology laser, which back-propagates through the light path, to have a stable polarization. This laser follows the full light path of GRAVITY and the VLTI before it is used to measure differential optical path differences above the primary mirror. This path difference is measured by using the interference between the light from the science and the fringe tracker beam. This measurement allows for phase references astrometry with GRAVITY but is extremely sensitive. To get the best possible contrast in the interference pattern, the polarization of the metrology beam is kept stable in the VLTI, and therefore, the polarization of the starlight in GRAVITY rotates with the field. The rotation correction described in Sect. 3.3 and Eq. (17) must still be applied to a polarization measurement. After the K-mirror and the HWP follow the tip-tilt, piston, and pupil control before the light is split into science and fringe-tracker and fed into optical fibers (Fig. 10). In the fibers, it passes the fiber control unit. Part of this unit are fibered polarization rotators, which rotate the polarization direction in the fiber and are used to match the polarization of all baselines (GRAVITY Collaboration 2017). From the optical fibers, the light is fed into the integrated optics system (Jocou et al. 2014; Perraut et al. 2018) and finally passed into the spectrometers (Straubmeier et al. 2014). In the spectrometers, there are Wollaston prisms, which can be put into the light path to allow for a polarimetric measurement (see Fig. 11) and split up the light into two polarizations (P1 and P2), with a 90° polarization angle between them. The polarization P1 is horizontally polarized in the VLTI lab frame or aligned with V (in the general Paranal coordinate system (V, W), where V is horizontal and W vertical in the lab; Gitton & Puech 2009). GRAVITY also includes a calibration unit (Blind et al. 2014), which can be used to test and calibrate the instrument. For this, it creates artificial stars in all beams. The calibration unit also includes the option to use a linear polarization filter to fully polarize the artificial light sources.

There are some field rotations in the light path of GRAVITY. However, the field rotations do not change over time. The GRAVITY fiber-coupler is aligned to ensure that a horizontal polarization on the calibration unit corresponds to one of the polarization directions on the detector. For this alignment, the linear polarizer in the calibration unit is used. With the linearly polarized light from the calibration unit, the fibered polarization rotators in the fiber control unit of GRAVITY are optimized to get a fully illuminated P2 spectrum on the detector, and no light on the P1, showing that the polarization vector is aligned with the vertical axis on the detector. The field rotations inside GRAVITY are, therefore, compensated by this alignment and do not have to be taken into account in this model. The fibered polarization rotators are made out of standard, not polarization-maintaining, fibers that can be rotated to adjust the polarization angle. The effect of this rotation is only to rotate the polarization. The fibers themselves are weakly birefringent, and no effect of the rotators on birefringence has been measured (see Perrin et al. 2024). The polarization effects of this system are, therefore, included in the measurement of GRAVITY and do not change with time.

Fig. 10 Schematic side and front view of the optical design of the GRAVITY fiber coupler. From GRAVITY Collaboration (2017).

Fig. 11 Optical design of the spectrometer. Most important for the polarimetric mode is the Wollaston prism, which can be moved in or out of the light path. From GRAVITY Collaboration (2017).

5.2 Measurements of polarization effects

To measure the polarization effects of GRAVITY, we performed two individual experiments. The first one was done with the same polarimeter as the VLTI measurements. We put the polarimeter in front of one of the beams and used the metrology laser as a light source. The second experiment uses the light of the calibration unit, which can be used with or without a linear polarization filter, producing fully polarized or unpolarized light. For the experiment with the calibration unit, we measured the signal on the detector, which can only measure one Stokes parameter at a time. To get the full linear polarization, we therefore needed to rotate the HWP between two exposures. This assumes that the HWP behaves as it should. This is where we used the polarimeter experiment to confirm that this is the case. A sketch of the two experiments is shown in Fig. 12, and we discuss them in detail in the following.

5.2.1 Polarimeter measurement

In the first measurement, we used the metrology laser as a light source. The metrology laser in GRAVITY is split into three parts: two low-power parts (with less than 1% of the laser power), which follow the full light path, and one high-power part fed into the light path only after the fiber optics. We used this last part, the so-called carrier beam, and measured it with the PAX polarimeter outside GRAVITY. The carrier beam follows some mirror optics and then passes the HWP and the de-rotator before leaving GRAVITY (as indicated in the left part of Fig. 12). We used the polarimeter to measure the outgoing light for each beam and rotated the HWP. The measured Stokes Q and U values are shown in Fig. 13.

The carrier beam we used in this measurement is fully polarized, with a fixed polarization state. However, the polarization direction was only coarsely aligned during the integration of GRAVITY, so we only have a vague idea about the polarization angle. Furthermore, the carrier’s light is fed into the light path at the fiber coupler, and we therefore do not have a measurement of some parts of the instrument, mainly the fibered and integrated optics. we could not use this measurement to characterize GRAVITY fully, but we could use it to characterize the HWP. The measured Stokes parameters are shown in Fig. 13 as a function of HWP rotation angle. As solid lines, the response of a perfect HWP is shown. The measured Stokes parameters show a very good agreement with the theoretical expectation, with an average discrepancy of 0.005. Given that small value in comparison to the other uncertainties in the polarization calibration, we decided to treat the HWP as an ideal component and will describe it with the following Mueller matrix: $$M=\left(\begin{array}{cccc}1& 0& 0& 0\\ 0& 1& 0& 0\\ 0& 0& -1& 0\\ 0& 0& 0& -1\end{array}\right)$$(21)

For all other elements of GRAVITY, we refer to the results from the second measurement.

Fig. 12 Sketch of the two polarimetric measurements for GRAVITY. Left: in the first measurement, the light from the metrology laser is used and measured with a polarimeter. Right: the light from the calibration unit is used and recorded on the spectrometer. In both cases, the light direction is indicated by arrows. The light path is the same as shown in Fig. 10, with the reflection for the Tip-Tilt and Piston control just shown as a change in direction in the light path.

Fig. 13 Measurement of the polarization in GRAVITY with a rotating HWP. The points show the data points for Stokes Q (red) and Stokes U (gray). The black lines show the corresponding values for an ideal optical element.

5.2.2 Calibration unit measurement

In the second experiment, we used the light from the calibration unit of GRAVITY, which is a Quartz Tungsten Halogen lamp with the possibility to add a linear polarization filter. Using this lamp as a light source, with the polarization filter, we get a constant and linearly polarized input source. We then rotated the K-mirrors of each beam from its initial position to 360°. At each location of the K-mirror, we took detector frames with the HWP at 0 and 22.5 deg. The mechanical rotation of 22.5 deg corresponds to a rotation of the polarization angle by 45 deg. At each position, we extracted the measured flux from the detector, split up into P1 and P2 by the Wollaston prism. Following Eq. (3) allows the linear stokes parameter to be measured. With this experiment, we measured the polarization on the GRAVITY science detector for a linear input polarization with a rotating polarization angle, shown in the top panels of Fig. 14. As there is no quarter-wave plate in GRAVITY, there is no possibility to measure Stokes V, but only the linear polarization parameters Q and U . As this leaves us with an incomplete measurement for the Mueller matrix of GRAVITY, we repeated the full measurement, but without the linear polarization filter of the calibration unit. The input light is then unpolarized and allows us to measure the IP of GRAVITY. The data are shown in the bottom panels of Fig. 14.

GRAVITY has two different observing modes, which come with slightly different optical paths. The first, the off-axis mode, is the mode in which two different objects are observed as science (SC) and fringe-tracker (FT) objects. In this mode, the field is spatially split into two fields by a roof prism in GRAVITY. The second mode is the on-axis mode, in which the same object is used for SC and FT. In this case, the light is split by a beam splitter in GRAVITY. Apart from the roof and the beam splitter, the light path is identical for both modes in GRAVITY. However, as these two elements are distinctly different, one has to assume that this affects polarization. Therefore, we did the full calibration unit test twice, once for the on- and for the off-axis mode. This leaves us with four measurements in total: the polarized and unpolarized input in both on- and off-axis mode. All four measurements are shown in Fig. 14. We repeated these measurements for all four GRAVITY beams, and as the results were very similar, they are shown as averaged data points in Fig. 14.

First look at calibration unit measurement. Looking at the data with polarized input light (top row in Fig. 14), one sees that GRAVITY in first order behaves as expected: the light is almost entirely linearly polarized, and the polarization direction rotates with the rotation of the K-mirror. For the unpolarized input (bottom row in Fig. 14), the measured polarization is much lower but not zero. Here, we see the main difference between the off-axis and the on-axis mode. In off-axis, the measured linear polarization for unpolarized input light is between 0 and 8%. In on-axis mode, we measure 8 to 16% for the same input light. We therefore have a much higher IP in the on-axis mode, which can be explained by the fact that beam splitters often show nonideal polarization properties.

Fit to the data. To fit a model to the test data, we separately fit a Mueller matrix for the K-mirror and GRAVITY. The HWP is assumed to be perfect, as before: $$\begin{array}{l}{M}_{\text{GRAVITY}}=M_{\text{GR}}\cdot R\left({\Theta }_{\text{HWP}}\right)\cdot M_{\text{HWP}}\cdot R\left(-{\Theta }_{\text{HWP}}\right)\\ R\left({\Theta }_{\text{KM}}\right)\cdot M_{\text{KM}}\cdot R\left(-{\Theta }_{\text{KM}}\right),\end{array}$$(22)

where M_GR is GRAVITY without the HWP and the K-mirror. As discussed before, the polarization properties of the on-axis and off-axis modes of GRAVITY are different. One therefore has to choose M_GR,on or M_GR,off here, depending on the observing mode. For simplicity, we use M_GR as a general name for it. Furthermore, R (Θ_HWP) and R (Θ_KM) are the rotation matrices for the HWP and the K-mirror. The K-mirror consists of three individual mirrors, which are fixed in one mount that moves them together. We can describe them as one Fresnel reflection according to Eq. (15). This leaves only two fitting parameters (diattenuation and retardation). The fit for the remaining optics of GRAVITY is more complicated. There are a large number of field rotations inside GRAVITY. They are all fixed and do not change with time, but this prohibits us from using a simple formula as for the K-mirror. We could split GRAVITY into individual mirror groups, as we did for the VLTI, with a rotation relative to each other. However, this would not lead to a decrease in fitting parameters. Additionally, we did not need the freedom to fit the parameters individually, as the rotations do not change. So for GRAVITY, we just fit a full Mueller matrix (according to Eq. (6)) with only the first component being fixed to 1 as we looked at normalized Stokes vectors. All other values can take values between −1 and 1.

In summary, we are using the data from Fig. 14 and fit Eq. (22) to it by fitting each component of M_GR , as well as the diattenuation and retardation M_KM , while not fitting M_HWP , but using Eq. (21) for it. As an initial guess, we took a perfect mirror for the K-mirror, and we calculated the expected Mueller matrix for the rest according to the mirror positions and materials in GRAVITY. The fit was then done twice, once for the off-axis mode and once for the on-axis mode.

The result of the fit is shown in Fig. 14, and the resulting Mueller matrices can be found in Appendix C. The comparison to the data in Fig. 14 shows that the fit worked well, and the instrument is well described by the given matrices.

Fig. 14 Measured polarization with GRAVITY in the different observing modes. The left column shows the measurement in the off-axis mode, and the right column in the on-axis mode. In the top row, the linear polarization filter is used for the input light source; in the bottom, it is not. In all plots, the Stokes Q data points are shown in red and orange and the Stokes U in gray. The data are the average over all four GRAVITY beams. The results from the fitted model are shown in black lines.

6 Full calibration of VLTI and GRAVITY

With the results of the GRAVITY measurement, we have all the information in hand to calculate the complete polarimetric response for an observation with GRAVITY and the VLTI. From the work presented in Sect. 4 and especially Eq. (20) we obtain the Mueller matrix of the VLTI, which depends on the elevation and azimuth of the telescope and is entirely defined by the parameters in Table 2. Together with the Mueller matrix of the K- mirror and the rest of GRAVITY, which are listed in Appendix C, the complete polarization response of GRAVITY and the VLTI is given by the following equation: $$\begin{array}{l}{M}_{\text{ALL}}=M_{\text{GR}}\cdot R\left({\Theta }_{\text{HWP}}\right)\cdot M_{\text{HWP}}\cdot R\left(-{\Theta }_{\text{HWP}}\right)\cdot \\ R\left({\Theta }_{\text{KM}}\right)\cdot M_{\text{KM}}\cdot R\left(-{\Theta }_{\text{KM}}\right)\cdot M_{\text{VLTI}}(Az,El,pa),\end{array}$$(23)

where R is the usual rotation matrix (Eq. (18)) with the position of the HWP (Θ_HWP) and K-mirror (Θ_KM). For M_GR, one has to choose the on- or off-axis one, matching the observing mode.

7 Polarimetric measurements with GRAVITY

To observe with GRAVITY in a polarimetric mode, one has to put the Wollaston prism into the light path, which can be selected in the preparation of the observing block. This gives measurements in two polarizations, P1 and P2, as discussed in Sect. 5. From the interferometric signal in both polarizations, one can calculate the measured source intensity and the first linear Stokes parameter, q: $$I=I_{\text{P}1}+I_{\text{P}2},q=I_{\text{P}1}-I_{\text{P}2}.$$(24)

To retrieve the second linear stokes parameter, u, one has to rotate the HWP by 22.5°, which rotates the polarization axis by 45°: $$u=I_{\text{P}1}^{22.5}-I_{\text{P}2}^{22.5}.$$(25)

The rotation of the HWP can also be selected in the preparation of the observation and is added as an offset to the nominally calculated HWP position. This means that the HWP still follows its calculated position (see Sect. 5.1) but in the second case with an added offset of 22.5 degrees.

With this q and u, one usually builds up a measured stokes vector and corrects it with the instrument matrix. However, this only works if the HWP is the last optical element in the beam. In the case of GRAVITY, the HWP sits at the entrance of GRAVITY, which means several optical elements after the HWP, which introduce birefringence. Therefore, what the detector sees is not the U state of the light but a Q with an altered instrument setup. We therefore assumed that we measured Q twice, once with the HWP at 0^° and one at 22.5^° . In normalized Stokes parameters, the measurement is then given by $$Q_1=\frac{I_{\text{P}1}^0-I_{\text{P}2}^0}{I_{\text{P}1}^0+I_{\text{P}2}^0},Q_2=\frac{I_{\text{P}1}^{22.5}-I_{\text{P}2}^{22.5}}{I_{\text{P}1}^{22.5}+I_{\text{P}2}^{22.5}}.$$(26)

Usually, Q₁ and Q₂ are taken in two subsequent exposures. One then needs a separate Mueller matrix for both measurements. This matrix is given by Eq. (23), which then includes the VLTI and GRAVITY, with all field rotations and polarimetric effects. It depends on the Azimuth and Elevation of the telescope, on the parallactic angle, and the position of the HWP and K-mirror. All this information is taken from the header of a normal GRAVITY fits file, with the keywords listed in Table 3. MALL has to be calculated twice for the two exposures used for Q₁ and Q₂ in Eq. (26). Using the header values for the two subsequent exposures gives two Mueller matrices, M_ALL,1 and M_ALL,2 (The HWP offset in case two is already taken into account in the HWP values in Table 3).

With M describing the response of the system in S_i = M_i ⋅ S_sky and the components of M from Eq. (10), the measured Q_i can be written as $$\begin{array}{l}{Q}_i={(}_I\cdot I_{\text{sky }}+{(}_Q\cdot Q_{\text{sky }}\\ +{(}_U\cdot U_{\text{sky }}+{(}_I\cdot V_{\text{sky }}.\end{array}$$(27)

Furthermore, we can define the polarization vector on sky as S Sky = (1, Q, U, 0). The I was set to 1 as we were assuming normalized Stokes parameters. We could set V to zero as sources in the near-infrared usually do not show significant circular polarization. This changes the given expression to $$Q_i={(I\to Q)}_i+{(Q\to Q)}_i\cdot Q_{\text{sky }}+{(U\to Q)}_i\cdot U_{\text{sky }}.$$(28)

Separating the IP and the crosstalk then allows a system of equations to be constructed, $$\begin{array}{l}\left[\begin{array}{c}{Q}_1-{(I\to Q)}_1\\ Q_2-{(I\to Q)}_2\end{array}\right]=\left[\begin{array}{cc}{(Q\to Q)}_1& {(U\to Q)}_1\\ {(Q\to Q)}_2& {(U\to Q)}_2\end{array}\right]\cdot \left[\begin{array}{c}{Q}_{\text{sky }}\\ U_{\text{sky }}\end{array}\right],\\ Q_{\text{vec }}=A\cdot {\left[Q_{\text{sky }},U_{\text{sky }}\right]}^T\end{array}$$(29)

which can be solved using linear least squares: $${\left[Q_{\text{sky }},U_{\text{sky }}\right]}^T={\left(A^{T}A\right)}^{-1}{A}^T{Q}_{\text{vec }}.$$(30)

This gives a measurement of Q and U on sky for every two exposures. In case the source polarization is assumed to be constant over time, Eq. (29) can easily be extended with more than two Q measurements (see van Holstein et al. 2020). Another solution to this approach is to forward model the polarimetric property of the source (GRAVITY Collaboration 2020b) to recover the complete polarization information. In this case, M_ALL can just be taken from Eq. (23).

Fig. 15 Polarization measurements of IRS 16C between 2019 and 2022. The top panel shows the polarization angle for each night, and the bottom panel shows the polarization degree. The average over all nights is shown as a horizontal line

8 Application to data

To test the calibration model, we used data taken with GRAVITY of the Galactic Center. A description of the observation of the Galactic Center black hole and its surrounding stars can be found in, example, GRAVITY Collaboration (2018a,b, 2020a). As shown in those papers, the science target SgrA* is not a good test target as it has a variable polarization state. However, the observations in the Galactic Center are done in the off-axis mode. In this mode, a close-by single star is used as a phase-reference source (for more information, see GRAVITY Collaboration 2017). In the case of the Galactic Center, the phase reference source is the star IRS 16C, which has a brightness of m_K = 9.55 and is well known to be slightly polarized. It was, for example, observed in Ott et al. (1999), Witzel et al. (2011), and Buchholz et al. (2013), and the studies found a consistent polarization of 4% (see Table 4) due to foreground polarization from dust.

In the off-axis mode, GRAVITY acts as two individual interferometers, which are referenced to each other via a metrology system. While the goal of the measurement was to measure the polarization of Sgr A*, we equally get an interferometric observation of the phase reference target. As the rotation of the HWP affects this target in exactly the same way as it does the science target, we applied the full analysis to the measured fluxes of the fringe-tracking star to measure its polarization. As this is a single star, one can simply take the coherent flux calculated by the GRAVITY data reduction software (DRS). We used the flux per polarization channel, took the average over the spectrum and all telescopes, and calculated Q_i for each exposure following Eq. (26). We then get Q_sky and U_sky by using Eqs. (29) and (30).

As IRS 16C is the usual fringe-tracking target for Galactic Center observations, there are several datasets with polarization observations of the star. Figure 15 shows the polarization angle and degree for different nights between 2019 and 2022. The polarization state in each night is well measured with an average scatter of the polarization angle of 5° and 0.4% for the polarization degree per night. This shows that the polarization calibration works well for different telescope orientations and that we can measure reliable polarization states with GRAVITY. With the complete analysis, we measure a polarization degree of (4.3 ± 0.4)% at an angle of (19.9 ± 4.6)°. These values perfectly agree with the other values from Table 4.

The fact that we can measure the expected polarization of the fringe-tracking star verifies our model and the polarimetric capabilities of GRAVITY. It also confirms the results from GRAVITY Collaboration (2020b), where we already used the presented calibration model and could show clear and varying signals in the polarization of SgrA*.

The way to extract the flux in the different polarization states from GRAVITY data depends on the observed target. While for single targets, such as IRS 16C, the DRS outputs can be used, a more complicated source structure might involve a different analysis. For the case of SgrA*, for example, we applied a multisource fit, which fits a model of the central black hole and several surrounding stars to the data, including the flux ratio between the stars and the black hole. By doing this for both polarizations individually, we could measure the flux ratio with respect to nearby stars in each polarization (for more details, see GRAVITY Collaboration 2020a). The flux measurement in the different polarization states will, therefore, depend on the science target. After that, the full calibration is comparably easy. The values for telescope azimuth and elevation, as well as the position of the GRAVITY K-mirror and HWP, are given in the header of each file (see Table 3). With this information, the Mueller matrix, as given in Eq. (23), can be calculated. To use the presented calibration model for polarized observations with GRAVITY, we put all the information into a small python package VLTIPOL, which is publicly available³. The package includes tools for reading out all the necessary information from a file header and calculating the Mueller matrix of a specific observation from these observations. This makes the polarimetric observing mode easy to calibrate and fully available to the community.

9 Differential effects

Apart from the absolute calibration of IP, another important question is whether the differential birefringence between the telescopes causes errors in astrometric measurements. This topic was initially addressed by Lazareff et al. (2014a). However, in their analysis, they used small random perturbations of the telescopes to estimate the phase error and get the best alignment of the optical components in GRAVITY. We could then extend their analysis, as we did not have to work with random perturbations but have the measurements of the IP for the four UTs.

9.1 Conventions – Jones formalism

For the analysis of differential effects between the telescopes, we needed a description of the propagated phase, which is not possible with the Stokes formalism. For this case, we instead used the Jones formalism, which ultimately allowed us to estimate the phase errors introduced by differential birefringence. In this formalism, the state of polarization of an electric field is described by a complex Jones vector: $$j=\left(\begin{array}{c}{E}_x\\ E_y\end{array}\right)=\left(\begin{array}{c}{A}_x\cdot e^{i{\varphi }_x}\\ A_y\cdot e^{i{\varphi }_y}\end{array}\right).$$(31)

A change in radiation is again described by a matrix, the 2x2 complex Jones matrix: $$j_{\text{out }}=J\cdot j_{\text{in }}\text{, }$$(32)

with the input Jones vector jin and the output vector j_out. One of the main disadvantages of the Jones formalism is that Jones vectors always represent fully polarized light. To be able to deal with partially polarized light, one has to use the Hermitian coherence matrix of the electric field (see Born & Wolf 1999): $$C=\left(\begin{array}{cc}\langle E_x{E}_x^{*}\rangle & \langle E_x{E}_y^{*}\rangle \\ \langle E_y{E}_x^{*}\rangle & \langle E_y{E}_y^{*}\rangle \end{array}\right),$$(33)

with E^* being the complex conjugate of the field. The coherence matrix has real values on the diagonal elements, corresponding to the total intensity in the x and y directions. The trace of the matrix gives the full intensity of the field: $$I=\mathrm{Tr}(C)=\langle E_x{E}_x^{*}\rangle +\langle E_y{E}_y^{*}\rangle .$$(34)

The off-diagonal elements are complex and describe the correlation between the x and y components of the electric field.

The DOP of an electric field described by a coherence matrix is (Born & Wolf 1999; Gil 2004) $$\text{DOP}=\sqrt{\frac{2\cdot \mathrm{Tr}\left(C^2\right)}{\mathrm{Tr}{(C)}^2}-1}.$$(35)

Applying Eq. (32) to the electric field, one can see that the Jones matrix can be applied to the coherence matrix in the following way (Hamaker 2000): $$C_{\text{out }}=J\cdot C_{\text{in }}\cdot J^{*}\text{. }$$(36)

The coherence matrix as described by Born & Wolf (1999) was then used to describe the effects of polarization in an interferometer by Hamaker et al. (1996), Hamaker (2000), and Smirnov (2011). For the combinations of two telescopes of an interferometer, m and n, the coherence matrix can be written as follows: $$C_{m,n}=\left(\begin{array}{cc}\langle E_{m,x}{E}_{n,x}^{*}\rangle & \langle E_{m,x}{E}_{n,y}^{*}\rangle \\ \langle E_{m,y}{E}_{n,x}^{*}\rangle & \langle E_{m,y}{E}_{n,y}^{*}\rangle \end{array}\right).$$(37)

The response of an interferometer to an electromagnetic signal E depends on the polarization properties of the two telescopes in the baseline. When one describes the polarization properties with Jones matrices J_m and J_n the coherence matrix propagates as follows: $$\begin{array}{l}{V}_{m,n}=2\langle J_m\cdot E_m{\left(J_n\cdot E_n\right)}^H\rangle \\ =2J_m\cdot \left(\begin{array}{cc}\langle E_{m,x}{E}_{n,x}^{*}\rangle & \langle E_{m,x}{E}_{n,y}^{*}\rangle \\ \langle E_{m,y}{E}_{n,x}^{*}\rangle & \langle E_{m,y}{E}_{n,y}^{*}\rangle \end{array}\right)\cdot J_n^H\\ =2J_m\cdot C_{m,n}^{in}\cdot J_n^H.\end{array}$$(38)

Here, ^H indicates the hermitian matrix. This approach of describing the influence of polarization effects on the interfering electric fields is taken from the radio interferometer measurement equation (RIME), which was introduced by Hamaker et al. (1996). This quantity V_m,n in the RIME formalism has different names across the literature. It has been defined as a “coherency matrix” by Hamaker (2000), a “visibility matrix” by Smirnov (2011), and a “cross-coherence matrix” in Lazareff et al. (2014a); here we adopt visibility matrix as the matrix indeed contains the complex visibilities (as outlined by Smirnov 2011). The visibility matrix describes the correlation of two electric fields ${E^{\prime }}_i$ , which have been modified by the properties described by the Jones matrix: ${E^{\prime }}_i=J_i\cdot E_i$ and the individual elements describe the correlation in two individual feeds, in our case for an orthonormal xy-basis.

As shown by Smirnov (2011) one can extend the RIME formalism to include further instrumental effects such as bandwidth smearing or can take an arbitrary brightness distribution into account, which will ultimately lead to a formulation of the Van Cittert-Zernike theorem out of the RIME formalism. However, in our case, we wanted to use this approach to describe the effect that two nonideal telescopes have on the interferometric signal and did not include any further effects. We therefore only used the formalism as given in Eq. (38). Even more so, we could directly take into account the fact that GRAVITY can only measure two linear polarization states ata time. The two polarization states are split by a Wollaston prism, which results in an angle of 90° between the two states. When we aligned our coordinate system with one of the two states, we could directly use the x and y components of Eq. (38). The useful quantities are then $V_{m,n}^{x,x}$ and $V_{m,n}^{y,y}$

As indicated before, these two quantities describe the complex visibility at the detector level of a (potentially polarized) point source, with only the polarization effects in the two light paths taken into account. From the complex visibility, we could then calculate more useful quantities, such as the correlated flux, $$F_{m,n}^x=\left|V_{m,n}^{x,x}\right|,$$(39)

and the total photometric flux, given by the sum of the flux arriving at the individual telescopes: $$I_{m,n}^x=\frac{1}{2}\left(V_{m,m}^{x,x}+V_{n,n}^{x,x}\right).$$(40)

With these two quantities, we could then calculate the maximum visibility or fringe contrast, which is just the fraction of flux that is coherent: $$v_{m,n}^x=\left|V_{m,n}^x\right|/I_{m,n}^x.$$(41)

Since at this point we are ignoring all effects except the IP and birefringence in the two light paths, this quantity gives the loss in coherence due to the polarization effects. When the two light paths are identical v_m,n equals one, and it reduces when we lose coherence due to differential polarimetric effects. We look at the coherence loss in the VLTI in Sect. 9.3.

Similarly, we could also calculate the fringe phase: $${\Phi }_{m,n}^x=\mathrm{Arg}\left(V_{m,n}^{x,x}\right),$$(42)

where Φ_m,n equals zero for identical light paths and deviates from zero otherwise. A deviation from zero fringe phase will be corrected by the fringe tracker and will not impact the coherence. However, when the science and fringe tracking target have a different polarization state, they will also have a different fringe phase. This will add an error to the measured science phase as it is referenced to the fringe tracking phase. The phase error for typical GRAVITY targets is studied in Sect. 9.5.

The equations shown here for the x component are the same for the y component. With this concept, we have everything in hand to calculate the interferometric response to polarized targets, taking into account the polarization properties of the individual telescopes. We continued fitting the VLTI model to get the Jones matrix of each telescope and then studied the differential effects and their impact on observations.

9.2 Fitting individual telescopes

Similar to the approach in Sect. 3.2, we again fit the model to the calibration data. But this time, we treated each telescope individually. The main concept of the model does not change: we reused the model as described in Eq. (20). In this case, however, we exchanged the Mueller matrix with a Jones matrix. Using Fresnel calculus, the Jones matrix that describes the linear retardation, as well as the reflection, can be written as follows: $$J=\frac{r_p}{2}\left(\begin{array}{cc}{r}_s\cdot \exp (i\delta )& 0\\ 0& r_p\end{array}\right).$$(43)

All our measurements are given as normalized Stokes vectors. One can convert them into Jones vectors with the following formula: $$j=\frac{1}{\sqrt{2}}\left(\begin{array}{c}\sqrt{1+Q}\\ \sqrt{1-Q}\cdot \exp (-i\arg (U+iV))\end{array}\right)=\left(\begin{array}{c}{A}_x\\ A_y\cdot e^{i\delta }\end{array}\right).$$(44)

As we did not have the full phase information in such a vector, we had to allow for an additional phase for each input state. In the end, we fit $$j_{\text{out }}=J\cdot \left(e^{i\varphi }\cdot j_{\text{in }}\right),$$(45)

where ϕ is an arbitrary phase factor that we ignored. For the Jones description of a mirror, we had to go back to three parameters per mirror group instead of the two parameters from Eq. (15). So, for each of the three mirror groups, we individually fit a reflection coefficient in the p and s directions and a phase difference between the two reflections. Furthermore, we split our dataset into individual telescopes, which divides the amount of available data for each fit by four. While the fit is very similar to the previous one in Sect. 4, we needed to refit it here as we are now treating the different telescopes individually.

The first result for the individual fits is shown in Fig. 16, where the reflectivity in s and p direction, as well as the phase difference between s and p, is shown. From this first look, we can conclude that in general, all the values are similar. This confirms the initial assumption that the birefringence in the individual telescopes is on the same order. Furthermore, the newly fitted values are close to the ones from the combined fit in Sect. 4. We see, however, some scatter in the reflection coefficients of the different telescopes. This means we have some degree of differential attenuation in the VLTI, which leads to some loss in fringe contrast. The phase difference, however, is very stable over all telescopes. This is a good sign, as it indicates that we have very little differential retardance, which would show up as a phase error in the observations. We look into both effects in more detail in the following.

9.3 Fringe contrast

With the response of each individual telescope, we can ask the question of how much the different IP in each light path influences the interferometric observations. For this, we used the matrix representation of the visibilities (Sect. 9.1). First, we looked at the loss in fringe contrast. We assumed that our fringetracking object is unpolarized. Nevertheless, IP may introduce a small DOP in the incoming light. As shown in Eq. (41) one can calculate the fringe contrast by the quotient of correlated and total flux, basically asking the question of how much of the total incoming light interferes. As the IP depends on the telescope’s position, we calculated the fringe contrast for a grid of telescope positions. This grid reaches from 0 to 90° in elevation and 0 to 180° in azimuth. In azimuth, the telescopes can rotate between 0 and 360°, but the polarization signal repeats after 180°, so it is sufficient to calculate the values in this range.

The fringe contrast for all telescope positions is shown in Fig. 17. For all baselines, the fringe loss is always well below 1%. This means we have a fringe contrast of above 99% for unpolarized fringe tracking targets. As the average on-sky fringe contrast for bright calibrators is on the order of 90%, we can conclude that the fringe loss due to polarization is insignificant.

Fig. 16 Comparison of the fitting values for the individual telescopes. The three plots show the reflectivity in the s and p directions, as well as the phase difference for each group of mirrors. The data points of all telescopes are grouped for one mirror group, and the fitting value for the combined fit is shown as a black horizontal line. rs , rp and δ are defined in Eq. (43).

9.4 Eigenvectors

Following Lazareff et al. (2014a), we used polarization eigenvectors to describe the polarization properties of the mirror train. A polarization eigenvector is defined as a linear input polarization that results in a linear output polarization. Lazareff et al. (2014a) have shown that each Jones matrix describing the VLTI has two distinct eigenvectors. They are generally not orthogonal to each other but close to orthogonal. With the now-measured Jones matrices for each telescope, we confirm both findings. To minimize the phase error Lazareff et al. (2014a) suggested that one should align one of the polarization directions, measured on the detector, orthogonal to one of the output eigenvectors.

The effect of phase errors due to the differential birefringence depends on the direction of the measured polarization on the detector in the instrument. We therefore use GRAVITY as an example for now, but the findings apply to all instruments.

As discussed in Sect. 5, the propagation of the polarization direction through GRAVITY is calibrated with the help of fibered polarization rotators. We could therefore assume that the two detector polarizations correspond to the same directions in the VLTI lab. For the proposed alignment of the eigenvector with the polarization axis on the detector, one could use the HWP. In the current instrument setup, the HWP rotates with the K-mirror (see Sect. 5), which does not achieve this alignment. In the following, we look into the phase errors due to the nonalignment of the eigenvectors and the detector and estimate how much they could be reduced with the optimal alignment.

Fig. 17 Fringe contrast in the six baselines shown for all possible telescope positions.

9.5 Phase errors

The measured visibility phase is referenced to the fringe-tracking object. This can be the same as the science target or a nearby star. As discussed in Lazareff et al. (2014a), this can lead to a measurement error if the fringe-tracking and the science target have a different polarization state. In the case of two different targets, the fringe-tracking object will most likely be a star and, therefore, unpolarized. In some cases, the fringe-tracking source could be slightly polarized due to foreground dust (see, for example, Buchholz et al. 2013) or intrinsic polarization of, for example, a dusty giant (Haubois et al. 2019). However, this should only be a few percent and is irrelevant here as the science target can have a much higher polarization. In the case of an (almost) unpolarized fringe-tracking target and a highly polarized science target, the polarimetric response of the mirror train is different for the two targets. This different IP will introduce a phase error, as the science phases are referenced to the fringe-tracking phases.

For the following tests, we assumed a science target with a linear polarization of 30%. This value is chosen, as it is a likely value for the Galactic Center supermassive black hole SgrA* in its flaring state (see, e.g., Genzel et al. 2010), which is one of the most extreme levels of near-infrared polarization known in celestial bodies. The following coherence matrix represents such a polarization state (see Eq. (33)): $$C=\left(\begin{array}{ll}0.65\hfill & 0\hfill \\ 0\hfill & 0.35\hfill \end{array}\right).$$(46)

As discussed in Sect. 9.1, the intensity of the electric field is given by the diagonal elements, and the DOP is calculated following Eq. (35), which gives a 30% polarization degree for this matrix. Following Eq. (41) the measured phase of a target is just the argument of the visibility matrix. We calculated this for the unpolarized fringe-tracking object and the slightly polarized science object and subtract the two phases from each other to take the phase referencing into account. As for the fringe contrast, we again calculated this for each baseline and each telescope position and show the results in Fig. 18. Theoretically, there is another degree of freedom, which is the orientation of the polarization vector on sky. This is given by the intrinsic polarization of the source as well as the parallactic angle. However, as this is just a rotation, it is redundant with the telescope azimuth and would shift the pattern in Fig. 18 to the left or right. We therefore ignore this for now.

As shown in Fig. 18, there is a small phase error that depends on the telescope position. On average, the error is on the order of 0.3°, with maximal values of 1.1°. However, the two polarizations show a somehow opposite pattern. If one averages the two polarizations, as one would probably do it for astrometry measurements anyway, to increase the signal-to-noise ratio, the phase error reduces. For the average value, the mean phase error is 0.2°, with a maximum value of 0.8°. In the simplest case, a single point source, the phases relate to the position on sky with the following formula: $$\Phi =2\pi s\cdot B/\lambda ,$$(47)

where Φ is the measured phase, s the measured position on sky and B the baseline length. Inverting this formula and using a baseline length of 100 m and a wavelength of 2.2 μm a phase error of 0.8° corresponds to an astrometric error of 10 μas.

It was already shown by Lazareff et al. (2014a) that this error can be improved if the output eigenvector of the telescope is aligned with the axis of the polarization measurement on the detector. We can confirm these previous results that each telescope always has two eigenvectors, which are roughly, but not exactly, 90° apart. If one aligns one detector polarization with one eigenvector, the astrometric error of this measurement drops to 0. However, the second polarization still shows a significant phase error. From simulating all the different options, we find that the lowest overall phase error can be achieved if one uses the average of the two eigenvector angles: $${\overline{\varphi }}_{EV}=\frac{1}{2}\left({\varphi }_{EV1}+{\varphi }_{EV2}-\frac{\pi }{2}\right).$$(48)

We could align the detector polarization with this vector by rotating the HWP by this angle. If we do so for each telescope individually, we reach the phase errors as shown in Fig. 19. One can see that the phase error has been reduced in comparison without the HWP rotation in Fig. 18. The mean phase error in this case is significantly decreased to 0.1° with maximum values up to 0.9°. Again, we could further improve this by averaging the two polarizations to mean values below 0.1° and a maximum error of 0.7°. We therefore see that while such an alignment improves the situation, it is only a small improvement if we work with the mean phase. In Fig. 19, the most significant values can be seen in all baselines with UT3. This is due to the fact that for UT3 the two eigenvectors are less orthogonal than for the other telescopes. Averaging the two eigenvector angles adds a slightly higher phase error than for the other baselines.

To reach the smallest possible phase error, as shown in Fig. 19, one would need to track the eigenvector with the HWP. This tracking angle is a pure telescope property and, therefore, does not depend on the polarization on sky or the parallactic angle. It could be implemented as a look-up table based on the derived values. The angles are very similar for the four telescopes but not identical. This means one must align each telescope individually to get the smallest phase error. This might not be desired as then each telescope would have a different orientation of the polarization axis on the detector. Furthermore, tracking of the HWP could introduce systematic effects into the metrology measurement. We find that the phase error without tracking is comparably small and does not dominate over other systematic effects.

In conclusion, one does not need to align the detector with the eigenvectors as long as one uses the mean phase for the highest resolution astrometry but the individual polarizations for the polarization measurement and the imaging.

Fig. 18 Error in the visibility phases due to differential birefringence for all telescope positions. The left two columns show the error for the first polarization P1 and the right two for the second polarization P2.

Fig. 19 Same as Fig. 18, but this time the phase error if the GRAVITY HWPs track the eigenvector of the individual telescopes.

10 Conclusion

In this paper we have presented the first complete polarization study of the VLTI and GRAVITY. As expected, both the observatory and the instrument itself show polarization effects. In both cases, we have characterized the effects and built a calibration model to calibrate for the instrumental effects in observations. We have outlined the necessary steps in observation, data reduction, and calibration to execute polarimetric measurements with GRAVITY. We then demonstrated the capabilities by remeasuring the polarization properties of the Galactic Center star IRS 16C, which is in excellent agreement with the literature.

We have also shown that differential birefringence between the light paths of the VLTI UTs is not a dominant error source; the four light paths were constructed with great care to minimize differential birefringence. For a typical observation of a calibrated source, our studies have shown a phase error due to differential birefringence of below 1°. Even for the extreme case with a 30% polarized science target, this results in an astrometric error of only around 10 μas. This error can be further reduced when using the average of the two polarizations for astrometry and the individual signal for polarimetry. The fringe contrast in such a case is only reduced by around 1%.

We have therefore demonstrated that observations with GRAVITY do not suffer from strong effects due to birefringence and that GRAVITY can be used for polarimetric observations. This can be done in very different ways. With IRS 16C, we have shown that GRAVITY can measure the polarization quantities of even slightly polarized targets with very good precision. One can also study the temporal evolution of polarized targets, as done in GRAVITY Collaboration (2018b, 2020b) for a bright SgrA* flare. Furthermore, the polarization information can be extracted not only for the full intensity, but also for the measured intensity in each spectral channel. One can therefore even map polarization changes over the spectral range of GRAVITY. As all of this comes together with the unprecedented resolution of GRAVITY, opening up a wide range of possibilities for polarimetry in the near-infrared.

Acknowledgements

We are very grateful to our funding agencies (MPG, ERC, CNRS [PNCG, PNGRAM], DFG, BMBF, Paris Observatory [CS, PhyFOG], Observatoire des Sciences de 1’Univers de Grenoble, and the Fundação para a Ciência e Tecnologia), to ESO and the Paranal staff, and to the many scientific and technical staff members in our institutions, who helped to make GRAVITY a reality. F.W. has received funding from the European Union’s Horizon 2020 research and innovation programme under grant agreement No 101004719. A.A., and P.G. were supported by Fundação para a Ciência e a Tecnologia, with grants reference SFRH/BSAB/142940/2018, UIDB/00099/2020 and PTDC/FIS-AST/7002/2020.

Appendix A Modeling of the coated mirrors

For the coated mirrors, we had to modify the Fresnel equations introduced in subsection 3.2 to account for an additional thin layer on top of the bare metal mirror. This is usually done with techniques from ellipsometry, a method for determining the thickness of a thin film by measuring its effect on the polarization of a reflected beam. A more complete overview of this topic can, for example, be found in Jellison (1999) or in Chapter 29 of Goldstein (2003).

Fig. A.1 Measured and modeled phase difference for the protected silver mirrors. The black points show the measured data, and the solid red line shows the model with the thickness that best fits the data. For illustration, two more models with higher thicknesses are shown in dashed and dotted lines.

For this analysis, we started with equations 12 and 13, which give the reflectance for a single layer. Here we could not assume that one of the refractive indices is equal to one, but we had to assume that both indices can be arbitrary. By using Snell’s law, we could modify equations 12 and 13 to the following: $$R_s=\frac{n_2\cos \left({\Theta }_1\right)-n_1\cos \left({\Theta }_2\right)}{n_2\cos \left({\Theta }_1\right)+n_1\cos \left({\Theta }_2\right)},$$(A.1) $$R_p=\frac{n_1\cos \left({\Theta }_1\right)-n_2\cos \left({\Theta }_2\right)}{n_1\cos \left({\Theta }_1\right)+n_2\cos \left({\Theta }_2\right)},$$(A.2)

We then needed to evaluate these equations twice, once for the transition from air to the thin film (R_1,s and R_1,p) and once for the transition from the thin film to the base layer (R_2,s and R_2,p). With those four values, we could then get the reflectance for the full mirror with the thin coating: $$R_{f,s}=\frac{R_{1,s}+R_{2,s}\cdot \exp (-ib)}{1+R_{1,s}\cdot R_{2,s}\cdot \exp (-ib)},$$(A.3) $$R_{f,p}=\frac{R_{1,p}+R_{2,p}\cdot \exp (-ib)}{1+R_{1,p}\cdot R_{2,p}\cdot \exp (-ib)},$$(A.4) $$b=\frac{4\pi d{n}_f\cos {\Theta }_f}{\lambda },$$(A.5)

where d is the thickness of the layer, n_f the refractive index of the layer, and Θ_f the angle within the layer, which is calculated with Snell’s law from the incidence angle: sin Θ_i = n_f sin Θ_f.

The quantities of the diattenuation and the phase difference are then calculated as before.

We used the thin film formulas for all silver mirrors in the VLTI, as they have a protective coating. We have measurements for the phase difference of one of the mirrors available (F. Delplancke, private communication), but we do not know specifics about the coating itself. To model the mirrors, we assumed a protective coating out of Al₂ O₃, a reasonably typical coating for silver mirrors. To estimate the thickness of the coating, we fit the calculated phase difference ∆ = arg R_f,p − arg R_f,s to the measured data and fit for the thickness. We find that a thickness of 210 nm represents the data reasonably well, as shown in Fig. A.1. We adopted this value for all silver mirrors in the initial model, but it is unclear if it is a good representation for all mirrors. However, as the diattenuation and retardance are fitted to data in a later step, it should not influence the calibration model.

Appendix B Wavelength dependence of the model

Fig. B.1 Effect of different wavelengths on the polarization model of the VLTI light path. The top panel is similar to Fig. 7, showing the measured polarization angle (corrected for field rotation) for all four telescopes and all linear input states in gray. The model is shown in four wavelengths (1900, 2000, 2250, and 2500 nm). As the model repeats after 180°, only the range from 0 to 180 degrees is shown. The bottom panel shows a scaled-up view of the model for state 4 without the data.

The data we used to calibrate our model was taken with a laser at a wavelength of 1908 nm. This is slightly lower than the science wavelength of GRAVITY, which is between 2000 and 2500 nm. As the refractive index of the mirrors we model changes with wavelength, the difference in wavelength introduces a small error into the model. However, for the used materials, the change in refractive index over the wavelength range is ±(0.1 + 2i). This is comparably small and should not introduce a big error. To test the effect, we modeled the VLTI light path at four different wavelengths: at 1908 nm, where the laser is, and at the minimum, mean, and maximum wavelength of the GRAVITY spectrum (2000, 2250, and 2500 nm). The model at different wavelengths is shown in Figure B.1. The difference between the four models is small, with a mean absolute difference of 0.3 between the laser wavelength and the mean science wavelength (2250 nm) and a maximum difference of below 1°, depending on the telescope position. This is in the same order as the repeatability of the experimental input states and smaller than the uncertainty of the fitted phase shifts (section 4). We conclude that the wavelength of the laser is not a dominant error in our model. Especially for longer observations, it will average out, and for shorter observations, it will add an uncertainty of below one degree, dependent on the telescope position.

Appendix C Mueller matrix of GRAVITY

The measured Mueller matrices for GRAVITY are as follows:

For the K-mirror, $$M_{KM}=\left(\begin{array}{cccc}1& -0.0146& 0& 0\\ -0.0146& 1& 0& 0\\ 0& 0& -0.8949& -0.4461\\ 0& 0& -0.4461& -0.8949\end{array}\right).$$(C.1)

For the HWP, we assumed a perfect matrix: $$M_{HWP}=\left(\begin{array}{cccc}1& 0& 0& 0\\ 0& 1& 0& 0\\ 0& 0& -1& 0\\ 0& 0& 0& -1\end{array}\right).$$(C.2)

For the remaining optics in off-axis mode, $$M_{off-axis}=\left(\begin{array}{cccc}1& 0.0036& 0.0250& 0.0110\\ -0.0237& -0.7143& -0.6685& -0.1674\\ 0.0261& -0.6743& 0.7041& 0.1918\\ 0.0294& 0.4122& 0.2723& -0.8742\end{array}\right)\text{, }$$(C.3)

and in on-axis mode, $$M_{on-axis}=\left(\begin{array}{cccc}1& 0.0612& 0.0496& 0.0141\\ -0.0901& -0.7457& -0.7160& -0.0836\\ 0.0870& -0.6884& 0.7768& 0.1066\\ -0.0246& 0.3469& -0.0438& -0.8325\end{array}\right)\text{. }$$(C.4)

References

All Tables

All Figures

	Fig. 1 Schematic layout of the VLTI platform. The positions of the four UTs are shown as big circles, and the possible stations for ATs as small circles. The delay lines and the VLTI Lab are indicated in the center of the platform.
In the text

Fig. 2 Schematic model of the VLTI light path. The light path is shown for two UTs, including GRAVITY, in red and blue. The rounded arrows indicate the possible rotations in the VLTI – for changes in azimuth (Az), elevation (El), and parallactic angle (Pa) – and in GRAVITY, at the field de-rotator and the half-wave plate (HWP). The straight arrows indicate the movement of the delay lines. The location of the star separator is only shown for the blue beam but is in the same place for both telescopes. The fiber optics and elements within the GRAVITY beam combiner are only indicated here and further explained in Sect. 5.

In the text

Fig. 3 Simplified version of the VLTI light path from Fig. 2 to show the modeling and the experimental setup. The black rectangles show where the laser is launched and where the polarimeter is mounted. The names of the mirrors used in the text are given. The color of the mirror number shows the grouping that was used for the fitting. Gray mirrors are not fitted in our calibration model. Each colored group is located in one common plane: M4−M8 are in one vertical plane, and M10−M18 are in one horizontal plane.

In the text

	Fig. 4 Diattenuation (top) and relative phase shift (bottom) of the individual VLTI mirrors in our model, shown as a function of the reflection angle for a wavelength of 2250 nm. The worst offenders are labeled in both plots. The colored lines show the continuous values for the three used materials: gold, silver, and aluminum. M9 does not lie on a line, as it is a dichroic mirror.
In the text

	Fig. 5 Properties of the output polarization for the input state S = (1, 0.3, 0.3, 0)T and all possible telescope positions. The plots show, from left to right, the DOP, the DOLP, and the polarization angle. The input state has a DOP and DOLP of 42.5% and a polarization angle of 22.5°. These values are shown as triangles in the color bars.
In the text

	Fig. 6 Expected DOP for unpolarized light at all telescope positions.
In the text

Fig. 7 Measured polarization angles for the different telescope positions, only corrected by the geometric rotation of the field. The four panels show the data for each linear input state, with the different colors showing the data from the four UTs. The data from 50 and 70° elevations are shown together. The gray horizontal line shows the input values for the polarization angles, the solid black line the value expected from the model for 70° elevation, and the dotted black line the model for 50° elevation.

In the text

Fig. 8 Measured DOLP for the different input states. Left and middle: linear input states, where the input DOLP is 100%. Right: circular states, where the input is fully circularly polarized, so 0% linear polarization. The data points are shown in different colors for the different states. The data from 50 and 70° elevations are shown together. The lines in the same color show the model prediction for each state (solid line for 70° elevation and dotted for 50° elevation). The dashed gray line shows the input values in all plots.

In the text

Fig. 9 Test data calibrated with the fitted polarization model. The data are the same as in Figs. 7 and 8 but now calibrated. In both plots, the data are shown as dots, and the input states as dashed lines. The different colors show the different input states. The left plot shows the polarization angle, and the right plot shows the DOLP. Since all states are fully polarized, the recovered degree of circular polarization is (100 - DOLP).

In the text

	Fig. 10 Schematic side and front view of the optical design of the GRAVITY fiber coupler. From GRAVITY Collaboration (2017).
In the text

	Fig. 11 Optical design of the spectrometer. Most important for the polarimetric mode is the Wollaston prism, which can be moved in or out of the light path. From GRAVITY Collaboration (2017).
In the text

Fig. 12 Sketch of the two polarimetric measurements for GRAVITY. Left: in the first measurement, the light from the metrology laser is used and measured with a polarimeter. Right: the light from the calibration unit is used and recorded on the spectrometer. In both cases, the light direction is indicated by arrows. The light path is the same as shown in Fig. 10, with the reflection for the Tip-Tilt and Piston control just shown as a change in direction in the light path.

In the text

	Fig. 13 Measurement of the polarization in GRAVITY with a rotating HWP. The points show the data points for Stokes Q (red) and Stokes U (gray). The black lines show the corresponding values for an ideal optical element.
In the text

Fig. 14 Measured polarization with GRAVITY in the different observing modes. The left column shows the measurement in the off-axis mode, and the right column in the on-axis mode. In the top row, the linear polarization filter is used for the input light source; in the bottom, it is not. In all plots, the Stokes Q data points are shown in red and orange and the Stokes U in gray. The data are the average over all four GRAVITY beams. The results from the fitted model are shown in black lines.

In the text

	Fig. 15 Polarization measurements of IRS 16C between 2019 and 2022. The top panel shows the polarization angle for each night, and the bottom panel shows the polarization degree. The average over all nights is shown as a horizontal line
In the text

	Fig. 16 Comparison of the fitting values for the individual telescopes. The three plots show the reflectivity in the s and p directions, as well as the phase difference for each group of mirrors. The data points of all telescopes are grouped for one mirror group, and the fitting value for the combined fit is shown as a black horizontal line. rs , rp and δ are defined in Eq. (43).
In the text

	Fig. 17 Fringe contrast in the six baselines shown for all possible telescope positions.
In the text

	Fig. 18 Error in the visibility phases due to differential birefringence for all telescope positions. The left two columns show the error for the first polarization P1 and the right two for the second polarization P2.
In the text

	Fig. 19 Same as Fig. 18, but this time the phase error if the GRAVITY HWPs track the eigenvector of the individual telescopes.
In the text

	Fig. A.1 Measured and modeled phase difference for the protected silver mirrors. The black points show the measured data, and the solid red line shows the model with the thickness that best fits the data. For illustration, two more models with higher thicknesses are shown in dashed and dotted lines.
In the text

Fig. B.1 Effect of different wavelengths on the polarization model of the VLTI light path. The top panel is similar to Fig. 7, showing the measured polarization angle (corrected for field rotation) for all four telescopes and all linear input states in gray. The model is shown in four wavelengths (1900, 2000, 2250, and 2500 nm). As the model repeats after 180°, only the range from 0 to 180 degrees is shown. The bottom panel shows a scaled-up view of the model for state 4 without the data.

In the text