© The Authors 2024

1 Introduction

In the last two decades, many exoplanets have been detected using radial velocity (RV) measurements, transit photometry and spectroscopy, direct imaging, and other techniques and a lot has been learned about exoplanet frequency, orbital parameters, masses, radii, composition, and surface structures. However, almost no progress has been made in the observational characterization of the surfaces and atmospheres of cold exoplanets, including potentially habitable planets. Cold planets are very faint sources in the mid-infrared, and they produce only scattered light in the optical to near-infrared range. Therefore, their investigation with direct imaging requires very deep contrast limits and, as of yet, there exists no successful detection of a cold exoplanet with direct imaging.

High-contrast imaging is a very successful technique for the investigation of young, self-contracting, self-luminous giant planets that are hot (T_eff ≈ 1000 K) and rather luminous in the near-infrared range (Nielsen et al. 2019; Vigan et al. 2012, 2017; Bowler 2016). Typical examples for current near-infrared techniques provide contrasts of I_p/I_★ ≈ 10^−3.24 for an angular separations of about ρ ≈ 0.5″ to ≈ 10⁻⁴ for ρ > 2″ for the L′ band (3.8 μm) with the Very Large Telescope (VLT) instrument NACO (Cugno et al. 2023) or of I_p/I_★ ≈ 10⁻⁵ for ρ ≈ 0.5″ to ≈10⁻⁶ for ρ > 2 for the H band (1.6 μm) with VLT/SPHERE/IRDIS (Langlois et al. 2021).

Cold planets are much fainter and harder to detect because they have no strong intrinsic energy source and only reprocess stellar irradiation. Therefore, their luminosity is proportional to the irradiation $L_{\mathrm{p}}\propto L_{\star }{R}_{\mathrm{p}}^2/d_{\mathrm{p}}^2$, where R_p is the planet radius and d_p the orbital separation. The planet radiation is partly emitted as scattered light at the same wavelengths as the stellar emission, and partly as thermal radiation peaking for cold planets T_eq ≤ 300 K at 10 μm or even longer wavelengths. The factor $R_{\mathrm{p}}^2/d_{\mathrm{p}}^2$ is very small, only ≈2.3 × 10⁻⁷ for a Jupiter-sized planet (R_p = R_J) at a separation of d_p = 1 au. The contrast is less demanding for a smaller physical separation dp where, however, the inner working angle for the high-contrast imaging becomes an issue (e.g. Milli et al. 2013). Furthermore, the achievable contrast is better for larger angular separations ρ, but for large d_p the planet signal is weak. These two conditions – small d_p and at the same time large ρ − limit the direct imaging search to planets in nearby systems within 5 pc to 10 pc (Lovis et al. 2017; Kasper et al. 2021).

This paper presents high-contrast imaging observations of ε Eri taken in 2019 and 2020 with VLT. This is a continuation of the SPHERE RefPlanets guaranteed time observation (GTO) programme which uses the SPHERE instrument (Beuzit et al. 2019) with the Zurich IMaging POLarimeter (ZIMPOL) subsystem (Schmid et al. 2018) for the search of polarization signals from the scattering of stellar light by planets. Thereby, polarimetry serves as a powerful differential imaging technique because the polarization signal can be distinguished in the point-spread-function (PSF) halo from the unpolarized light of the much brighter star (Schmid et al. 2006a). The expected fractional polarization of a planet is at the level of about p_p ≈ 5%-50%, it depends on the orbital phase angle α, and it strongly constrains the properties of planets as described in Seager et al. (2000), Stam et al. (2004), Schmid et al. (2006a), and Buenzli & Schmid (2009).

The first results of the RefPlanets programme are presented in Hunziker et al. (2020) who observed six targets: α Cen A, α Cen B, Sirius A, Altair (α Aql), ε Eri, and τ Ceti. This programme was executed similar to a blind search and did not target known giant planets or a potential planet candidate. Typically, integration times between 1.5 and 3.5 h were obtained and the achieved 5σ signal-to-noise ratio (S/N) contrast limits for the polarized flux contrast is about C_P = (p_p · I_p)/I_★ ≈ 10⁻⁷ at an angular separation of ρ = 0.5″ and about ≈10⁻⁸ at ρ = 1.5″. Thus, Jupiter-sized planets were only within reach for the nearest targets α Cen A and α Cen B, where 1 au corresponds to ρ = 0.7″. The indicated limits also require that the planet has a high albedo, that it produces a scattering polarization of about 20%, and that it is located at the right orbital phase during the observations. Therefore, a blind search can easily miss a detectable planet if the observations are carried out during an unfavourable orbital phase. For other systems, which are further away than α Cen, even a planet with R_p ≈ R_J would have been too faint to be detected. An important result of this programme was the demonstration that the contrast limits at separations larger than 0.5″ reach the photon noise limit and will therefore improve with the square root of the integration time (Hunziker et al. 2020).

The ε Eri observations presented in this work were collected to achieve a point source contrast limit of C_P < 10⁻⁸, which is significantly deeper when compared to the study of Hunziker et al. (2020). This requires the combination of more than 30000 integrations of 3 or 5 s for a total exposure time of more than 38 h. For this, one has to combine data from different runs and consider that a possible planet moves on its orbit around the star by several pixels per week or many pixels per month. Therefore, one needs to combine the whole time series with a prediction for the Keplerian orbit of possible targets in the post-processing as described in Nowak et al. (2018), Le Coroller et al. (2020) and Dallant et al. (2023).

The system ε Eri is very interesting for such a deep search and performance test with SPHERE/ZIMPOL, because of strong indications for the presence of a RV planet with a semi-major axis of about 3.5 au or a separation of about ρ ≈ 1″ (e.g. Llop-Sayson et al. 2021). Further, ε Eri shows strong thermal dust emission in the infrared including a component peaking at 20 μm from warm dust (Backman et al. 2009), which should produce an extended polarization signal from light scattering by dust within the ZIMPOL field of view. As ε Eri is the nearest and brightest single solar-type star, it is important for the investigation of extra-solar planetary systems (Backman et al. 2009).

This paper is organized as follows. Section 2 summarizes the parameters for the planet candidate ε Eri b and the dust near the star from the literature and provides predictions for the possible polarization signal from the planet and dust. In Sec. 3 we describe the observations, the data reduction including the specific post-processing procedures for the search of a faint point source, and signatures from circumstellar dust. In Sects. 4 and 5, we present the final detection maps and explain how we derived the sensitivity limits of our planet search and the search for dust scattering. In Sec. 6, we discuss our findings and give our conclusions in Sec. 7.

2 ε Eridani

The star ε Eri is a single star at a distance of 3.2 pc (Gaia Collaboration 2020), with a spectral type of K2V, T_eff ≈ 5040 K, a mass of 0.82 M_⊙, a luminosity of L_s = 0.32 L_⊙ (e.g. Baines & Armstrong 2012), and an apparent brightness of m_V = 3.7 mag, m_R = 3.0 mag, m_I = 2.5 mag. It is a very active solar type star with a high level of chromospheric activity, with activity cycles of about 13 and 3 yr (Metcalfe et al. 2013) and a rotation period of 11.67 days (Donahue et al. 1996). The system is young and might belong to the 500 Myr old Ursa Major association (Fuhrmann 2004), but its spatial motion is one of the most deviant compared to the mean motion of the group. Therefore Janson et al. (2008) assume a large age range of 200–800 Myr.

Significant RV variations ±20 m s⁻¹ were found by Campbell et al. (1988). They are at least partly caused by chromospheric activity, but there seems to exist also a periodic RV signal from the reflex motion introduced by a giant planet. Initial estimates for the RV period were 10 yr (Walker et al. 1995), and then 6.9 yr (Cumming et al. 1999; Hatzes et al. 2000). Longer time series established quite firmly an orbital period of about 7.3 yr (Anglada-Escudé & Butler 2012; Mawet et al. 2019; Llop-Sayson et al. 2021) while Zechmeister et al. (2013) did not find this periodicity. In Mawet et al. (2019) and Llop-Sayson et al. (2021) it is strongly suggested that this signal is caused by a planet on an orbit with low eccentricity e < 0.1 producing a RV semi-amplitude of about 11 m s⁻¹ (Anglada-Escudé et al. 2010; Mawet et al. 2019; Llop-Sayson et al. 2021). However, with the available data it is hard to rule out the possibility of another stellar activity cycle introducing such a RV signal.

We adopt for this work the interpretation of a RV-planet with a period of 7.3 yr (2671 days) and base much of our signal predictions and parts of the data interpretation on this assumption. This solution predicts a_p = 3.53 au for the semi-major axis of the planet orbit and that the planet was further away from us than the star around T = 2019.5 (JD 2458 666), at RV-phase φ = 0.5 according to Llop-Sayson et al. (2021) while quadrature phase was around 2021.3. Our data were taken between these two dates which correspond roughly to an orbital phase of ϕ = 0°, when we expect maximum intensity for a reflecting planet, and ϕ = 90° when the fractional polarization of the planet should be highest (see Sec. 2.1). These orbital phase estimates could be affected significantly by RV uncertainties introduced by the chromospheric activity.

The presence of a planet in ε Eri is supported by astrometric measurements taken with the HIPPARCOS satellite (Reffert & Quirrenbach 2011), the HST Fine Guidance Sensor (Benedict et al. 2006), and by the search of proper motion anomalies combining older astrometry with results from the Gaia early data release 3 (e.g. Kervella et al. 2022; Benedict 2022; Makarov et al. 2021). All these studies indicate for ε Eri systematic deviations from a constant proper motion vector but the remaining uncertainties are still quite large. The results seem to be compatible with the presence of a planet with mass of ≈1 M_J as measured by radial velocity, on a prograde orbit (N over E), and that it should be located about north-east of the star for our observations from 2019 and 2020 (Benedict 2022). It is expected that accurate astrometric data for ε Eri, for example from the Gaia mission, will provide in the near future, strong, new constraints on the orbit and the mass of the planet.

Radial velocity and astrometric measurements have been combined to obtain more detailed parameters for ε Eri b (Benedict et al. 2006; Reffert & Quirrenbach 2011; Llop-Sayson et al. 2021; Benedict 2022). However, these studies are based on extra assumptions, in particular it is not considered that additional planets in the system could contribute to the measured reflex motion. Therefore, we use the derived orbital parameters only as possible values.

If the ε Eri b planet is on a circular orbit, then the measured RV semi-amplitude K = 10.34 m s⁻¹ and period P = 2671 days from Llop-Sayson et al. (2021, Table 3) constrain the radius of the stellar orbit as = (K · P)/(2π · sin(i)) = 2.54 × 10⁻³ au/sin(i) = 0.78 mas/sin(i) and this defines the minimum mass for the planet m_p sin(i) = (a_s/a_p) · M_s = 0.651 M_J. For small i the planet mass and therefore the stellar reflex motion measurable by astrometry would be significant larger than the minimum values.

An estimate for the equilibrium temperature for ε Eri b, using a_p = 3.53 au, L_s = 0.32 L_⊙, and Bond albedo A_B = 0.3 gives T_eq ≈ (L/L_⊙)^1/4 (a/au)^−1/2 T_Earth = 112 K similar to Jupiter. If there is no strong internal energy source, then the scattered intensity will dominate in the visual and near-infrared range up to at least 3 μm and the expected planet to star intensity contrast is to first order wavelength independent at a level of about Δm ≈ 20 mag to 21 mag. At 10 μm and perhaps also 5 μm the planet brightness could be dominated by thermal radiation, from Kelvin-Helmholtz contraction in case ε Eri is younger than estimated.

Several high-contrast imaging searches for faint point sources around ε Eri have been carried out and we give here an incomplete selection of reported intensity contrast limits from the literature for a separation of 1″: about Δm ≈ 10 mag in the N-band at 11 μm (Pathak et al. 2021), Δm ≈ 13 mag in Ms at 4.7 μm (Mawet et al. 2019), Δm ≈ 13.5 mag in Lp at 3.8 μm (Mizuki et al. 2016), Δm ≈ 14 mag in H at 1.6 μm (Janson et al. 2007), Δm ≈ 14 mag in the RI-band at 0.75 μm (Hunziker et al. 2020). These contrast limits are more than 6 mag away from the detection of the reflected intensity of a planet. For polarimetric imaging Hunziker et al. (2020) report in the RI-band a polarimetric contrast of Δm_p ≈ 18 mag which is only about 4 mag from the expected planet signal considering that the polarized intensity is about 1.5 mag fainter than the intensity (see Sec. 2.1).

The ε Eri system is also well known for its infrared excess of L_dust/L_★ ≈ 1.0 × 10⁻⁴ caused by the thermal emission of circumstellar dust (Gillett 1986; Decin et al. 2003; Backman et al. 2009). This infrared-emission is composed of emission from an outer circular ring of cold dust with an inclination of about 30^° and T_dust ≈ 50 K located at r ≈ 65 au or ρ ≈ 20″ (e.g. Greaves et al. 2014; Chavez-Dagostino et al. 2016).

Also a region of warm dust has been resolved and it extends to about r = 14 au or ρ ≈ 4″ and a temperature of about 100 K (Greaves et al. 2014). According to the SED analysis of Backman et al. (2009) and Su et al. (2017), the warm dust has two components, one at 20 au and one at a separation of around 3 au producing an infrared-bump at 20 μm from dust with T ≈ 120 K. This latter ‘warm’ dust component with an infrared excess of L_warm/L_★ = 3.3 × 10⁻⁵ (Backman et al. 2009) is most interesting for this study, because a signal of scattered light is expected in the field of view of our SPHERE/ZIMPOL observations. Therefore we also search for extended emission of polarized light in our data. Even finding only upper limits could be of interest as ε Eri is the closest debris disk system to the Sun and therefore a prototype and important test case for the theoretical modelling (e.g. Reidemeister et al. 2011; Su et al. 2017). With HST intensity data Wolff et al. (2023) could not find dust scattering outside of the inner working angle of 1′′.

2.1 Predictions for the signal of ε Eri b

The results from the RV and astrometry data are very helpful for predicting the expected signal of the light reflection from ε Eri b. The intensity contrast C_I = I_p/I_★ and the closely linked polarized intensity (or polarization) contrast C_P can be calculated according to

$$C_{\mathrm{P}}=p(\alpha )C_{\mathrm{I}}=p(\alpha )f(\alpha )\frac{R_{\mathrm{p}}^2}{d_{\mathrm{p}}^2},$$(1)

where R_p is the planet radius and d_p the star to planet separation (Stam et al. 2004). The reflectivity f(α) and the fractional polarization p(α) of the scattered light are both functions on the planet scattering angle α and depend on the wavelength.

