Imaging low-mass planets within the habitable zone of α Centauri

Giant exoplanets on wide orbits have been directly imaged around young stars. If the thermal background in the mid-infrared can be mitigated, then exoplanets with lower masses can also be imaged. Here we present a ground-based mid-infrared observing approach that enables imaging low-mass temperate exoplanets around nearby stars, and in particular within the closest stellar system, α Centauri. Based on 75–80% of the best quality images from 100 h of cumulative observations, we demonstrate sensitivity to warm sub-Neptune-sized planets throughout much of the habitable zone of α Centauri A. This is an order of magnitude more sensitive than state-of-the-art exoplanet imaging mass detection limits. We also discuss a possible exoplanet or exozodiacal disk detection around α Centauri A. However, an instrumental artifact of unknown origin cannot be ruled out. These results demonstrate the feasibility of imaging rocky habitable-zone exoplanets with current and upcoming telescopes.


Artifact Modeling I. Persistence Stripes
The most significant systematic artifact is detector persistence accumulated as a result of photons collected during the chopping sequence. These line-like features rotate with the field of viewresulting in reduced sensitivity throughout the entire image and greater false-positive potential in the areas directly between the stars. Either an empirical or a physically motivated model is needed to correct for this effect−similar to the charge trapping and delayed release models recently developed for HgCdTe detectors (e.g., 1,2). We constructed two persistence models, which are shown in Supplementary Fig. 1. The first is a simple empirical model that includes a superposition 10 of lines convolved with the instrumental PSF. The second is a higher-fidelity model in which we simulated the chopping kinematics and charge accumulation. In this model, we assumed a persistence response that is proportional to the signal in the PSF and included a free parameter to scale the intensity of the persistence pattern of α Centauri B with respect to that of α Centauri A to remove the major contributions of non-linearity. We translated the PSFs of α Centauri A and B 15 across a simulated image plane for both strokes of the chopping sequence, including a sub-pixel offset and corresponding re-centering behind the coronagraph at the end of each stroke. We do not account for a potential time decay of the signal and retain all accumulated intensity in the pattern. We subtracted the persistence model of the second stroke from that of the first, and then subtracted the combined model from the individual frames. 20 To fit the parameters, we measured the positions of the off-axis PSFs of α Centauri A and B in each individual image and performed a grid search with velocity (degenerate with intensity, which we did not vary) and the offset from the center of the coronagraph as free parameters. We assumed that the offset for the second stroke is opposite that of the first to mimic a small error in the chopping angle. For each set of parameters within the grid, we constructed and subtracted the 25 persistence model from each individual image prior to PSF subtraction, and repeated the remaining steps of the data reduction-primarily including high-pass filtering and PSF subtraction via the KLIP algorithm. We adopted the parameters that minimized the standard deviation of pixel values within the final combined image in an annulus from 10−45 pixels (∼0.45−2 arcsec).
Both models are able to qualitatively reproduce the persistence pattern observed in the data 30 (see Supplementary Fig. 1). Following model subtraction, a greater speckle density is present in the areas covered by persistence. The speckle intensity distribution in this region is clearly nonuniform, with several maxima present primarily within the persistence pattern from α Centauri A. Therefore, we refrain from characterizing sources in the areas impacted by persistence (between the dashed lines in Fig. 4). Additionally, the residual noise in the areas not impacted by persistence 35 is improved following subtraction of the simulated model prior to PSF subtraction. This is because the persistence rotates with the field of view and biases the PSF subtraction. Accordingly, we utilize the images generated by the removal of the simulated persistence model prior to high-pass filtering and PSF subtraction for the proceeding analysis. The detection and properties of the candidate C1 are not affected by the choice in persistence model (including without subtracting 40 any model).

Artifact Modeling II. Optical Ghosts
The second most significant artifact that appears in the final images are two dark arcs−likely bright features in the off-axis PSF of α Centauri A−that come within 1" of the coronagraphic residual.
These features are likely faint optical ghosts that are introduced by the spectral filter, and which rotate with the parallactic angle to create arc-like features in the final image. The presence of these arcs weakens the overall image sensitivity. Similar to our handling of the detector persistence stripes, we created an empirical model of the dark arcs as circular segments of negative flux and subtracted this from the images before PSF subtraction. We repeated the process to optimize the parameters of the arcs: including radius, brightness, and width, while keeping the center fixed to 5 the off-axis position of α Centauri A. Following the removal of the arcs, the image background within ∼1" of the coronagraphic mask is nearly flat, with the exception of the persistence stripe residuals (see Supplementary Figs. 1-3).
Plate Scale and True North Calibration 10 We utilized the orbit of the binary as reported in (3)