For our simple estimate we consider circular orbits because then the factor $R_{\mathrm{p}}^2/d_{\mathrm{p}}^2$ does not change with orbital phase ϕ. Its value is 1.85 × 10⁻⁸ for R_p = R_J and d_p = 3.53 au. The relation between α and ϕ depends on the inclination i of the orbital plane

$$\alpha =\arccos (\sin (i)\cos (\varphi )),$$(2)

where ϕ = 0^° is the phase when the planet illumination as seen by the observer is maximal. For a face-on orbit (i = 0^°) the scattering angle is constant at α = 90^°. For an edge-on orbit (i = 90°), there is α = ϕ and we define α = 0° for the back-scattering configuration when f(α) is maximal and f(0°) identical to the geometric albedo of the planet, while f(180°) = 0. The p(α) phase curve has typically its maximum close to right-angle scattering p_max ≈ p(90°).

Figure 1 illustrates the expected prograde orbital motion of ε Eri b and the time of our four observing epochs. The inserted PSFs are splitted left/right and show the relative signal strengths along an orbit for the intensity in red and the polarized light in blue. The orientation and inclination of the plotted orbit is not well known and therefore tentative.

The movement of the planet is small enough to allow the combination of data from consecutive nights of an epoch. The maximum separation of the planet is 1.1″ and an orbit takes 2671 d (Llop-Sayson et al. 2021). If the orbital inclination is 0° and perfectly circular, then the orbital speed is constant: (2πr)/P = 2π · 1.1″/2671 = 2.585 mas d⁻¹ = 0.718 px d⁻¹. The projected speed is equal or slower for inclined orbits and particularly small near quadrature phase for high i.

The reflectivity f(α) and the fractional polarization p(α) of a planetary atmosphere depends on the atmospheric structure, which is hard to predict. There exist detailed model calculations for the reflected intensity and polarized intensity of giant planets (e.g. Stam et al. 2004; Bailey et al. 2018) and simple parameteric models (e.g. Seager et al. 2000; Buenzli & Schmid 2009; Madhusudhan et al. 2011). For our estimate we just pick a model for a Rayleigh scattering planet from Buenzli & Schmid (2009), with an optical depth τ_sc = 2 for a Rayleigh scattering layer, with a single scattering albedo of ω_R = 0.95, above a cloud layer with a Lambertian surface with an albedo of A_S = 1. This is the same model as shown in Fig. 1 of Hunziker et al. (2020). This model produces a polarized reflectivity p(α)f(α) of about ≈0.055 for quadrature phase ϕ = α = 90°.

In Fig. 2 we show the expected intensity contrast C_I and polarized intensity contrast C_P (Eq. (1)) for i = 78° and 30°. For inclined orbits p(α)f(α) and therefore also C_P is small or very small for orbital phases ϕ ≈ 120°−240° when the planet is closer to us than the star so that we can hardly see the illuminated hemisphere. For i ≥ 30° the favourable orbital phases are around ϕ ≈ ±70°, when the planet is further away from us than the star so that we see a substantial fraction of the illuminated hemisphere, but still near a scattering angle which is strongly polarizing. The polarized reflectivity is then at max(p(α)f(α)) ≈ 0.07 or about 25% higher than at quadrature phase. For conjunction ϕ = 0° or maximum illumination the polarized reflectivity is still high for low inclination (i ≈ 30°), but has a dip for i ≈ 40°−70°, while there is only a weak signal for i > 70°, because then the back-reflection produces only a very small scattering polarization. We show the expected contrast as a function of time: in the upper panel for the intensity contrast $C_{\mathrm{I}}=f(\alpha )\cdot R_{\mathrm{p}}^2/d_{\mathrm{p}}^2$ and in the lower image for the polarized intensity contrast C_P (Eq. (1)). The black curve illustrates the signal strength for the orbit as shown in Fig. 1 with an inclination of 78^° and the blue curve for an alternative inclination of 30^° closer to the inclination of the outer dust. The alternative curve for 30^° only shows the effect of a different illumination and scattering of the same hypothetical planet. This means that this curve does not account for the fact that the planet with an orbital inclination of 30^° has, according to RV measurements and the resulting factor m sin(i) a higher mass. This higher mass could lead to a different planetary radius which ultimately would change the signal strength overall. The dots in the colour of the curves show our four observing epochs in 2019 and 2020, while the octagons stand for earlier SPHERE/ZIMPOL observations in Hunziker et al. (2020). Trusting the available RV data, our fourth observing epoch is at the ideal timing with a large angular separation and a scattering angle close to 90° for a strong polarization signal.

The curves in Fig. 2 are only rough estimates, because the reflectivities f(α) and p(α)f(α) depend on atmospheric parameters (Buenzli & Schmid 2009). For orbits with an eccentricity of e ≈ 0.2 the separation d_p(ϕ) can be 20% larger or smaller during the orbit. Also the adopted planet radius R_p = R_J could be larger or smaller by 10% to 20% as follows from the distribution of measured giant planet radii (Thorngren et al. 2019). We conclude that the expected planet signal could be up to a factor two larger or a factor of a few smaller than the estimates given above. Circumplanetary rings similar to Saturn could also boost the signals of the reflected light by a factor of two or even more (Arnold & Schneider 2004; Dyudina et al. 2005).

Fig. 1 Illustration for a possible orbit (i = 78°) and signal strength of ε Eri b. The blue colour indicates the normalized strength of the polarized intensity and the red colour is analogously for the intensity. The four observing epochs are marked with black dots and the points in the orbit with a minimal, zero, and maximal RV are shown in grey.

Fig. 2 Planet model contrast as a function of time for the intensity (upper panel) and the polarized intensity (lower panel) for the orbit with an inclination of 78° (black curves) as in Fig. 1 and for an inclination of 30° (blue curves). The dots show our four observing epochs, while the octagon stands for earlier ZIMPOL observations in Hunziker et al. (2020).

2.2 Predictions for the scattered light from warm dust

The spectral energy distribution of ε Eri shows a multi-component infrared excess with a ‘cold’ dust component in the far-infrared and a warm dust component peaking around 20 μm. The warm component has a relative luminosity of L_warm/L_★ ≈ 3.3 × 10⁻⁵, and it is estimated that this dust is located at a separation of roughly 3 au from the star according to Backman et al. (2009) or between 2.5 and 6 au according to Su et al. (2017), depending a lot on the adopted grain properties. This dust should also produce a scattered light signal within the field of view of SPHERE/ZIMPOL which covers about r = 5 au.

Based on this infrared-emission of the warm dust, we can make estimates about the polarization Q_ϕ, which depend strongly on the adopted dust scattering properties and the disk geometry. We use for a rough estimate of Q_ϕ simple axisymmetric and optically thin models as illustrated by a disk ring in Fig. 3.

The thermal emission of the dust is radiated isotropically, while the disk integrated polarized emission Q_ϕ/I_* depends strongly on the disk inclination and the polarized scattering phase function f_φ(θ) of the dust, where θ is the scattering angle measured as angle of deflection. The function f_φ(θ) = f_I(θ) × p(θ) is described by a Henyey-Greenstein function f_I = f_HG(θ, g) with asymmetry parameter g for the distribution of the scattered intensity, while the fractional polarization has the same angle dependence as for Rayleigh scattering p(θ, p_max) = p_max(sin²(θ))/(1 + cos²(θ)), but with a scaling factor of p(90°) = p_max which accounts for the reduced scattering polarization produced by large, compact dust particles. The adopted parameters are roughly representative for the zodiacal dust in the solar system (Leinert 1975).

The integrated polarization Qϕ/I_* of an optically thin, axisymmetric disk depends on the total cross section of the scattering dust σ, and on a disk averaged polarized scattering phase function ⟨f_φ(i,g, p_max)⟩ (Schmid 2021). One can relate ⟨f_φ(i, g, p_max)⟩ to the intensity phase function for isotropic scattering ⟨f_I(g = 0)⟩ and with the total dust absorption cross section κ also to the isotropically emitted infrared excess according to

$$\frac{Q_{\varphi }}{I_{\star }}=\frac{\sigma }{\kappa }\frac{⟨f_{\phi }(i,g,p_{\mathrm{m}\mathrm{a}\mathrm{x}})⟩}{⟨f_{\mathrm{I}}(g=0)⟩}\frac{L_{\mathrm{w}\mathrm{a}\mathrm{r}\mathrm{m}}}{L_{\star }}.$$(3)

We obtain for i = 60°, g = 0.6 and p_max = 0.25 the disk phase function ⟨f_φ(60°, 0.6, 0.25)⟩/⟨f_I(g = 0)⟩ = 0.085, and with σ/k = ω_d/(1 − ω_d) = 1.0 the intrinsic, disk integrated polarization signal of Q_ϕ/I_* = 0.085 · L_warm/L_★ = 2.8 × 10⁻⁶ for the warm dust in ε Eri.

For the comparison with our observations we can now calculate the Q_ϕ signal with the dust parameters given above and investigate the expected surface brightness for the polarized signal SB_p for different disk inclinations and radial distributions of the dust. Figure 4 shows four cases with a mean ring radius r_c = 1.25″, one with a small width of Δr_narrow = 0.1 · r_c and i = 60°, and three broader disk rings Δr_broad = 0.4 · r_c, with i = 30°, i = 60° and i = 75° to investigate the inclination dependence. All four models produce the same amount of infrared excess L_warm/L_★ = 3.3 × 10⁻⁵ and they are convolved with the mean PSF of the ε Eri data for a prediction of the expected observational signal Q_ϕ/I_★ as given in Table 1 (Col. 5).

Most important for dust detection is the resulting surface brightness of the polarized intensity SB_p which depends strongly on the disk geometry. It is clearly visible in Fig. 4, that a much higher surface brightness SB_p is obtained for inclined disks and for the narrow dust ring when compared to the wide ring. In our simple models with constant dust emissivities ε(r) = ε₀, the intrinsic surface brightness is anti-correlated with the disk widths SB_p ∝ 1/Δr. After convolution the ratio of the peak S B_p between the narrow and wide ring models is about a factor of three, as the convolution degrades more the peak surface brightness flux of narrow structures.

The three models for the wide disks illustrate that SB_p is lower for low i because only a small amount of light is scattered towards the observer by dust with strong forward scattering. For high i the intensity surface brightness SB_I of the front side increases strongly because of the strong forward scattering. The polarization SB_p has a dip on the disk front for i = 75° because small scattering angles θ < 30° produce only a weak polarization p(θ). The locations with the maximum SB_p-signal shifts for higher i towards the apparent major axis of the projected disk, where the scattering angle is close to θ ≈ 90° and therefore p(θ) ≈ p_max.

The models provide expected maximum surface brightness contrasts ΔSB_p which is ΔSB_p = SB_p − m_* for the brightest disk section SB_p measured relative to the central star m_*. For the narrow disk model (i = 60°) the intrinsic contrast is ΔSB_p = 13.3 mag arcsec⁻², or ≈13.9 mag arcsec⁻² if we also consider the signal degradation by the PSF convolution. Compared to the stellar PSF peak SB(0), the brightest disk region is about SB_p − SB(0) ≈ 20.3 mag arcsec⁻² fainter because there is SB(0) − m_* ≈ −6.4 mag arcsec⁻² between the peak and the total brightness of the central star. This disk polarization to PSF peak contrast is about seven times higher (≈7 × 10⁻⁹) than the point source contrast calculated for the RV planet in Sect. 2.1. In addition, the two brightest regions of the disk on opposite sides of the ring are roughly 0.125″ wide and 0.4″ long and cover each a surface area of about 50 × 10⁻³ arcsec², while a point source has only a size of about 1 × 10⁻³ arcsec². Thus, the signal-to-noise ratio for the predicted disk polarization signal is about a factor of ~50 higher than for the RV planet, if the data would only be limited by photon noise.

The situation is less favourable for wider disks because the peak contrasts are only ΔSB_p ≈ 15.6, 15.1, and 14.4 mag arcsec⁻² (Table 1) for i = 30°, 60° and 75°, respectively. The signal is distributed over a larger area and might therefore become visible with strong pixel binning.

The intensity signal of the scattered radiation has for inclined disks a strong maximum on the disk front side because of the strong forward scattering (Fig. 3). The highest surface brightness contrast for the narrow disk with i = 60° is about ΔSB_I = 10.5 mag arcsec⁻² or a factor of 22 brighter than for the peak polarization signals. Despite this, we expect for SPHERE/ZIMPOL observations a higher disk detection sensitivity with polarimetry because the polarized disk signal can be distinguished from the strong steller intensity halo. Nonetheless, a bright disk detected with imaging polarimetry could also produce a detectable intensity signal in the data. Figure 3 shows also the expected surface brightness distribution for the thermal radiation of the dust in the mid-infrared. All these model results depend strongly on the adopted disk parameters. Therefore, a search should also consider signals which could be a factor of a few higher or lower than the predictions for the selected disk models shown in Fig. 4.

Fig. 3 Relative flux distribution for a narrow dust ring model in ε Eri with r = 4 au based on the warm component in the infrared excess. Top: thermal emission. Middle: expected scattered intensity 1disk for g = 0.6. Bottom: polarized intensity Qϕ for g = 0.6. The model images are convolved with the PSF of SPHERE/ZIMPOL, normalized to their flux maximum and aligned with the sky coordinate axes.

Fig. 4 Predicted polarized flux Qϕ [ct/(s · px)] produced by scattering in a narrow and a broad dust ring model for ε Eri with an inclination of i = 60° (top row), and for broad rings with i = 30° and i = 75° (bottom row). All models would produce the measured infrared excess of Lwarm/L★ = 3.3 × 10−5.

3 Observations and data analysis

3.1 SPHERE/ZIMPOL instrument

The planet ε Eri b was searched with direct imaging using the SPHERE instrument (Beuzit et al. 2019) at the Nasmyth focus of the VLT unit telescope UT3 of the European Southern Observatory (ESO). This instrument consists of an extreme adaptive optics (AO) system, an image derotator, stellar coronagraphs and three focal plane instruments (Fusco et al. 2006, 2014; Petit et al. 2014; Sauvage et al. 2014) including the Zurich Imaging Polarimeter (ZIMPOL) used for this programme. ZIMPOL works in the visual spectral regime 500–900 nm and is tuned for the search of reflecting planet and circumstellar disks using fast-modulation imaging polarimetry with a modulation frequency of about 1 kHz to ‘freeze’ the speckle variations (Schmid et al. 2018). The SPHERE AO system achieves under good observing conditions regularly a Strehl ratio of about 40% in the I-band (Fusco et al. 2015) and a resolution of about 25 mas full width at half maximum (FWHM). ZIMPOL has a detector field of view of 3.6″ × 3.6″ and a pixel scale of 3.6 mas × 3.6 mas and is equipped with different coronagraphs, filters, and instrument calibration components. For polarimetry, the telescope and instrument polarization are compensated and calibrated with rotatable half-wave plates and further calibration components (Bazzon et al. 2012). ZIMPOL has two arms, each of them with one camera and its own filter wheel, and data are taken simultaneously in both arms. The two CCD detectors are operated in frame transfer mode and provide a high gain, fast read-out mode with small detector overheads for high flux applications.

The deep planet search with SPHERE/ZIMPOL is based on the high-contrast imaging, which provides at a separation of about 1″ a raw contrast at the level of 10⁻⁴ using AO and coronagraphy. Combining this with polarimetric differential imaging (PDI) and angular differential imaging (ADI) gives an additional polarimetric contrast improvement of about 10⁻⁴. This provides a total contrast at the level of about C_P ≈ 10⁻⁸ and this can be improved further by reducing the photon noise with sufficiently long integrations for bright targets (Schmid et al. 2006a; Thalmann et al. 2008; Hunziker et al. 2020).

3.2 Observations

The observations of ε Eridani were taken in visitor mode with an observing strategy similar to the observations of Hunziker et al. (2020). Three runs were executed between Oct. 2019 and Jan. 2020 and one run in Nov. 2020. This planning considered the planet motion from run to run, and the combination of the data with a Keplerian motion prediction. The observation log and basic information on exposure times and observing conditions are listed in Table 2.

The ZIMPOL instrument setup was optimized for the highest possible throughput, using the Very Broad Band (VBB) filter (λ_c = 735.4 nm, Δλ = 290.5 nm) for both camera arms, polarimetry in fast modulation mode and a classical Lyot coronagraph (V_CLC_MT_WF) with a spot radius of ρ = 77.5 mas. The mask has a transmission of around 0.1% so that, under good conditions the stellar spot is visible behind the Lyot mask and can be used for centring. Approximately once per hour we were offsetting the star from the coronagraph and added the ND2 neutral density filter (reduces flux by a factor of about 100) to measure the flux throughput and the PSF shape for an improved beam-shift correction (see Appendix A.3).

Observations are taken in polarimetric cycles which consist of a sequence of Q⁺, Q⁻, U⁺ and U⁻ images taken with four different half wave plate (HWP) orientation each with nDIT subintegrations. The Q⁺ and Q⁻ data provide the total intensity for the Q-measurement Iq = I₀ + I₉₀ and Stokes Q = I₀ − I₉₀, which is already corrected for the instrument polarization thanks to the switch between the Q⁺ and Q⁻ measurement. This does not correct for the telescope polarization introduced before the HWP switch. The equivalent parameters are obtained for Stokes U = I₄₅ − I₁₃₅ and I_U = I₄₅ + I₁₃₅ with I_Q = I_U = I. The detector integration time DIT = 3 s was initially used for the first two runs, but then extended to DIT = 5 s for a higher count level with typically >100 counts, or photo electron numbers of N_e > 1000 e⁻(gain factor of 10.5 e⁻/ct) for ρ ≳ 1.2″ (see intensity PSF in Fig. 5). This level is required to achieve photon noise limited observations, with (N_e)^1/2 larger than the read-out noise level of N_ron ≈ 20 e⁻ (2 counts).

The ZIMPOL P1 derotator mode was used which is optimized for polarimetry because the derotator and all other components after the rotating and switching HWP2 with their corresponding instrument polarization effects are constant during the night. Therefore the sky rotates in the image and very importantly this allows to use ADI together with PDI to correct better for quasi-static speckles and other instrumental effects of SPHERE and ZIMPOL. However, in P1 mode also the telescope pupil rotates, most notably the M2 spider pattern, but with a different rotation law than the sky image.

3.3 Data reduction

The data were mainly reduced with the IDL-based sz (SPHEREZIMPOL) software developed at the ETH Zurich. Basic steps include signal extraction for the two simultaneously measured polarization modes, bad pixel cleaning, bias subtraction, flat-field correction, calibration of the polarimetric modulation efficiency, and the polarimetric combination of the four frames of each cycle. Additionally the frame transfer smearing was corrected in the intensity frames by subtracting the average row level multiplied by 56 ms/DIT to account for the illumination during the frame transfer for the detector mode used for fast modulation polarimetry. Important for our analysis are the corrections for the differential polarimetric beam shifts (Schmid et al. 2018; Hunziker et al. 2020), which can be tricky to define from coronagraphic data. Therefore we also use the PSF images as described in Appendix A.3. It is important for polarimetry to consider the telescope polarization which adds a fractional polarization p_tel with a position angle δ_tel to the signal. This is rotating with the parallactic angle of the telescope θ_para so that the uncorrected polarization of ε Eridani rotates along a circle in the Q/I-U/I-plane. The correction for the telescope polarization and the second order radial dependence thereof are described in Appendix A.2, while particular detector corrections are discussed in Appendices A.1 and A.4.

For the search of scattered light from a planet and from an extended cloud of dust, we use the azimuthal Stokes parameters Q_ϕ and U_ϕ with respect to the central star. The following formulas relate the Stokes Q and U to the azimuthal polarization Q_ϕ = −Q cos(2ϕ) − U sin(2ϕ) and U_ϕ = −Q sin(2ϕ) + U cos(2ϕ), where ϕ is the position angle with respect to the central star measured from north over east (Schmid et al. 2006b; Monnier et al. 2019). The scattered light from a planet or a pole-on disk is mostly azimuthally polarized and therefore should be visible in the Q_ϕ image only, while U_ϕ is zero.

In Fig. 5 we show one complete polarimetric cycle from night 9 after all the data reduction steps mentioned above. The plot is splitted to emphasize the large difference in residual signal between the total intensity I in the top half and Stokes Q_ϕ for the lower half. In the intensity half, one can clearly see the bright speckle ring at around 0.45″ up to which the AO system corrects well the atmospheric seeing. Outside of 0.5″ the residual intensity halo of the star decreases steadily with separation. Additional features are the four cross shaped stripes from the secondary mirror mount (telescope spiders) which are partially visible, dominant PSF features left and right of the star in the speckle ring, and small black astrometric dots from the coronagraph mask, for example at (y = +1″, x = −1″, 0″, +1″). In Schmid et al. (2018) full plots of many PSFs for SPHERE/ZIMPOL are shown and described. The differential polarization signal Q_ϕ is much smoother and with significantly lower count numbers than the intensity image (note colour scales in Fig. 5). The exposure time of this polarimetric cycle is NDIT × DIT × four half wave plate positions = 12 · 5 s · 4 = 240 s.

Fig. 5 One polarimetric cycle texp = 240 s corrected for the telescope polarization: the top half shows the intensity I and bottom half the polarized intensity Qϕ with corresponding colour scales in ct/(px· DIT). The counts inside 0.9″ are scaled down for a better visibility by a factor of six for the intensity and a factor of three for Qϕ.

Fig. 6 All data of night 9 (texp = 9360 s), derotated, averaged, and corrected for telescope polarization: the top half shows I and the bottom half Qϕ with corresponding colour scales in ct/(px · DIT) as in Fig. 5. Artificial point sources are inserted at different separations and contrasts CI and CP are changed with position angle as indicated. Counts inside 0.9″ are scaled down for better visibility with a factor of six for I and a factor three for Qϕ.

4 Planet search

4.1 Pushing the limits

Additional steps are necessary for the search of a reflecting planet because the expected signal is much weaker than the stellar halo in the coronagraphic intensity image or the noise in the differential Q_ϕ image. The aim is essentially to reduce the relative noise in the halo by averaging a large amount of data such as in Fig. 6 where all observations of night 9 are rotationally aligned, so that north is up and east to the left, and combined to enhance the S/N for the detection of a faint companion. We have injected fake planets at different separations and azimuthal locations where the points for a given angle have all the same contrast value as indicated. A source with a given contrast value is easier to detect at larger separations and much easier for a polarized contrast value C_P in the Q_ϕ image than the corresponding intensity contrast C_I in the intensity image. The PSF peak flux of the star is roughly 6.5 × 10⁶ ct/(px · DIT) and therefore a planet with a contrast of 10⁻⁶ has 6.5 ct/(px · DIT).

The averaging reduces the noise for Q_ϕ to a level of about +0.1 ct/(px · DIT) depending on separation and this would allow a detection of a point source with a contrast of about C_P ≈ 10⁻⁷ outside a separation >0.9″ as demonstrated by the artificial point sources inserted in the image. The rotational alignment smooths the stellar halo in the intensity image and averages out localised features fixed to the instrument such as the telescope spiders or the dark points from the coronagraphic mask. However, higher contrast limits require the subtraction of the strong PSF halo in the intensity data and a few improvements for the polarimetric data as described below.

The data from ε Eri, such as those shown in Figs. 5 and 6, allow rough estimates of the photon noise in the images compared to a planet signal with a contrast of C_I = 4 × 10⁻⁹. Such a planet would produce an intensity signal with a PSF peak flux of 0.004 ct/(px · s) because the corresponding ε Eri peak flux is about 10⁶ ct/(px · s). This competes at ρ ≈ 1″ with a count intensity for the stellar halo of about I_ct ≈ 40 ct/(px · s) or the total number of photons collected during 38.5 h of n_γ ≈ 6 × 10⁷ px⁻¹. This considers the detector gain of 10.5 e⁻/ct. The corresponding relative photon noise limit (n_γ)^−1/2 = 1.3 × 10⁻⁴ can be expressed as statistical noise limit per pixel I_ct · (n_γ)^−1/2 ≈ 0.005 ct/(px · s) for the entire data set around 1″.

This is comparable to the expected PSF peak of the planet with C_I = 4 × 10⁻⁹ or a S/N ≈ 1 for a pixel near the PSF peak of the planet. Because the PSF has a FWHM of about 6 pixels the S/N for the whole planet PSF would be roughly at a level ≈5 with respect to the photon noise of the stellar halo.

The estimated polarimetric signal of the planet is lower, about C_P ≈ 0.25 C_I = 1 × 10⁻⁹ and one needs to take Q and U measurements, if the position angle of the polarization is not known. Therefore, twice the measuring time is required to collect 2 n_γ to reach a required relative photon noise limit of (n_γ)^−1/2 for the signal in the Q_ϕ image. However, only in polarimeteric imaging the detection limit is close to the photon noise, while in intensity imaging the detection limits are clearly above the photon noise because of the strong speckle noise (Sect. 4.6).

4.2 Post-processing

To push the noise levels further we combine in our postprocessing ADI with fitting and subtraction of the PSF halo in the intensity imaging or of the large scale pattern in the residual differential polarization for polarimetric imaging. We also apply additional steps to suppress special instrument features in the data.

The subtraction of the structure of the stellar PSF halo is rather simple for the calibrated polarimetric imaging data of ε Eri because the residual signal of a polarimetric cycle is very weak (see Fig. 5). For ρ > 0.6″ the systematic structures are smaller than the photon noise except for residuals at the position of the eight astrometric spots of the coronagraph mask and some along the telescope spider. Residuals from the telescope spider remain because they rotate in our data with another rotation law than the sky field and this introduces small systematic alignment errors in the combination of the polarimetric data. The features from the coronagraphic spots in the polarized intensity are introduced by the beam shift calibration which corrects the sky image for the differential polarimetric shifts introduced by the inclined mirrors of the VLT and the SPHERE instrument. This correction produces artifical polarimetric beam shifts for all image features introduced by components located after the inclined mirrors, such as the astrometric spots of the used coronagraph or dust on optical components of the science cameras. The position of these residual spot features vary during the night after centring the star because of drifts in the alignment between the star and the focal plane coronagraph. Therefore a simple subtraction procedure gives unsatisfactory results in the residual images. Consequently, we mask these localized features in all Q_ϕ images with pixels having ‘not a number’ (nan) values. This masking leads to a loss of around 6% of the photons in return for a much smoother residual image. We then used for the subtraction of the residual halo structure in all Q_ϕ-images of a night the median of the masked, non-derotated data of that night.

The situation is more delicate for the intensity images, because the stellar intensity halo is strong and quite variable due to AO performance variations with strong short-lived speckles and the overall PSF halo structure changes. For this reason, the stellar halo is fitted for each image using a principal component analysis (PCA) (Amara & Quanz 2012) and then these fits are subtracted from the data. This procedure was carefully investigated to avoid the introduction of spurious point-like features or possible self-subtraction of real point sources in the resulting data residuals, which are used for the search of a planet. Best results are obtained using about 20 principal components in the PCA analysis. The astrometric spots from the coronagraph (one such spot is visible in Fig. 5, 1″ above the centre), and the telescope spiders were not masked because the use of nan – values introduces strong spurious effects in the averaging of data with variable flux levels.

Finally, we compared different combination methods, mean, median, or noise-weighted mean (Bottom et al. 2017), for the derotated, differential polarimetric and intensity images of each night. Best results for both types of data are obtained with noise-weighted means (in case of the polarized intensity with masked spiders and astrometric spots).

The resulting data for night 9 after all these post-processing steps are shown in Fig. 7. The PSF halo subtraction for the intensity image yields as final data sets difference images with a mean value of zero and deviation at the level of a few counts for separation >0.9″ and deviations of the order 10 counts inside <0.9″. This allows to spot artificial point sources which are ten to 30 times fainter than what could be seen in the averaged imaging data of Fig. 6. Dominant noise sources in the halo subtracted intensity image are the speckle noise at small ρ and the residuals from the diffraction pattern of the telescope spider for larger ρ. The spider effect is a result of the ZIMPOL P1 mode because ZIMPOL offers no polarimetry with pupil stabilization. The spider effect is small, less than 1% of the halo intensity but it might still be beneficial to correct for this in future investigations.

The post-processing of the Q_ϕ image removed most instrument features and shows for night 9 a very clean pattern of Gaussian noise on top of which the artificial sources with a contrast of C_P = 10^−6.5 and C_P = 10⁻⁷ are clearly visible. This illustrates that the differential polarization data reach a contrast very close to the photon noise limit for the search of point sources.

Fig. 7 All data of night 9 (texp = 9360 s), derotated, PSF subtracted, masked regions for the Qϕ, averaged and corrected for telescope polarization: the top half shows I and the bottom half Qϕ, with corresponding colour scales in ct/(px · DIT) as in Fig. 5. Artificial point sources are inserted at different separations and contrasts CI and CP are changed with position angle as indicated. Counts inside 0.9″ are scaled down for better visibility with a factor of six for I and a factor three for Qϕ.

4.3 Contrast curves and detection maps