Background-Limited Sensitivity
To compute the background-limited sensitivity, we measured the standard deviation of pixel values in a region of the detector far from α Centauri A and B, and compared this to the flux measured in 20 a four-pixel radius aperture centered on the off-axis PSF of α Centauri A. We converted this ratio to a 3-σ sensitivity estimate by multiplying by a factor of (3 × √16π) to account for the aperture area, and arrived at 5×10 −7 contrast with respect to α Centauri A, or about 65 μJy (the N-band flux of α Centauri A is 129 Jy, 4). The background near the stars is higher (by a factor of two at 1 arcsec) due to the increased photon noise and the glow of the AGPM (5).

Simulated Planetary Contrast vs. Separation Curves
We first calculated the radiative equilibrium temperature of each simulated planet on a circular orbit as a function of its orbital semi-major axis (a), its Bond albedo (AB), the stellar radius (R), 30 and stellar temperature ( !"#$ ). This is given by Supplementary Equation (1), which is obtained by setting the planet's luminosity to be equivalent to its total input of irradiant energy from the star (6).
We assumed !"#$ = 5800 K and R = 1.2 R⊙ for A, and !"#$ = 5250 K and R = 0.86 R⊙ for B (48). We then included additional heat as a fraction of the planet's equilibrium temperature via Supplementary Equation (2), and computed Blackbody spectra for the star and planet. 40 We then created simulated photometric measurements by integrating the spectra between 10-12.5 µm in steps of 0.1 µm. These synthetic photometric measurements were used to calculate planet-star contrast ratios. We repeated the analysis for planets of different radii, Bond albedo, 45 orbital semi-major axes, and additional heating conditions to compute the simulated contrasts that are utilized to estimate the properties of C1 and also for comparisons in Figs. 1, 4, & 5.
Despite its relative simplicity, this approach is sufficient to gain an understanding of the sensitivity of our single-band photometric measurements. Specifically, uncertainties in planetary 5 properties (e.g., molecular content and sources of additional heating) preclude the usefulness of a more sophisticated model. Our model does not include a wavelength-dependent model of molecular absorption and emission processes. Instead, absorption processes are incorporated in bulk via the Bond albedo. For Earth's atmosphere, the dominant absorbers in the λ~8-13 µm range are O3 at λ < 10 µm and CO2 at λ > 12.5 µm (7). Therefore, the overall level of absorption for 10 Earth-like planets is expected to be low within the λ = 10-12.5 µm range of the NEAR filter (8). For Neptune, the dominant feature within this range is C2H6 emission at 12.2 µm, which is an order of magnitude above the thermal continuum (9). For Jupiter, NH3 absorption at ~11 µm is present at a comparable equivalent width to the also present C2H6 emission (10). Such planetary atmospheres at ~1 au would have comparable equilibrium temperatures to Earth, which will 15 change the overall ratios of these constituents and their spectral contributions. However, more detailed modeling of such atmospheres also suggests that blackbody spectra are reasonable approximations for λ = 10-12.5 µm (11). This approach also does not account for specific processes that could be responsible for producing additional heating or non-thermal emission and includes the sum of all of these in a simple fractional parameter, %2(*) . If the solar system's giant 20 planets are an accurate indication of the larger population, then these sources of additional heating can be expected to raise the thermal infrared planetary fluxes by up to an order of magnitude compared to the predictions of radiative thermal equilibrium (corresponding to %2(*)~0 .8, which is above the maximum of %2(*)~0 .5 considered here; see 12, 13). 25 While these models carry a significant uncertainty in the maximum luminosity, they establish useful lower limits in cases of low additional heat and Bond albedos representative of the solar system planets (AB ~ 0.3 for Earth, 6). Finally, we verified the predicted contrasts by comparing to the Earth's brightness at λ ~ 10 µm (7). A planet around α Centauri A would experience a similar irradiance as the Earth if it were at a distance of ~1.2 au, or ~0.9 arcsec. Scaled to a radius of 1.7 R⊕ and shifted to the distance of α Centauri (1.3 pc), such a planet would appear with a contrast of ~5×10 -7 . Our corresponding model of a super-Earth with %2(*) = 0.1, with the same planetary radius (1.7 R⊕), and at the same separation from the star (0.9 arcsec) has a contrast of ~4×10 -7 . This serves as a useful consistency check of our models and assumptions.

Simulated Planet Injection and Retrieval Tests
To compute the image sensitivity, we performed extensive point source injection and recovery tests (see Figs. 3 & 4). For the injected PSF, we utilized the off-axis PSF of α Centauri A and included a scaling correction for the transmission of the AGPM (14). The PSFs were injected prior to high-pass filtering, after which we completed the remaining data processing steps (PSF 40 subtraction, etc.). We measured the signal within an aperture of two pixels in diameter centered on the location of the injected source. This diameter corresponds to ~0.3×FWHM, which is smaller than typically used. This provides a larger number of apertures that are non-contaminated by persistence. To measure the noise, we blocked the regions of the image affected by detector persistence and a 12-pixel wide box centered on the injected source and calculated the mean and 45 standard deviation of fluxes within all non-overlapping apertures at the same separation. We converted these measurements into a signal-to-noise ratio (SNR) estimate (15), adopting SNR=3 as a detection threshold, as this resulted in all of the automatically retrieved sources also being identified by eye in the images (see Fig. 3). We caution that as a result of the residual structures in the images, the statistics are non-Gaussian, and probabilities derived from Gaussian statistics should not be applied to these measurements. We repeated the analysis with injected sources along 5 ten different position angles and from separations of 0.3 to 1.5 arcsec. We adopted the average sensitivity level over the tested position angles as the sensitivity for each separation in the contrast curves shown in Fig. 4. For the completion analysis (see Fig. 5 and following subsection), we adopted the maximum brightness of the recovered sources as the sensitivity for each separation. 10

Completeness Analysis
We converted the sensitivity analysis into completion estimates using a Monte Carlo simulation to draw randomly sampled orbital parameters. For different sets of uniformly sampled planet radius and semi-major axis, we sampled the argument of periapsis (ω ∝ unif[0 , 360]; this parameter is degenerate with the time of closest approach for the circular orbits considered here), angle of 15 ascending node (Ω ∝ unif[0 , 360]), and inclination (i ∝ sin i) and then determined the planet's separation and position angle on the sky. The radiative equilibrium temperature was determined assuming a Bond albedo of 0.3 and then increased by 10% and 50% to account for different amounts of internal/additional heat from the planet. The temperature was converted to an N-band contrast following the previous subsection. This contrast was then compared to the maximum 20 point-source brightness at that separation that would have been detected based on simulated pointsource injections (see above). The completeness maps are shown in Fig. 5. We also tested restricting the plane of the planetary orbits to the orbital plane of the binary, and also allowed for eccentric orbits, and found that these did not significantly alter the completeness maps.

25
Properties of The Candidate Detection 'C1' Based on a similar analysis to the simulated planet injection and recovery tests, C1 is detected with SNR~3.1 (using apertures with diameter of two pixels, resulting in 27 noise measurements not contaminated by the persistence artifact). Using apertures with diameter of 1×FWHM, corresponding closely to the definition of SNR in (15), gives SNR~3.5, although with only ten 30 independent noise measurements after rejecting seven due to persistence contamination. Again, we caution that probabilities derived from Gaussian statistics should not be applied to this SNR estimate. To measure the brightness and position of C1, we performed negative point source injections over a grid in radius, position angle, and contrast to minimize the residual squared flux in the total image area (see the bottom-left panel of Supplementary Fig. 4). We injected the point 35 sources in the data processing prior to high-pass filtering and PSF subtraction and performed a similar analysis with both pipelines. The findings were consistent, and we averaged the results to arrive at ρ = 0.85 ± 0.05 arcsec, θ=228.9º ±3.3º E of N, and FC1/FαCenA = 9.3 × 10 −6 ± 3.1×10 −6 , corresponding to FC1=1.2±0.4 mJy, after correcting for the ∼20% throughput of the coronagraph, and by estimating the astrometric uncertainties as the FWHM divided by the SNR. Assuming that 40 the emission comes from a planet with a Bond albedo of ≤0.3 (similar to Earth), and that the planet's temperature is between 100−150% of its radiative equilibrium temperature (Teq~300 K) to account for the possibility of additional sources of heat (such as greenhouse heating or residual heat from the planet's formation), and assuming blackbody spectra, this brightness and projected separation correspond to a radius of ∼3.3−11 R⊕. This is a plausible radius and level of additional 45 heating compared to the ice− and gas-giants within our own solar system (e.g., 12, 13).