To characterize the sensitivity of our search, we inject fake planets in the images before the post-processing and then evaluate how well we can retrieve the planet. The template for the fake planets for a given night is the median PSF from unsaturated images of ε Eri taken with the ND2 filter. The counts are scaled to the flux of the coronagraphic images used for the planet search considering the wavelength dependent telescope and instrument transmission for the broad VBB filter with and without ND2, the atmospheric transition, and the ε Eri spectrum. The derived scale factors are 174 for nights 1 and 2 when the data were taken with 3 s DITs, and 290 for the remaining nights taken with 5 s DITs, respectively.

We then injected many such template PSFs multiplied with a small contrast factor (e.g. 10⁻⁶) as fake planets in the science image. The fake planets were distributed in a spiral pattern with sufficient separation to avoid cross talks in the signal extraction. We then changed the contrast factor until the planets could be retrieved with a 5 σ_{$\mathcal{N}$} Gaussian significance detection. This procedure was repeated six times with the planet spiral pattern rotated each time by 60° so that for each radial separation the planet was injected at six different angle positions. Finally the mean of the six measurements was taken. This whole process was repeated with planets at radii between the previous radii to get a denser radial sampling. For the planet apertures a radius of 5 pixel (or 18 mas) was used as this is the optimal size to extract as much point source signal as possible while keeping at the same time the background noise low. Around each aperture an annulus of 4 px width and 2 px separation to the aperture was used to subtract the local background. This background subtraction does not improve the contrast, but reduces the standard deviation between the six injected planets for a given radius.

We use for the detection metric the false positive fraction (FPF) as described in Mawet et al. (2014). As explained in Bonse et al. (2023), Gomez Gonzalez et al. (2017) and Christiaens et al. (2023) the S/N for the student’s t-distribution should not be confused with the Gaussian sigma significance. In short, we first calculate the S/N as defined in Mawet et al. (2014) (based on two sample t-test): $T=(P-B)/\left(s_{\mathrm{B}}\sqrt{1+1/n}\right)$, where P is the planet signal measured in the aperture, B the mean and s_B the standard deviation for the background apertures, and n the number of background apertures. The same T value corresponds to different FPF for different separations λ/D. T follows a student t-distribution with n − 1 degrees of freedom and one can calculate for each radius the $\mathrm{F}\mathrm{P}\mathrm{F}={\int }_T^{\mathrm{\infty }}p(T=t\mid H_0)\mathrm{d}x$ (Bonse et al. 2023) for a found value of T. For better readability we express the FPF in terms of the quantiles of the standard normal distribution, so that the 5 σ_{$\mathcal{N}$} Gaussian significance corresponds to a FPF of 2.87 × 10⁻⁷ and 3 σ_{$\mathcal{N}$} to a FPF of 1.35 × 10⁻³. Because n increases with separation, a measured S/N value of T = 5 corresponds to a 4.76 σ_{$\mathcal{N}$} at a separation of 500 mas (4.87 σ_{$\mathcal{N}$} at 1000 mas and 4.92 σ_{$\mathcal{N}$} at 1500 mas). It was checked that the behaviour of the noise was close to Gaussian to fulfil the statistical assumptions, as expected for planet signals located at separations of ρ > 6(λ/D).

We derived with this method for ε Eri azimuthally averaged radial contrast curves shown in Fig. 8 for 5 σ_{$\mathcal{N}$} Gaussian significance for the intensity and the polarization for each individual night and for each epoch. The twelve nights are spread over four epochs of a few consecutive nights as indicated in Table 2. The PSF changes significantly from one night to another night and therefore also the contrast curves. Because the expected motion of a planet around ε Eri is less than 0.7 pixels per night (Sec. 2.1) we can combine the nightly results from one epoch to a time weighted ‘mean epoch’ data set.

We also produce detection maps for the Q_ϕ-data to search for a significant signal of a polarized point source in the entire field of view (Fig. B.1). For this we use for an individual night the averaged, residual Q_ϕ frame and treat each pixel as central pixel of an aperture and calculate the Gaussian sigma detection values. Detection maps for the epochs are based on time weighted Q_ϕ averages of the individual nights. A source with a significant azimuthal polarization Q_ϕ should then show up in these maps as bright spot with a significance of >5 σ_{$\mathcal{N}$}.

Fig. 8 Sensitivity expressed as contrast curves with 5 σ$\mathcal{N}$ Gaussian significance. Top: individual nights for the intensity CI = Ip/I★ (upper curves), the polarized light CP = pp · Ip/I★ (lower curves). Bottom: the four epochs in colour in comparison to the individual nights in grey.

4.4 Results: Individual nights

The sensitivity of our survey is illustrated with 5 σ_{$\mathcal{N}$} contrast curves in Fig. 8 for all the individual nights for the intensity (upper curves) and the polarized light Q_ϕ (lower curves) as a function of separation ρ from the star.

Polarized intensity. For ρ = 1″ the mean polarized contrast of the twelve nights is (5.57 ± 1.85) × 10⁻⁸. Between the best night 10 ((2.90 ± 0.37) × 10⁻⁸) and the least sensitive night 8 ((9.0 ± 1.42) × 10⁻⁸) is roughly a factor of three independent of the separation. The uncertainties reflect the standard deviation of the contrast values derived for the six fake planets with the same separation but inserted at different position angles. At ρ = 0.6″ the mean contrast is (11.7 ± 3.81) × 10⁻⁸ and at 1.6″ it is (2.79 ± 1.04) × 10⁻⁸. The contrast is less good at small ρ because of strong, short lived speckles, larger halo flux, and less field rotation in absolute pixel values. For ρ > 1″ no contrast curve improvement could be obtained after PDI and ADI with different post-processing methods such as median or PCA halo subtraction probably because PDI and ADI already achieved the photon noise limit.

Comparing the nights 5 and 6 exemplifies well the effect of the observing conditions on the achieved limits. The integration time was similar, but the conditions and the PSF were much better in night 6 with better seeing and longer coherence time τ₀ = 6.1 ms instead of only 2.5 ms in night 5 (Table 2). The achieved contrast limits are approximately 2.5 times better for night 6 than for night 5, also because 20 min of integration time was unusable in night 5 because the AO system was not stable.

The Gaussian significance detection maps for the polarized light Q_ϕ of the individual nights are shown in the top panel of Fig. B.1. Out of the twelve images, each consisting of 723 736 pixels, there are 37 pixels with a significance larger than 5 and these 37 pixels belong to 10 points of neighbouring pixels (see Table 3). As the Gaussian significance corresponds to a FPF of 2.87 × 10⁻⁷ one would expect roughly 2.5 pixels with a σ_{$\mathcal{N}$} > 5 in the twelve images assuming perfect Gaussian noise.

One should note that most points with more than 5 σ_{$\mathcal{N}$} signal have a contrast significantly brighter than our expectation of about 10⁻⁹ for ε Eri b (see Sec. 2.1). Furthermore some points are at a separation ρ larger than expected from the RV orbit. In particular, the predicted separation for the fourth epoch should be close to 1.16″, making detections at much smaller or larger separations unlikely. The noise properties of the U_ϕ residual images, in which we expect no signal from a planet, look indistinguishable from the Q_ϕ images.

Total intensity. We get for the intensity contrast curves for the 12 nights in Fig. 8 a 5 σ_{$\mathcal{N}$} mean contrast of (6.15 ± 2.90) × 10⁻⁶ at ρ = 0.6″, (1.58 ± 0.83) × 10⁻⁶ at ρ = 1″, and (2.55 ± 1.15) × 10⁻⁷ at ρ = 1.6″. In the intensity residual images there is more systematic noise left from short-lived speckles and the artefacts related to the telescope spiders. We expect that some contrast improvement might be possible with future, more sophisticated post processing methods.

4.5 Results: Epochs

Polarized intensity. The 5 σ_{$\mathcal{N}$} Gaussian significance contrast curves for the epochs consisting of 2–4 consecutive nights are illustrated in Fig. 8. The mean contrast for the Q_ϕ polarization for the four epochs is (7.04 ± 2.14) x 10⁻⁸ for ρ = 0.6″, (3.29 ± 1.01) × 10⁻⁸ for 1″, and (1.60 ± 0.54) × 10⁻⁸ for 1.6″. Epoch 4 is most interesting because the planet separation is at a maximum, the scattering angle is close to 90° and ideal for a strong Q_ϕ signal from the planet, and epoch 4 is the longest with the best contrast: (5.02 ± 0.76) x 10⁻⁸ for 0.6″, (2.11 ± 0.43) x 10⁻⁸ for 1″, and (1.10 ± 0.24) x 10⁻⁸ for 1.6″. The three nights of epoch 4 have an individual mean contrast of approximately 4 x 10⁻⁸ for ρ = 1″ and each night is about 5 hours long. The combined epoch sensitivity is approximately $\sqrt{3}$ times higher than for one night or an improvement similar to $\sqrt{t_{\text{exp }}}$. Most notably, the achieved 5 σ_{$\mathcal{N}$} contrast is (1.22 ± 0.32) × 10⁻⁸ for the expected maximum planet separation of 1.16″ occurring around epoch 4.

In lower panel of Fig. B.1 we show the Gaussian significance epoch maps. For perfect Gaussian noise we would expect 0.8 pixels with a significance greater than 5 and we find 4 pixels with a σ_{$\mathcal{N}$} > 5 attributed to two locations (see Table 3): In epoch 3 with σ_{$\mathcal{N}$} = 5.04 at ρ = 472 mas at an angle (north over east) of 101°. This would correspond to a contrast of (2.21 ± 0.46) × 10⁻⁷ which is at least one order of magnitude higher than the expected planet signal. The other spot is in epoch 4 and consists of three neighbouring pixels with a σ_{$\mathcal{N}$} > 5 with the central pixel at ρ = 686 mas, 337° and significance of 5.1, which corresponds to a contrast of (4.32 ± 0.66) × 10⁻⁸. The separation is smaller than expected at that time and the contrast is a rather high value. Another point in epoch 4 which might be worth mentioning is at ρ = 994 mas, 30° with σ_{$\mathcal{N}$} = 4.33 and contrast 1.96 ± 0.46 × 10⁻⁸.

Total intensity. For the search of the intensity signal of a companion the mean 5 σ_{$\mathcal{N}$} contrast for the four epochs is (3.06 ± 1.36) × 10⁻⁶ for 0.6″, (6.05 ± 1.73) × 10⁻⁷ for 1″ and (1.21 ± 0.46) × 10⁻⁷ for 1.6″. For most separations epoch 1 is the most sensitive although for example epoch 4 contains more than twice more observing time. Probably collecting more individual images (DIT = 3 s instead of DIT = 5 s) offers an advantage in post-processing and PSF halo subtraction. This consideration does not take into account the read noise at larger separation as the fake planet injection uses a high signal PSF template multiplied by the small contrast number. For the PSF subtraction technique in the post-processing, the used PCA with 20 principal components gives approximately a factor two contrast improvement when compared to a median PSF subtraction. For small separations this factor is slightly larger and for larger separations a bit smaller. Going through the full parameter space and combining the best PCA contrast curves with different number of components (see e.g. contrast curve documentation of applefy package Bonse et al. 2023), using different settings of aperture sizes depending on the night or using more advanced PSF subtraction techniques, for example, (see e.g. Cantalloube et al. 2020; Gebhard et al. 2022) could lead to an additional improvement in contrast.

4.6 Advantage of using PDI

The ε Eri data were taken with angular differential imaging (ADI) and polarimetric differential imaging (PDI) with the use of the ZIMPOL P1 mode. This offers an ideal opportunity to compare the performances for the high-contrast searches of the intensity signal and the polarization signal of point sources because the analysis can be based on data taken with the same instrument and under the same observing conditions.

Thus, we compare observing and data analysis methods for the total intensity planet search using ADI, plus PSF fitting and subtraction with a PCA method (Amara & Quanz 2012) with the polarimetric planet search using PDI, ADI and residual pattern subtraction. This should provide a very reliable assessment because we are using state of the art procedures for the observations, the data reduction and post-processing. However, it should be noted, that the Strehl ratio provided by the SPHERE AO system is about 40% for the short wavelength range of ZIMPOL, which is significantly lower than for the near-infrared (Fusco et al. 2015).

Polarimetry is a very efficient high-contrast technique because it provides the differential signal of the opposite polarization modes simultaneously and this minimizes very strongly the temporal variability effects of the speckle noise (see e.g. Schmid 2022).

The advantage of polarimetry is clearly visible from the much better polarized contrast curves C_P when compared to the intensity contrast curves C_I in Fig. 8. The ratio C_P/C_I based on the curves for the mean contrast limits derived from the full data set as function of separation is shown in Fig. 9.

A planet detection is easier with polarimetic imaging for planets with a fractional polarization p > C_P/C_I, and easier with intensity imaging for planets with p < C_P/C_I. Especially for small ρ, searches in polarized light offers a strong advantage because of the very efficient suppression of unpolarized speckles.

Of course, the polarization signal C_P from a planet is strictly smaller than the intensity signal C_I, typically by a factor of 3 to 20 for a scattering angle of about α ≈ 90°. The measured fractional polarization p_p of solar system gas planets in the R-band are for the Rayleigh scattering planets Uranus and Neptune about p_p ≈ 20% (Schmid et al. 2006a; Buenzli & Schmid 2009). Reflections by clouds produce less p_p while atmospheric haze can produce a very high polarization p_p of up to 50%. Therefore the integrated polarization of Jupiter is about 10% (Smith & Tomasko 1984) for equatorial sight lines and about p_p ≈ 15% for polar sight lines when the reflection from the polar haze is well visible. For Saturn p_p is only about 5% (without disk) because of the predominant atmospheric clouds Tomasko & Doose (1984).

For the selection of the best observing strategy for future planet searches with SPHERE/ZIMPOL one should carefully consider the ratio curve C_P/C_I in Fig. 9 together with the expected fractional polarization of the planet. This curve is based on a large and representative data set and it varies only slightly between observations taken under good or bad atmospheric conditions. This can be inferred from the C_I and polarization C_P contrast curves in Fig. 8, which are for given nights both going up or down in step with the observing conditions. The C_P/ C_I curve derived for SPHERE/ZIMPOL gives a useful benchmark for the design of future instruments, but one should also consider that the overall performance depends on many instrument parameters for the AO-system and the differential imaging concept.

Fig. 9 Ratio between polarization CP and intensity CI contrast limits for the full data set, illustrating the advantage of PDI for the speckle suppression. Planets with a fractional polarization p > CP/CI should be easier to detect in the polarimetric data.

4.7 Results of K-Stacker orbital search

Combining consecutive individual nights into a set of four different epochs is straightforward because the expected orbital motion from day to day is less than 3 mas or one ZIMPOL detector pixel. Improving the detection limits by the combination of the data from different epochs requires much more care. Adopting for ε Eri an orbital semi-major axis of a ~ 3.5 au around a 0.82 M_⊙ star gives an orbital period of about 7 yr, or an astrometric motion of 3 au yr⁻¹, or 1″ yr⁻¹ for the distance of 3.2 pc. This converts to about 80 mas per month for a face-on, circular orbit, which is more than 3 times larger than the width (FWHM) of the PSF of SPHERE/ZIMPOL. Thus, the epochs must be properly combined considering the orbital motion of any putative planet. To do so, we used the K-Stacker algorithm (Nowak et al. 2018; Le Coroller et al. 2020), which searches for potential companions in series of images along a grid of pre-determined orbital parameters. For this we use the Q_ϕ-maps of all nights, such as the one shown for night 9 in the lower half of Fig. 7.

4.7.1 Potential solutions

The overall process can be summarized in a few steps according to the detailed description of the algorithm in the reference publications (Nowak et al. 2018; Le Coroller et al. 2020):

1. The algorithm first calculates the planet position for each night for a set of orbital parameters p.

2. For each night t, it extracts a ‘signal’ value s_t(p) which is the integrated photometry in a circle of 6 px diameter, which corresponds to the size of the instrumental PSF, as well as a ‘noise’ value n_t(p) calculated from the distribution of signal values extracted along a circle whose radius corresponds to the separation at night t for orbit p.

3. The algorithm calculates the total estimated S/N for all points of an orbit p using: $(\mathrm{S}/\mathrm{N})(p)={\sum }_{\mathrm{t}}{s}_{\mathrm{t}}(p)/\sqrt{{\sum }_{\mathrm{t}}{n}_{\mathrm{t}}(p{)}^2}$. It then ranks all the orbits of the grid by order of decreasing S/N.

4. A subset of the best orbits (in our case, the best 70 orbits) are further optimized using a gradient-descent algorithm to allow for parameter values between the initial grid-points.

5. The calculation ends with a report of all the values calculated along the grid of orbital parameters, and the optimized results for the best 70 orbits.

To setup the grid of orbital parameters, we followed Nowak et al. (2018) and Le Coroller et al. (2020), and first determined the typical width of a maximum in the (S/N)(p) function along the different parameters. We then define the step sizes for each parameter as typically one fifth of the corresponding (S/N)(p) peak width, to ensure that K-Stacker would not miss any potential S/N maximum. The explored range of orbital parameters in our grid was defined based on previous studies of the planet ε Eri b, mainly on Llop-Sayson et al. (2021). The used parameter grid is given in Table 4. It restricts orbital periods according to (P[yr])² = (a[au])³/M_star[M_⊙] to the range between about P ≈ 5.5 yr and 10.9 yr and the orbits’ eccentricity to e ≤ 0.6. The grid search does not constrain the orbital inclination i and the orientation of line of nodes with respect to sky plane Ω, nor the argument of the periapsis ω. Thus, the search allows for prograde and retrograde planetary orbits.

Among the 6 × 10¹¹ planetary orbits explored by the algorithm, the best solution found reached a S/N(p) of up to 7.4. The re-optimization of the best 70 solutions with a gradient-descent algorithm leads to typical improvements of about 0.4 in the S/N, and it yields a distribution of best results along 3 potential orbits, represented in Fig. 10.

Solutions 2 and 3 obtained by K-Stacker (see Fig. 10) are actually very similar in terms of position for the individual epochs, with a difference of at most 50 mas for the last epoch, and about 15 mas for the other epochs. This suggests that the algorithm has actually caught some feature within the images. The emergence of two solutions arise from degeneracies in the orbital parameters, a consequence of the limited number of independent epochs available for our data. Interestingly, solution 2 shows orbital parameters which are largely compatible with the solution presented by Llop-Sayson et al. (2021), as shown in Table 5. It should be noted, though, that despite their use of astrometric data in the analysis, their orbital parameters are poorly constraining planet positions, with large error bars on i and Ω, which increases the chance of having “compatible solutions”. We note, however, that solution 1 found by K-Stacker corresponds to a completely different set of orbital parameters, and to different positions at the ZIMPOL epochs. This demonstrates that K-Stacker converges towards multiple potential solutions for this data set.

The fact that K-Stacker reports a solution compatible with the orbit of (Llop-Sayson et al. 2021) is interesting, and could suggest that this is indeed a true signal from a point-like source. However, the existence of solution 1 with a very similar (S/N)(p) tells us that solutions 2 and 3 cannot be taken as evidence of a detection: at most, they are only marginally better than the noise floor in K-Stacker.

We also notice, that all solutions find point signals around ρ ≈ 0.5″, which corresponds to the separation of the strong speckle ring in the ZIMPOL PSF (see Figs. 5 and 6), where the contrast performance for an individual night is only about C_P ≈ 4 × 10⁻⁷. A 1 R_J is expected to have a typical contrast of C_P ≈ 1.4 × 10⁻⁹ for low eccentricity orbit with a ≈ 3.5 au (Eq. (1)), and according to Llop-Sayson et al. (2021), ε Eri b is expected to have a mass of ~0.8 M_J. Therefore, even taking into account a gain of factor $\sqrt{12}$ by combining all the images, it seems unlikely that K-Stacker would be able to detect such a small planet in this data set. This is confirmed by the estimated detection limits shown in Fig. 11. Probably, the strong speckles at this separation are responsible for a non-Gaussian noise distribution which could result in the detection of high S/N noise features by K-Stacker.

A less likely explanations could be, that an unresolved dust cloud with a diameter of <0.1 au could produce a polarization signal with a detectable contrast at the level of C_P ≈ 2 × 10⁻⁷. This signal would be ten times fainter than the expected integrated signal of a dust responsible for the so-called warm 20 μm infrared excess. Thus, if the narrow dust ring model from Fig. 4 has a clumpy structure made of a few dozen components, then it seems quite likely, that the brightest one has a detectable contrast of ≈2 × 10⁻⁷. Small dust clouds were observed in scattered light for the debris disk in AU Mic (Boccaletti et al. 2015), and a body with a dust cloud could also explain the point-like scattering object Fomalhaut b (Kalas et al. 2013; Gaspar & Rieke 2020).

Fig. 10 Illustration of the 70 best orbits found by K-Stacker, after the re-optimization step. The markers correspond to the positions along each orbit at the epoch of the individual SPHERE/ZIMPOL observations. These 100 orbits are actually distributed along 3 main orbits, among which two (orbit 1 and 2) correspond to similar positions in the image.

Fig. 11 Detectability of planets as a function of the planet radius and semi-major axis for three different values of the inclination. The inclination has a significant impact on the detectability of planets, due to the influence of the phase angle on the polarization contrast.

4.7.2 K-Stacker detection limits

Determining the detection limits of K-Stacker is difficult for at least two main reasons. Firstly, as revealed by the two truly different solutions found with the grid of orbital parameters, the algorithm can converge to relatively high (S/N)(p) values (higher than the commonly used threshold of 5) without the presence of a true companion. Secondly, the algorithm takes into account the orbital motion of the planet, and therefore combines different separations and phase angle together, which makes the usual presentation of contrast as function of separation irrelevant. The true detection limits of K-Stacker can only be understood in orbital parameter space.

Nonetheless, to provide at least a rough estimate of the sensitivity achieved by combining all the available data, we created a set of “detection maps” as follows:

1. For each set of orbital parameters p of the K-Stacker grid, we calculate the amount of “missing signal” ΔS which would be required to reach a threshold of (S/N)_KS = 8, defined as $\mathrm{\Delta }S(p)=[8-(S/N{)}_{\mathrm{K}\mathrm{S}}(p)]\times \sqrt{{\sum }_{t}n(p{)}^2}$. This corresponds to the missing signal along orbit p to reach an (S/N)_KS = 8. This arbitrary threshold is chosen to be slightly higher than the K-Stacker ‘noise-floor’, which corresponds to the maximum (S/N)KS = 7.5 found by K-Stacker.

2. Injecting a companion at a known contrast C = 2 × 10⁻⁷ in the SPHERE/ZIMPOL data, and extracting the corresponding K-Stacker signal, we determine a “contrast-to-signal” conversion factor γ_t for each night.

3. From the set of orbital parameter p, we calculate the 3-dimensional positions at each SPHERE/ZIMPOL epoch t. From this, we extract both the distance to the central star d_t(p) and the phase angle α_t(p).

4. From these distances and phase angles, and using the same Rayleigh scattering as discussed in Sec. 2.1, we calculate the polarization contrasts C_P,t(p) of a 1 R_J planet using Eq. (1), and the corresponding total reference K-Stacker signal $S_{R_{\mathrm{J}}(p)}={\sum }_{\mathrm{t}}{\gamma }_{\mathrm{t}}{C}_{\mathrm{P},\mathrm{t}}(p)$.

5. The minimum detectable radius on orbit p is finally taken as $R_{min}(p)=\sqrt{\mathrm{\Delta }S(p)/S_{R_{\mathrm{J}}}}$.

In Fig. 11, we show, as a function of planetary radius and semimajor axis, the fraction of the orbits from the initial grid on which the planet could have been detected. Since the minimum detectable radius is strongly dependent on the inclination (through the phase angle), this map is presented for 3 different values of the inclination.

For i = 0^° orbits the scattering angle is always α_t(p) = 90^° and the contrast of reflecting planets is constant for e = 0 or varies with the separation $1/d_{\mathrm{p}}^2$. The fact that the planet radius for a given detectability line decreases slightly with decreasing semi-major axis indicates that the planet brightness increases a bit faster for smaller separation d_p than the detection limits. The probability for the detection of a planet changes very rapidly for i = 0°, because once the radius of a planet is large enough to be detectable then it is close to the detection limit along the whole orbit. Thus, for a = 3.5 au a planet needs to have a radius of about 2.5 R_J or a contrast of C_P ≈ 8 × 10⁻⁹ to be detected with a high probability of >95% with K-Stacker. For planets on inclined orbits, there are orbital phases where the planet is bright or faint (Fig. 2): for i ≈ 30° the planet is only bright during about 40% of its orbit and for i ≈ 70° only about 20% of its orbit and this leads to the strongly reduced detectability rates in Fig. 11 for a given planet radius and high i.

5 Search for dust

5.1 Method

The infrared SED of ε Eri indicates the presence of warm dust emitting thermal radiation within the field of view of ZIMPOL as described in Sect. 2.2, and therefore we search for the scattered, polarized radiation from this dust. We assume that the dust distribution around ε Eri did not change significantly during our four epochs of observations and we expect an azimuthal polarization signal with a positive signal for Q_ϕ(x, y) and zero signal in U_ϕ(x, y).

We process first the polarimetry of each night individually and then combine the Q_ϕ and U_ϕ maps of all twelve nights to achieve the highest possible signal to noise. Important is the correction for the radial dependence of the instrument polarization as described in Appendix A.2, because a residual instrumental signal could introduce a faint, extended Q_ϕ feature. Further, the processing of the individual nights does not include a subtraction of an extended residual Q_ϕ-pattern before derotating the individual images. Such a subtraction was applied for the search of a point source, but for an extended source this could cause significant self-subtraction, or introduce a spurious extended signal. However, the residual artefacts from the telescope spiders and the eight astrometry spots of the coronagraph are again masked as for the search of a point source. Further we can increase the sensitivity by pixel binning or image smoothing. A binning of 20 × 20 pixels (72 mas × 72 mas) would still resolve the dust structures of a narrow ring with a width of about 0.4 au.

5.2 Derived signal

The final Q_ϕ(x,y) and U_ϕ(x,y) maps for ε Eri based on 38.5 hours of SPHERE/ZIMPOL integration are shown in Fig. 12. Both maps, Q_ϕ and U_ϕ, show for ρ > 0.6″ a quite smooth wedge pattern, with alternating positive and negative signals. The morphology of the patterns is quite similar for Q_ϕ and U_ϕ, with wedges with a width of roughly ≈ 60^o but the location of the positive and negative wedges are at different position angles for Q_ϕ and U_ϕ. The Stokes Q and U maps show smooth structures with one dominant positive and one dominant negative wedge region in the field.

The patterns in Fig. 12 decrease in strength for larger separations while there is a noisy central region ρ < 0.6″ without clear structure. The maps for the fractional polarization Q_ϕ/I, U_ϕ/I, Q/I and U/I plotted in Fig. 13 show for ρ > 0.6″ to first order a wedge pattern which does not depend strongly on separation. This indicates that the polarization signals are roughly proportional to the intensity I(ρ, θ) of the stellar halo.

The wedge pattern is very weak, despite the fact that it dominates the final image. The differential polarization Q_ϕ and U_ϕ signals at ρ ≈ 1″ are of the order of ±0.01 ct/(s · px) or a surface brightness contrast of ΔSB_p = 13.6 mag arcsec⁻². The pixel count rates are about five times lower at the border of the field of view. The fractional polarization Q_ϕ(x, y)/I(x, y), U_ϕ(x, y)/I(x, y) relative to the stellar halo is for ρ > 0.6″ at the level of ±0.01%, while the noise in the central region has an amplitude of ±0.03%.

The rather constant signal in fractional polarization could indicate that a large part of the observed pattern is introduced by a cross-talk I(x, y) → Q(x, y), U(x, y) from the intensity halo of the stellar PSF. We find in each individual night roughly the same Q_ϕ(x, y)/I(x, y) and U_ϕ(x, y)/I(x, y) pattern for ρ > 0.6″. This is also true for the first half and the second half of a given night, where parallactic angles and altitudes for the telescope are different. This polarization pattern must therefore rotate together with the sky field on the detector. We expect, that dust scattering in the ε Eri system would only produce a ring-like or disk-like signal in the Q_ϕ image and none in U_ϕ and not a wedge pattern as in our data. Therefore, we must consider other potential sources for the obtained signal. Interstellar polarization could produce in our data a field independent polarization offset. However, for nearby stars the interstellar polarization is small, and indeed for ε Eri only a polarization of Q/I = +0.0028% and U/I = −0.0012% was measured by Cotton et al. (2017) and we obtained very similar values from our data for the integrated polarization signal (Appendix A.2). This signal is about ten times weaker than the amplitude of the measured polarization patterns shown in Fig. 13 and therefore we can exclude interstellar polarization as significant contributor to the residual pattern.

We think, that the obtained pattern results possibly from an instrumental effect, most likely related to the half wave plate HWP2, which is inserted in the beam for polarimetric observation. This component is rotating the polarization position angle from the sky synchronously with the field rotation into the sky coordinate system in the detector plane. A systematic effect introduced by this component would be rotating with the field, and not be averaged out by the image derotation applied in the reduction, and therefore also be constantly present in all data sets.

One might hope, that the origin of the unexplained polarization pattern can be understood and perhaps calibrated if more very deep polarimetric imaging data with SPHERE/ZIMPOL are taken and analysed.