Pre-imaging for Background Sources
We first ruled out the possibility that C1 may be a background source by utilizing K-band (∼2 μm) images taken in 2009 (3). The faintest source identifiable in these earlier images is a K∼13.6 star that approaches the trajectory of α Centauri B in 2021. There is nothing comparable in brightness 5 and proximity in mid-2019. Given that the NEAR campaign images are sensitive to N∼12−12.5, and that even the reddest stars have neutral mid-to-near infrared colors, we can exclude the possibility that C1 is a background star. If C1 is a background galaxy, it would likely have a rest peak wavelength of λ ~ 1 μm and would thus be redshifted by Z ~ 10. Such a galaxy would be among the most distant galaxies observed, and unlikely to appear projected on the sky within an 10 arcsecond of α Centauri A and with such an apparent brightness. In all of the HST legacy fields spanning >800 arcmin 2 , only a handful of such galaxies were detected (16), and all have near-IR magnitudes >25, which is much fainter than our observations would have been able to detect.
Exozodiacal Dust Disk Modeling 15 We simulated an exozodiacal disk at the central observing wavelength of ∼11.25 μm using zodipic (17). The code enables simulating both the disk and the star simultaneously, which, when convolved to the resolution of the observations, enables straightforward model flux calibration. We assumed a 79º inclination to be coplanar to the orbit of the binary (3) and fit the position angle and small central offset. Finally, we ran a grid search over the amount of dust in the model disk to 20 match the brightness of C1. We subtracted each model from the data prior to high-pass filtering and PSF subtraction and adopted the model that minimized the squared residuals in an aperture of 10−45 pixels from the center as the best-fit. The results are shown in Supplementary Fig. 4 in comparison with the residuals from the point source subtraction.