We would like to note, that we obtain a positive net Q_ϕ signal at the level of Q_ϕ/I_★ = (4.6 ± 1.0) × 10⁻⁶ when we integrate Q_ϕ from 0.6″ to 1.6″ and normalize this to the total flux of the star. This value is about 3.5 times higher than U_ϕ/I_★ = (1.3 ± 0.6) × 10⁻⁶ obtained for the same integration region. This net Q_ϕ/I_★ signal could originate from circumstellar scattering and it has the strength expected from the models presented in Sec. 2.2. The unexpected and strong wedge pattern seen for Q_ϕ(x, y) and U_ϕ(x, y) casts some doubts on this interpretation because the weak integrated Q_ϕ/I_★-signal could simply be a net effect of systematic errors.

A polarization offset, as introduced by interstellar or instrumental polarization would create constant offsets for the fractional Stokes maps Q/I and U/I shown in Fig. 13 and quadrant patterns in the Q_ϕ(x, y) and U_ϕ(x, y) images but only a very small net contributions for the integrated Q_ϕ and U_ϕ. Similar arguments hold for an uncorrected beamshift effect, which produces a gradient in the Stokes Q(x, y) and U(x, y) and for Q_ϕ(x, y) and U_ϕ(x, y), but without a significant impact on the integrated Q_ϕ and U_ϕ. Higher order effects are required to create offsets in the integrated Q_ϕ and U_ϕ signals but this needs still to be investigated.

One may suspect that unidentified instrumental effects may introduce a negative Q_ϕ-signal as likely as a positive Q_ϕ signal, and a strong U_ϕ signal (|U_ϕ| > |Q_ϕ|) as likely as a weak U_ϕ signal (|U_ϕ| < |Q_ϕ|). The fact that we obtain a positive Q_ϕ signal which is a few times stronger than the absolute value of |U_ϕ| as expected for a real circumstellar scattering signal should therefore attract our attention. This could be a real signal, but also an instrumental effect. The latter case would be an unfortunate coincidence, at the level of one out of four possibilities, that the introduced spurious signals Q_ϕ and U_ϕ behave as expected for a weak signal from circumstellar scattering.

Fig. 12 Final polarimetric maps Qϕ, Uϕ, Q, and U derived from the 38.5 h of the SPHERE/ZIMPOL integration of ε Eri. The colour scale gives the surface brightness signal in counts/(s · px), where one pixel is 3.6 mas × 3.6 mas.

Fig. 13 Final fractional polarization maps Qϕ/I, Uϕ/I, Q/I, and U/I derived from the 38.5 h of the SPHERE/ZIMPOL integration of ε Eri.

Fig. 14 Convolved disk models from Sec. 2.2 for ε Eri injected in the S/N map for the Qϕ observations. The noise per pixel is defined by the standard deviation along the 12 nights in the image cube.

5.3 Disk detection limits

We can search in the final Q_ϕ or Q_ϕ/I maps of ε Eri in Figs. 12 and 13 for a ring-like or disk-like structure from circumstellar dust scattering on top of the instrumental pattern. A careful inspection reveals no obvious such structure in these maps. However, we can estimate rough sensitivity limits by inserting the calculated polarized flux Q_ϕ(x, y) for the four disk models described in Sec. 2.2 into the data. Figure 14 shows S/N maps, where S is the average Q_ϕ(x, y) signal with the inserted disk model and N is the standard deviation between the 12 individual nights’ Q_ϕ(x, y) images. The standard deviation is large for separations close to the star and small for large angular separations. This type of map allows to inspect the full field of view with the same scaling.

For the narrow ring model the polarized intensity is easily visible in our data, despite the systematic Q_ϕ pattern. The presence of the disk is more difficult to recognize for the broad ring models. Only the high inclination systems with i = 60^° and i = 75° have disk regions with sufficient surface brightness to be recognized. The disk model with i = 30° is hardly visible despite the fact that the disk integrated Q_ϕ signal is 16% higher than for i = 60°, or 32% higher than for the disk with i = 75°. Important for a detection is the peak surface brightness, which is higher for more inclined disks, and the presence of sharp disk boundaries which are helpful for distinguishing between disk structure and the systematic wedge pattern. We consider the wide disk model with i = 60^° as rough disk detection limit. In the fractional polarization map the signal needs to be of the order of Q_ϕ(x,y)/I_★(x,y) ≈ +0.007% (Fig. 13) to be visible as distinct positive signal on top of the smooth ‘background’ pattern.

The detection limits can be quantified as a polarized surface brightness contrast ΔSB_p = SB_p − m_* This is obtained from the count rates per pixel ct/(s · px) using the total count rates (2.09 ± 0.07) × 10⁸ ct/s for ε Eri for an aperture with a diameter of 3″ and the pixel size of 1/77 160 arcsec². This yields contrast limits for the ε Eri data of ΔSB_p ≈ 15.2 mag arcsec⁻² at ρ = 1.25″, and ″14.7 mag arcsec⁻² at smaller (0.8″) or ≈15.7 mag arcsec⁻² at larger (1.5″) separation, respectively, as plotted in all panels of Fig. 15 in black. The comparison with the disk model predictions excludes the existence of the narrow ring model, because a radial cut through the brightest ring section plotted with an orange line in Fig. 15A, would introduce a detectable polarization signal. Also, any other compact, and therefore high SB_p dust cloud structure responsible for the measured 20 μm infrared emission, such as circumplanetary dust, can be excluded within the covered separation range. The nondetection is compatible with a low inclination disk with a large width Δr > 1 au, similar to the red lines in Fig. 15A, or dust located partly outside of the SPHERE/ZIMPOL field of view r > 5 au.

We noticed in Sec. 5.2 a weak positive net signal for Q_ϕ, if we integrate the Q_ϕ from 0.6″ to 1.6″. We can also derive an azimuthally averaged surface brightness contrast profile ΔSB_p(r) for based on the mean value for annuli with Δr = 20 px, which is shown with filled dark green circles in Fig. 15B. This signal would be compatible with the two i = 60° models, if their convolved surface brightness is also averaged in these concentric annuli. This might indicate that an extended disk with an integrated disk signal of Q_ϕ/I_★ ≈ 4.6 × 10⁻⁶ but without sharp edges could be present. The caveat is the unexplained, strong systematic pattern in the data which causes significant doubts about the nature of the integrated Q_ϕ-signal. The surface brightness contrast profile ΔSB_p(r) for the U_ϕ signal illustrates this uncertainty, as it should be zero for optically thin circumstellar scattering. Clearly, there is a lot of noise for the innermost region <0.6″. For r > 0.6″ the Q_ϕ(r) signal is typically a factor 2–3 larger than the absolute value for the average U_ϕ(r) signal, which is often used as noise indicator in imaging polarimetry of circumstellar disks.

More studies are required on the instrument polarization of SPHERE/ZIMPOL to clarify and improve the measurements for ε Eri. This is not easy, because the achieved limit of the presented data is already very deep. This is illustrated in panel C of Fig. 15, which compares the surface brightness contrast limits ΔSB_p(r) of the ε Eri observations with ΔSB_I(r) and the ‘typical’ ΔSB_p(r) of the bright circumstellar disk HD 169 142 (Tschudi & Schmid 2021) also observed with SPHERE/ZIMPOL. The contrast limit ΔSB_p(r) for ε Eri is more than 5 mag arcsec⁻² deeper than for the proto-planetary disk and therefore we are faced with previously not recognized systematic noise effects.

Fig. 15 Comparison of the observational contrast limits for the polarized surface brightness (black crosses and dotted line) with the disk model calculations. (A) Comparison with cuts through the dust disk with i = 60° for the intrinsic (dashed line) and convolved (full line) models for the narrow ring (orange) and broad ring (red) models. (B) Annuli comparison of the convolved dust models (dash dotted line) with respect to the averaged Qϕ and Uϕ observations. (C) Surface brightness comparison between ε Eri and a prominent protoplanetary disk such as HD169142 together with the intensity PSF profile of ε Eri.

6 Discussion

This work presents for the ε Eri system a very deep search for scattered light from a planet or from circumstellar dust using high-contrast imaging polarimetry in the visual spectral range with SPHERE/ZIMPOL. We achieve for ε Eri in two out of the four epochs 5 σ_{$\mathcal{N}$} polarimetric contrast limits for point sources at levels of C_P ≈ 5 × 10⁻⁸ at a separation of 0.6″, 2 × 10⁻⁸ at 1.0″, and 1 × 10⁻⁸ at 1.6″ (Fig. 8). The achieved detection limits for the surface brightness contrast for polarized light from an extended source of dust scattering is at a level of about SB_p ≈ 15mag arcsec⁻² at 1.25″. These limits are discussed in this section with respect to the expected signal for the ε Eri system considering also possible observational improvements and requirements towards a successful detection.

6.1 Search for a point source

The system ε Eri is an attractive target for pushing the detection limits for the search of reflecting planets because several studies postulate the presence of the planet candidate ε Eri b with an orbital period of about 7.3–7.6 yr based on RV and astrometric data (Mawet et al. 2019; Llop-Sayson et al. 2021; Benedict 2022). The RV measurements predicted for Nov. 2020 an orbital quadrature phase which is the best phase for a polarimetric search of reflecting planets because we can expect a larger orbital separation ρ ≈ 1.1″ and possibly the maximal polarization signal of C_P ≈ 1.2 × 10⁻⁹ assuming a giant planet with an atmosphere producing a lot of scattering polarization. For the 2019 observations the planet separation and polarization could be similar, if the orbit inclination i is low, but it could also be substantially less favourable (ρ ≈ 0.5″ and C_P ≈ 0.5 × 10⁻⁹) for a high i ≈ 80° (Sec. 2.1). Unfortunately the inclination of the ε Eri b orbit is not well known.

Our contrast limits of about C_P ≈ 1 × 10⁻⁸ can exclude for ε Eri b “exotic” models, such as a gas planet with a giant ring system R_ring ≈ 10 · R_J (e.g. Arnold & Schneider 2004), which would increase strongly the reflecting surface and therefore enhance the scattering polarization by a factor of ten with respect to the reflection from the planetary atmosphere. A strong signal from a giant circumplanetary disk requires also the right values for the disk inclination and dust scattering properties. Therefore, such a system is unlikely for the nearest, single, solar type star.

Our observations show first and foremost, that the contrast limits for the search of a faint point source get deeper by increasing the integration time by roughly $C_{\mathrm{P}}\approx 1/\sqrt{t_{\exp }}$. This was already shown in Hunziker et al. (2020) for t_exp up to 100 min. This trend continues when data from two to four nights from one observing epoch with t_exp of up to 15 h are combined, as shown in Fig. 8 and numbers given in Sec. 4.5.

We also investigated the possible improvement by combining data from different epochs using the K-Stacker software (Le Coroller et al. 2020), which combines the data based on an Keplerian orbit prediction. This method allows the combination of data from different epochs for an object with substantial orbital motion and further increases the contrast roughly according to the $C_{\mathrm{P}}\approx 1/\sqrt{t_{\exp }}\text{ law }$.

This study clarifies the possible improvements for future deep searches with SPHERE/ZIMPOL. Significant deeper contrast limits are achievable within a given observing time, if the observations are only taken under very good seeing conditions. For example the contrast is about a factor 2.5 better for night 6 with an average seeing of 0.66″ when compared to night 5 with very similar t_exp but a seeing of 1.12″. The typical seeing for our t_exp = 38.5 h of ε Eri was about 0.75″. It can be expected that the same contrast as in this work would be achievable within about half the exposure time if the average seeing of the observations would be 0.5. Also favourable for a deep contrast limit is the coverage of a large field rotation during one night, which helps to improve the averaging of the residual speckle noise in the data with angular differential imaging.

The required significance for claiming a detection is substantially reduced for follow-up observations if the position of the planet ε Eri b is known from accurate stellar astrometry of the reflex motion or from the direct detection of the planet with imaging. If the planet position is known to a precision of about 25 mas, then a 3 σ_{$\mathcal{N}$} detection is sufficient to claim a significant polarimetric signal. A 3 σ_{$\mathcal{N}$} detection limit with a polarization contrast of C_P ≈ 5 × 10⁻⁹ could be possible with the combination of all our data, if ε Eri b turns out to have a favourable orbit with a low inclination, low eccentricity, and the same high brightness in polarized light for all our epochs.

Knowing the astrometric orbit would also allow to obtain for each night only Stokes Q measurements with a position angle aligned with the Q_ϕ of the planet. This is possible with ZIMPOL and would save 50% of the measuring time as the corresponding Stokes U signal from the planet is expected to be zero.

Taking all these steps into account, a 3 σ_{$\mathcal{N}$} detection for the planet ε Eri b with a contrast of about C_P ≈ 1 × 10⁻⁹ would be possible with a well known planet orbit, using only observations taken under best seeing conditions during the best orbital phases of the planet, and measuring only Stokes Q || Q_ϕ with an integration time of about 200 hours with ZIMPOL. This seems to be technically feasible with VLT/SPHERE within the framework of an ESO large programme.

6.2 Search for extended emission from circumstellar dust

We searched for an extended polarization signal from circumstellar dust in ε Eri and achieve a contrast limit of about ΔSB_p ≈ 15 mag arcsec⁻² at a separation of 1.25″ (4 au). This limit is set by an unknown systematic noise effect. This is very unfortunate, because the statistical noise limit for our data is much lower. For a single pixel at a separation of about 1″ the photon noise of the whole data set is roughly ≈14 mag arcsec⁻² (0.005 ct/(s · px)) and this could be strongly pushed by pixel binning, for example by a factor 30 to ≈17.7 mag arcsec⁻² for an area of 30 × 30 pixels. This binning provides still a very good spatial sampling of 0.11″ × 0.11″ or 0.35 au × 0.35 au for the detection of extended dust scattering around ε Eri. At this contrast level the predicted polarization signal should clearly show up, even for very unfavourable spatial distribution and scattering properties of the inner dust in ε Eri.

Actually, we measure an integrated azimuthal polarization signal of about Q_ϕ/I_★ ≈ 4.6 × 10⁻⁶ for an annulus covering the separation range from 2 to 5 au. This is roughly at the predicted level of the simple scattering models for the warm dust in ε Eri. However, we hesitate to claim a detection without a better understanding of the dominant systematic Q_ϕ noise pattern. Clearly, a better correction or calibration of the systematic pattern shown in Figs. 12 and 13 would allow a major progress for the investigation of the inner dust in ε Eri. Also other bright targets with warm dust, such as Fomalhaut (e.g. Gáspár et al. 2023) could be investigated with high sensitivity in scattered light.

For this reason we have put quite some efforts to understand the systematic effects by using various types of alternative post-processing procedures. We also analysed deep polarimetric imaging data for α Cen A available in the ESO archive (ESO programme ID 2107.C-5008), which were taken with the same instrument mode as ε Eri. We find also for α Cen A a wedge pattern in the residual polarization image similar to Figs. 12 and 13 with roughly the same contrast ΔSB_p but a different wedge geometry (Tschudi 2023). The α Cen A data are too different to be used for a calibration of the systematic pattern in the ε Eri data. More studies are required to improve further the polarimetric sensitivity for extended sources with SPHERE/ZIMPOL.

6.3 ε Eri with the Roman Space Telescope

The Nancy Grace Roman Space Telescope, which is planned to become operational in a few years, includes a Coronagraph Instrument technology demonstrator (hereafter Roman-CGI) for the detection of the reflected visible light from cold planets and circumstellar dust (Kasdin et al. 2020; Mennesson et al. 2022; Bailey et al. 2023; Doelman et al. 2023). The ε Eri system will also serve for this instrument as important test case for the instrument performance (Douglas et al. 2022; Anche et al. 2023), because it is one of best planetary systems known for a successful detection (Carrión-González et al. 2021). Therefore, it is interesting to compare briefly the expectations of this space instrument using advanced coronagraphy with our SPHERE/ZIMPOL using a ‘basic’ stellar Lyot coronagraph and PDI for the suppression of the strong, residual speckle halo of a ground based AO instrument.

It is expected that Roman-CGI will achieve for ε Eri contrast levels of about 1 × 10⁻⁹ for the reflected intensity with an integration time of about 100 hours using the ‘wide’ field imaging mode for separations 0.45″−1.4″, in the I-band (λ_c = 825 nm) (e.g. Bailey et al. 2023). This is ideal for the planet ε Eri b for which one can expect a contrast between C_I ≈ 10⁻⁹ and 10⁻⁸ as illustrated in Fig. 2. Such data could provide the planet intensity phase curve for the I-band, an accurate orbit, and pin-down for SPHERE/ZIMPOL the best orbital phase and location for a polarimetric detection of ε Eri with follow-up observations as described above.

The Roman-CGI is also expected to be very sensitive for the polarimetric mapping of the scattered light from the warm dust in ε Eri, and it could provide a detection with an integration of about 10 min (Anche et al. 2023; Douglas et al. 2022). Similar to our ZIMPOL study a narrow disk ring would be easy to measure while a smooth, extended cloud will require deeper observations and more accurate calibrations. However, a detection of the dust in ε Eri seems to be straight forward with Roman-CGI (Anche et al. 2023), even for difficult circumstances, because of the excellent coronagraphic contrast and the high sensitivity. SPHERE/ZIMPOL could complement the Roman-CGI dust scattering polarimetry with data having significantly higher spatial resolution, provided the polarimetric calibration problem described above can be solved. Additionally, ZIMPOL multi-wavelength polarimetry could constrain dust properties based on the colour of the reflected signal.

7 Conclusions

This work demonstrates the potential of PDI with SPHERE/ZIMPOL despite the fact that we could not detect the planet ε Eri b or map the extended polarization signal from circumstellar dust with the combination of t_exp = 38.5 h of data from 12 nights spread over more than a year. For ε Eri, this provides the deepest contrast limits for a point source so far and unprecedented contrast limits for the polarimetric search for extended emission from dust scattering. This pilot study is therefore useful to understand the limits of SPHERE/ZIMPOL imaging polarimetry better and it clarifies strategies to optimize future searches of faint sources around bright stars in polarized light.

On the search of point sources. The presented imaging polarimetry of ε Eri reaches 5 σ_{$\mathcal{N}$} point source contrast limits at the level of C_P ≈ 10⁻⁸ at a separation of 1σ. It seems that the limits would just further improve with the square root of t_exp or the number of collected photons. Similar point source contrast limits were previously achieved with this instrument for the much brighter stars α Cen A and α Cen B by selecting data of a single night with perfect seeing (Hunziker et al. 2020). The ε Eri data prove that the combination of data from different nights can be used as a standard procedure to push planet detection limits using the service observation mode offered for the VLT telescopes. However, this observing strategy must consider the substantial orbital motion of potentially observable reflecting planets around nearby stars using a software that searches for the planet with a Keplerian motion prediction as described in this work for the K-Stacker software (Le Coroller et al. 2020; Nowak et al. 2018).

We also compared the chances for a successful detection of a reflecting planet with SPHERE/ZIMPOL in searching for a polarized signal or the intensity signal. We find that a planet is easier to find with polarimetry. At small separations, polarimetry is up to 30 times more efficient in suppressing the strong speckle noise than the search of the corresponding intensity signal. Because Earth (Stam 2008; Bazzon et al. 2013), Jupiter (Smith & Tomasko 1984), or Uranus and Neptune (Schmid et al. 2006b; Buenzli & Schmid 2009) show a scattering polarization of about 15% or higher in the R band, a polarimetric search of a planet with SPHERE/ZIMPOL is attractive.

The polarimetric contrast limits that were reached of about C_P ≈ 1 × 10⁻⁸ at ρ ≈ 1″ are still a factor of about ten above the expected signal for the RV-planet ε Eri b. Steps to reduce the gap between the planet signal and detection limit are discussed and we conclude that a 3 σ_{$\mathcal{N}$} contrast limit of C_P ≈ 1 × 10⁻⁹ would require about 200 hours of VLT integration time under good seeing conditions and a well-known planet orbit to optimize the observing strategy in a follow-up search. This is very demanding but technically feasible.

Contrast limits of C_P ≈ 1 × 10⁻⁹ at ρ ≈ 1″ are easier to achieve for planets around the nearest bright stars Sirius, α Cen A and B, Altair, and a few others. For these objects it is possible to measure a photon flux that is ten times higher with SPHERE/ZIMPOL without harmful detector saturation effects and to reach contrasts of C_P ≈ 1 × 10⁻⁹ within 20 hours. Moreover, because α Cen is so close, a contrast of ≈1 × 10⁻⁹ would allow the detection of a planet with a radius of roughly 0.4 R_J, because 1″ corresponds to a physical separation of only d_p = 1.3 au.

Such a programme with SPHERE/ZIMPOL (ESO programme ID 2107.C-5008) was approved following the announcement of a planet candidate in α Cen A by Wagner et al. (2021); however, because of scheduling issues during the corona pandemic, only observations were executed for 4 hours. We analysed these α Cen A observations very similar to the ε Eri data and reached a 5 σ_{$\mathcal{N}$} detection limit of C_P ≈ 8 × 10⁻⁹ at 1″ for these 4 h (Tschudi 2023). This is a similar contrast to the 38.5 h for ε Eri because of the higher photon flux and the better seeing conditions for the α Cen A data. We could not detect a point source in these data. Because the planet candidate around α Cen A has not been confirmed yet, a future SPHERE/ZIMPOL programme aiming for deeper observations would face the risks of a blind search, where the location and brightness of the planet for a given epoch are unclear and therefore could be very unfavourable if the observations are scheduled at the ‘wrong’ time.

On the search of extended emission. SPHERE/ZIMPOL is the only polarimeter regularly used for high-contrast imaging polarimetry in the visual 500–900 nm range (e.g. Schmid 2021). It is therefore very useful for the study of wavelength dependencies of the scattered radiation of circumstellar disks and shells (e.g. Ma et al. 2024; Haubois et al. 2023). For these applications the ZIMPOL performance is competitive with respect to the state-of-the-art near-infrared polarimetric modes of SPHERE/IRDIS (de Boer et al. 2020) or of GPI at Gemini (Perrin et al. 2015). For dust around bright stars R < 5 mag, ZIMPOL has the advantage that it was designed for the search of planets around very bright stars and the instrument can therefore fully exploit the photon collecting power of the VLT telescope in broad-band filters without harmful detector saturation effects.

For ε Eri we added 38.5 hours of integration for the search of extended polarized emission of dust. By averaging the signal over areas of 0.11″ × 0.11″, we could further reduce the statistical noise to theoretical contrast limits of ΔSB_p ≈ 17.7 mag arcsec⁻² at a separation of about 1″. This would be enough to see the expected dust scattering clearly from the warm dust in ε Eri inferred from the 20 μm infrared excess bump measured in Spitzer data (Backman et al. 2009). However, we find that the contrast is limited by a weak systematic noise pattern to about ΔSB_p ≈ 15 mag arcsec⁻². Even this contrast is much deeper than previous measurements of the extended polarized emission in high-contrast imaging with SPHERE and GPI. Typical contrast limits for one hour observations are about Delta;SB_p ≈ 8 to 10 mag arcsec⁻² at ρ ≈ 0.4″ for high inclination debris disks (e.g. Engler et al. 2017, 2018; Esposito et al. 2020) or ≈10 to 12 mag arcsec⁻² around 1″ for protoplanetary disks (e.g. Avenhaus et al. 2018; Tschudi & Schmid 2021). We are not aware of deep searches for polarized circumstellar emission that combined long integrations of about 10 hours or more for a very bright star to reach much fainter extended sources. The presented ε Eri data show that much deeper limits ΔSB_p ≈ 15 mag arcsec⁻² can be achieved with high-contrast imaging polarimetry and the potential of such observations should be used more often for the future investigation of faint circumstellar dust around bright stars.

Acknowledgements

C.T. and H.M.S. acknowledge the financial support by the Swiss National Science Foundation through grant 200020_181983. SPHERE is an instrument designed and built by a consortium consisting of IPAG (Grenoble, France), MPIA (Heidelberg, Germany), LAM (Marseille, France), LESIA (Paris, France), Laboratoire Lagrange (Nice, France), INAF – Osservatorio di Padova (Italy), Observatoire de Genève (Switzerland), ETH Zurich (Switzerland), NOVA (The Netherlands), ONERA (France) and ASTRON (The Netherlands), in collaboration with ESO. SPHERE was funded by ESO, with additional contributions from CNRS (France), MPIA (Germany), INAF (Italy), FINES (Switzerland) and NOVA (The Netherlands). SPHERE also received funding from the European Commission Sixth and Seventh Framework Programmes as part of the Optical Infrared Coordination Network for Astronomy (OPTICON) under grant number RII3-Ct-2004-001566 for FP6 (2004–2008), grant number 226604 for FP7 (2009–2012) and grant number 312430 for FP7 (2013–2016). A.Z. acknowledges support from ANID – Millennium Science Initiative Program – Center Code NCN2021_080. This research is based on observations obtained with ESO Telescopes at the La Silla Paranal Observatory under programme IDs: 0104.C-0178(A), 0104.C-0178(B), 0104.C-0178(C) and 106.2144.001.

Appendix A Advanced improvements for the SPHERE/ZIMPOL data reduction

This appendix describes improvements in the ZIMPOL data reduction for the very deep ε Eri observations presented in this work, which are new or go beyond the procedures described previously (Hunziker et al. 2020; Schmid et al. 2018). We had to improve our data analysis because this paper pushes the limits of the ZIMPOL performance to deeper limits.

Appendix A.1 Camera 1 readout issue

The detector of ZIMPOL camera 1 had issues with the electronics during the first eight nights of our observations. The analogue to digital converter of the left read out register produced wrong results for pixels with a count level of about 7000 ADUs as shown in Figure A.1. The speckle pattern in the coronagraphic image varies as result of the atmospheric turbulence and therefore the affected pixels changed from frame to frame because different regions had exposures levels near 7000 ADUs. Typically, there are around 1000 affected pixels located near the bright speckle ring or near the ‘coronagraph mask’. Fortunately, ZIMPOL takes data simultaneously with camera 1 and camera 2 and therefore the same image is taken twice with only a small scaling factor difference because the ZIMPOL beamsplitter sends a few percent more light to camera 2. This allowed a correction of the bad pixels in camera 1, identified with an outlier detection procedure, by replacing them with corresponding, scaled pixel values from camera 2. A corrected image is shown in Figure A.1.

This reconstruction could save all the affected camera 1 data from the first eight nights without producing spurious effects for the applied post-processing procedures. Not correcting and not including the affected data in our analysis would have reduced the effective exposure time of our programme by about 28% (56 % of the observing time is affected, however only one camera, therefore 28 % of the ‘photons’). Of course, the information of the affected ≈1000 pixels per camera 1 frame is lost and only duplicated by camera 2 data, but this loss corresponds to only ≈0.1 % of one camera 1 frame or only ≈0.05% of the total pixels for each integration registered with the two detectors. ESO has solved the issue before the 9^th night by replacing read-out electronics boards and two broken cooling fans.

Fig. A.1 Illustration of the ZIMPOL camera 1 readout issue affecting the left side of the image. Illustrated is one individual intensity image (two 5 s sub-integrations = 10 s) with readout problem before and after correction of the affected pixels.

Appendix A.2 Telescope polarization correction

For deep polarimetric observations with SPHERE/ZIMPOL it is important to apply a correction for the residual telescope polarization as described in Schmid et al. (2018). A polarization of about 4 % is introduced by the aluminium coated M3 mirror in the VLT, which is for ZIMPOL polarimetry compensated with a rotating half-wave plate and a ‘crossed’ M4 mirror. The compensation is not perfect, but the telescope effects are reduced to about p_tel ≈ 0.5 % or less. This residual polarization depends on the parallactic angle θ_para of the telescope and the measured values lie in the Q/I-U/I plane on a circle with radius p_tel with position angle θ_tel = θ_para + δ_tel (see e.g. Hunziker et al. 2020; Tschudi & Schmid 2021). The centre (q_m, u_m) of the circle can be offset from the origin (0,0) due to interstellar or intrinsic polarization of the star. We measure for ε Eri with the VBB filter in 2019 and 2020 p_tel = 0.26 ± 0.03 and δ_tel = 32.1 ± 1.5° and a centre of (q_m = 0.0032 ± 0.006 %, u_m = −0.0013 ± 0.004 %) with zero polarization as expected and in very good agreement with previous high precision measurements for this target ε Eri q_m: 0.00284 ± 0.00056 % and u_m: −0.00120 ± 0.00057 % and p ≈ 30 · 10⁻⁶, θ = 168.5° (Cotton et al. 2017). The expected polarization from light scattering by the circumstellar dust around ε Eri is much lower, less than 10⁻⁴ % (see Section 2.2). For such objects without strong intrinsic polarization, a good first order correction for the telescope polarization is obtained by the normalizations I₀ = I₉₀ and I₄₅ = I₁₃₅ for each cycle. This is equivalent to setting the integrated polarization to zero: Q = 0 and U = 0. Not correcting for the telescope polarization would introduce an I → Q, U cross talk and intensity speckles would be visible as polarized features with a relative strength at the level of the telescope polarization p_tel.