25
Astrometry and Orbital Modeling To measure the astrometry of C1, we combined back-to-back nights, and re-performed minimization of the image residuals by negative synthetic planet injection tests. The astrometric precision was estimated as 1/3 of a beam diameter, or ∼8 mas. The measurements are shown in Supplementary Fig. 5. To check consistency of the measurements with orbital motion, we 30 generated 10,000 sample orbits using the OFTI method (18) in the orbitize! package (19) with the default settings (i.e., uniform priors). One hundred randomly selected sample orbits are shown in Supplementary Fig. 5 along with the retrieved orbital parameters. The data are consistent with a co-planar and dynamically stable orbit interior to the α Centauri AB binary (a<3 au and i ∼70º). While the astrometric data are consistent with a Keplerian orbit, we caution that this does not 35 constitute proof of orbital motion, as the data are also consistent with no motion.
The measurements show a tentative change of ∆ρ ~ −0 .1 arcsec (i.e., toward the star) and ∆PA ∼ 8º from 2019-05-24 to 2019-06-11, which establishes a useful vector to check against the data taken in late June. Projected forward in time, C1 would likely be at a separation of ρ ∼ 0.6−0.7 arcsec on 2019 June 27. At this position, C1 would overlap the area of the image that is strongly 40 impacted by the persistence stripes, and thus would not likely be detected. While C1 is not detected on this night, forward modelling of the anticipated signal assuming zero motion suggests that such a source would likely have been detected (see Supplementary Fig. 6). The non-detection of C1 on June 27, compared to its identification in ∼85% of the single-night images taken prior to this date, tentatively supports the planetary hypothesis over those invoking exozodiacal dust or an unknown 45 instrumental artifact, which would remain in a stable position throughout the campaign.   Supplementary Fig. 3. NEAR data processed by two completely independent data processing pipelines. The primary data processing (a) implements removal of the known systematic artifacts, whereas the persistence stripes and dark arcs are present in b. Excluding the known instrumental artifacts, C1 is the only source with a repeated detection in both images. Supplementary Fig. 4. C1 subtracted via both an exoplanet and exozodiacal disk models. a shows the data without C1 subtracted for comparison. c shows the subtraction of C1 with a point source (i.e. exoplanet) model. b shows the model disk, which is inclined by 79º from the plane of the sky, extends from 0.8−1.2 au, includes 60 zodis of dust, and is offset by ∼0.3 au in the direction 5 of C1. d shows the image with the disk subtracted. Both models match the data reasonably well. Supplementary Fig. 5. Astrometric data and orbital fit for C1 generated from the NEAR campaign. The astrometric data (a, b) are consistent with a bound orbit with a semi-major axis of ~1.1 au and an inclination of 65º±25º (c). Error bars in a and b represent the estimated astrometric uncertainty of 1/3 of a beam diameter, or ∼8 mas. While the nature of the candidate is not yet