For ε Eri, we need to consider also second order effects of the telescope polarization. The two parameters describing the telescope polarization p_tel and δ_tel depend on the wavelength (Schmid et al. 2018) andthis is an issue for very deep polarimetry in the VBB filter with large bandwidth 590 – 890 nm and therefore significantly different instrument polarization for short and long wavelengths. This is shown in Tschudi & Schmid (2021) for the R′ and the I′ filter which cover roughly the short and long wavelength parts of the VBB filter, respectively. In addition, the SPHERE AO PSF for the VBB filter is a superposition of different radial profiles for the different wavelengths. This produces a 2-dimensional polarization effect which cannot be corrected with a single vector (p_tel, δ_tel) without leaving second order calibration errors.

Fig. A.2 Radial dependence of the telescope polarization parameters ptel and δtel for the Very Broad Band (VBB) filter. The plotted curves are the mean curve for all ε Eridani observations listed in Table 2.

The dominant feature for the radial dependence of the instrument polarization is the wavelength dependent location of the strong PSF speckle ring, which is defined by the control radius (20 δ/D) of the SPHERE AO system. This ring is for the R-band at a separation of ρ ≈ 0.35″ and for the I-band at ρ ≈ 0.45″ (Figure 11 in Schmid et al. 2018) and therefore the R-band instrument polarization dominates at smaller separations while the I-band polarization contributes strongly for ρ ≥ 0.45″. We could measure clearly a corresponding radial dependence of the telescope polarization p_tel(r) and δ_tel(r) with amplitudes of ±0.01 % and ±2° for the mean radial curve as illustrated in Figure A.2. The curves p_tel(r) and δ_tel(r) look similar in shape from night to night but show small variations Δp_tel < 0.02 % and Δδ_tel < 1° because of PSF variations introduced by different observing conditions. We don’t see a long term trend within the 13 months covered by our observations, but on longer timescales one should expect systematic changes because of the evolution of the coatings of the telescope mirror M3 and the first folding mirror M4 in SPHERE.

For the general case with a (q_m, u_m) offset from (0,0) we would determine the mean p_tel(r) and δ_tel(r) curves of each night to correct the telescope polarization of that night. From these profiles we construct two interpolated 2d maps to correct each polarization image pixel-wise depending on the parallactic angle of that data. In the special case of ε Eri with no offset polarization (q_m, u_m) = (0, 0), we can measure Q/I, respective U/I in every image individually and correct it to 0. We do this for each radial annuli separately to fully account for the radial dependence as described above. The advantage of this method is that PSF variations, which can happen within minutes, and higher order telescope polarization effects are also corrected.

This second order correction for the telescope polarization is not crucial for the search of a planet outside the speckle ring ρ > 0.6″ where the dependence of p_tel(r) and δ_tel(r) is smooth. At small separation (<0.5″), the residual pattern of strong speckles is slightly reduced. The second order correction should however improve the search of an extended weak signal from dust scattering. For example, the effect is clearly seen in the Q and U images of the debris disk observations of HIP 79 977 by Engler et al. (2017, Figure 3) as over-corrected central area r < 0.5″. At the time of that analysis the origin of this calibration problem was unknown.

Appendix A.3 Beamshift correction

It is important for polarimetric differential imaging (PDI) performance that the differential aberrations between the opposite polarization directions I_‖ and I_⊥ are very small so that unpolarized speckles and other PSF features cancel out by the subtraction of the two images and the polarized planet signal is easier to detect. The ZIMPOL design was optimized to reduce such differential aberrations and the two polarization states are for example recorded with the same detector pixels. However the inclined third mirror M3 of the telescope, the 45^° pupil tip-tilt mirror and the three image derotator mirrors mainly introduce a, wavelength and telescope pointing dependant, differential beamshift of up to 0.3 pixels (or ≈1 mas) between I_‖ and I_⊥ Schmid et al. (2018); Hunziker et al. (2020). For a given pointing of the telescope, meaning the same altitude and parallactic angle, the beamshift effect is the same and therefore we can sort the ε Eridani data and derive the beamshift parameter as a function of their local sideral time (LST). There are many beamshift parameters as the shifts are different in X and Y directions in the image, different for the polarization images Q⁺, Q⁻, U⁺, U⁻, different for camera 1 and camera 2 and different for the 0 and π phase images. To measure the beamshift accurately an unsaturated point source is required as available from the regularly taken non-coronagraphic observation using the ND2 filter. The used Lyot coronagraph (V_CLC_MT_WF) has a slightly transparent mask and under good conditions it is possible to see the stellar PSF peak through the mask and to measure the beamshift also in the corona-graphic observations. In the top panel of Figure A.3(A) the non-coronagraphic measurements of the ΔX(0 phase) beamshift parameters for Q⁺, Q⁻, U⁺, U⁻ are displayed. The same parameters measured for the coronagraphic data can be seen in the middle panel of Figure A.3(B). It is obvious that the dispersion is larger for the measurements with coronagraph, however there are many more images and no time gaps between the images, but sometimes the determination of the beamshift fails for the coronagraphic data (many of these points are far outside the showed y-axis in the middle panel of Figure A.3(B)). Unfortunately, there exists no model for SPHERE/ZIMPOL which could predict the beamshift effects. Therefore we use all the existing ε Eri data and derive the shift corrections from the best fit to the data. For this we calculate a Gaussian process fit, apply an outlier detection and removal method and iterate the procedure a few times (see bottom panel of Figure A.3(C)). The fitting procedure is particularly important for data taken under bad conditions (seeing > 1’’). All the images are visually checked after applying the beam shift correction.

Fig. A.3 Beamshift correction parameters as a function of the pointing (local sidereal time). Illustrated as example is the shift in X direction for camera 2 and 0 phase images. (A): Beamshift as measured in the unsaturated (ND2) PSF images. Different colour represent the different polarization images Q+, Q−, U+ and U−. (B): Analogous figure as measured in the coronagraphic images. The measurement accuracy is reduced especially for bad observing conditions. (C): Example for the iterative process of Gaussian process fitting and outlier detection to calculate the fit of this parameter for Q−.

Appendix A.4 Chargetrap correction

Both ZIMPOL detectors consist of alternating open and masked rows (512 rows each) with 1024 pixels. In the polarimetric ZIMPOL modes the charges created in the illuminated rows are shifted up and down in synchronization with the polarimetric modulation (frequency depends on fast or slow polarization mode) (Schmid et al. 2018). The final detector frame consists of an “even rows” subframe with one polarization state I_⊥, and an “odd rows” subframe for the opposite polarization state I_‖. In column direction the subframes are then interpolated in a flux conserving manner to create 1024x1024 images. Although the charge transfer efficiency of the CCDs are better than 99.9995 % there exist pixels which do not shift electron charges perfectly. For example a pixel can block one electron which is not down shifted but in the following up-shift it is transferred. Because of the fast (de)modulation such a trap can cause a hole of many electrons in one polarization state and a corresponding spike in the other. To get rid of this effect in ZIMPOL polarimetric modes there are always an even number of subintegrations and in every second subintegration the charge shifting is reversed with respect to the polarization modulation. The images from zero to π phase subintegration are then combined to create a double difference in which the charge trap effects are cancelled (Schmid et al. 2018; Gisler et al. 2004; Schmid et al. 2012). Unfortunately, for the intensity frame I = + I_⊥ + I_‖ derived from the polarimetric data the charge traps do not vanish with the combination of the subintegrations. The charge traps produce a negative-positive pattern with a very specific appearance in column direction. First a pixel with too many counts, then an interpolated ‘neutral’ pixel and then a pixel with too few counts or the other way around. Sometimes secondary and tertiary pixels are also slightly affected. We search for this pattern in each image and correct it by shifting counts from a ‘spike’ to a corresponding ‘hole’ in a flux conserving manner as illustrated in Figure A.4). More than 99 % of the charge traps can be recognized and corrected with this method. This correction is especially helpful for polarimetric p2 mode when the field is fixed and one has only a few long exposures.

Fig. A.4 Illustration of the charge trap correction for ZIMPOL intensity images [ct DIT−1 px−1] obtained in polarimetric mode. Article number, page 24 of 25

Appendix B Detection maps

Fig. B.1 Polarized light Qϕ detection maps showing the Gaussian significance of the individual pixels. North is up, east is to the left, and the image field of view is the same as in Figure 7. Top: Detection maps of all the individual nights. Bottom: Detection maps of the four epochs.

References

All Tables

All Figures

	Fig. 1 Illustration for a possible orbit (i = 78°) and signal strength of ε Eri b. The blue colour indicates the normalized strength of the polarized intensity and the red colour is analogously for the intensity. The four observing epochs are marked with black dots and the points in the orbit with a minimal, zero, and maximal RV are shown in grey.
In the text

	Fig. 2 Planet model contrast as a function of time for the intensity (upper panel) and the polarized intensity (lower panel) for the orbit with an inclination of 78° (black curves) as in Fig. 1 and for an inclination of 30° (blue curves). The dots show our four observing epochs, while the octagon stands for earlier ZIMPOL observations in Hunziker et al. (2020).
In the text

	Fig. 3 Relative flux distribution for a narrow dust ring model in ε Eri with r = 4 au based on the warm component in the infrared excess. Top: thermal emission. Middle: expected scattered intensity 1disk for g = 0.6. Bottom: polarized intensity Qϕ for g = 0.6. The model images are convolved with the PSF of SPHERE/ZIMPOL, normalized to their flux maximum and aligned with the sky coordinate axes.
In the text

	Fig. 4 Predicted polarized flux Qϕ [ct/(s · px)] produced by scattering in a narrow and a broad dust ring model for ε Eri with an inclination of i = 60° (top row), and for broad rings with i = 30° and i = 75° (bottom row). All models would produce the measured infrared excess of Lwarm/L★ = 3.3 × 10−5.
In the text

	Fig. 5 One polarimetric cycle texp = 240 s corrected for the telescope polarization: the top half shows the intensity I and bottom half the polarized intensity Qϕ with corresponding colour scales in ct/(px· DIT). The counts inside 0.9″ are scaled down for a better visibility by a factor of six for the intensity and a factor of three for Qϕ.
In the text

Fig. 6 All data of night 9 (texp = 9360 s), derotated, averaged, and corrected for telescope polarization: the top half shows I and the bottom half Qϕ with corresponding colour scales in ct/(px · DIT) as in Fig. 5. Artificial point sources are inserted at different separations and contrasts CI and CP are changed with position angle as indicated. Counts inside 0.9″ are scaled down for better visibility with a factor of six for I and a factor three for Qϕ.

In the text

Fig. 7 All data of night 9 (texp = 9360 s), derotated, PSF subtracted, masked regions for the Qϕ, averaged and corrected for telescope polarization: the top half shows I and the bottom half Qϕ, with corresponding colour scales in ct/(px · DIT) as in Fig. 5. Artificial point sources are inserted at different separations and contrasts CI and CP are changed with position angle as indicated. Counts inside 0.9″ are scaled down for better visibility with a factor of six for I and a factor three for Qϕ.

In the text

	Fig. 8 Sensitivity expressed as contrast curves with 5 σ$\mathcal{N}$ Gaussian significance. Top: individual nights for the intensity CI = Ip/I★ (upper curves), the polarized light CP = pp · Ip/I★ (lower curves). Bottom: the four epochs in colour in comparison to the individual nights in grey.
In the text

	Fig. 9 Ratio between polarization CP and intensity CI contrast limits for the full data set, illustrating the advantage of PDI for the speckle suppression. Planets with a fractional polarization p > CP/CI should be easier to detect in the polarimetric data.
In the text

	Fig. 10 Illustration of the 70 best orbits found by K-Stacker, after the re-optimization step. The markers correspond to the positions along each orbit at the epoch of the individual SPHERE/ZIMPOL observations. These 100 orbits are actually distributed along 3 main orbits, among which two (orbit 1 and 2) correspond to similar positions in the image.
In the text

	Fig. 11 Detectability of planets as a function of the planet radius and semi-major axis for three different values of the inclination. The inclination has a significant impact on the detectability of planets, due to the influence of the phase angle on the polarization contrast.
In the text

	Fig. 12 Final polarimetric maps Qϕ, Uϕ, Q, and U derived from the 38.5 h of the SPHERE/ZIMPOL integration of ε Eri. The colour scale gives the surface brightness signal in counts/(s · px), where one pixel is 3.6 mas × 3.6 mas.
In the text

	Fig. 13 Final fractional polarization maps Qϕ/I, Uϕ/I, Q/I, and U/I derived from the 38.5 h of the SPHERE/ZIMPOL integration of ε Eri.
In the text

	Fig. 14 Convolved disk models from Sec. 2.2 for ε Eri injected in the S/N map for the Qϕ observations. The noise per pixel is defined by the standard deviation along the 12 nights in the image cube.
In the text

Fig. 15 Comparison of the observational contrast limits for the polarized surface brightness (black crosses and dotted line) with the disk model calculations. (A) Comparison with cuts through the dust disk with i = 60° for the intrinsic (dashed line) and convolved (full line) models for the narrow ring (orange) and broad ring (red) models. (B) Annuli comparison of the convolved dust models (dash dotted line) with respect to the averaged Qϕ and Uϕ observations. (C) Surface brightness comparison between ε Eri and a prominent protoplanetary disk such as HD169142 together with the intensity PSF profile of ε Eri.

In the text

	Fig. A.1 Illustration of the ZIMPOL camera 1 readout issue affecting the left side of the image. Illustrated is one individual intensity image (two 5 s sub-integrations = 10 s) with readout problem before and after correction of the affected pixels.
In the text

	Fig. A.2 Radial dependence of the telescope polarization parameters ptel and δtel for the Very Broad Band (VBB) filter. The plotted curves are the mean curve for all ε Eridani observations listed in Table 2.
In the text

Fig. A.3 Beamshift correction parameters as a function of the pointing (local sidereal time). Illustrated as example is the shift in X direction for camera 2 and 0 phase images. (A): Beamshift as measured in the unsaturated (ND2) PSF images. Different colour represent the different polarization images Q+, Q−, U+ and U−. (B): Analogous figure as measured in the coronagraphic images. The measurement accuracy is reduced especially for bad observing conditions. (C): Example for the iterative process of Gaussian process fitting and outlier detection to calculate the fit of this parameter for Q−.

In the text

	Fig. A.4 Illustration of the charge trap correction for ZIMPOL intensity images [ct DIT−1 px−1] obtained in polarimetric mode. Article number, page 24 of 25
In the text

	Fig. B.1 Polarized light Qϕ detection maps showing the Gaussian significance of the individual pixels. North is up, east is to the left, and the image field of view is the same as in Figure 7. Top: Detection maps of all the individual nights. Bottom: Detection maps of the four epochs.
In the text