Improved Companion Mass Limits for Sirius A with Thermal Infrared Coronagraphy Using a Vector-apodizing Phase Plate and Time-domain Starlight-subtraction Techniques

Thermal IR is challenging to observe from the ground due to the rapidly varying background illumination from the sky and noncryogenic optics. Unlike stellar speckle noise, background counts are best modeled in the original pixel frame, before image alignment. We used the sky frames we recorded to construct a principal-components basis. To model and subtract the background signal in a given science frame, we first identified pixels within a certain radius of the twin gvAPP-180 PSFs or a number of Clio-specific artifacts and excluded them (Figure 1). To mitigate impacts on sensitivity, only the small subset of pixels with negligible light from Sirius A were used to estimate the background in the reconstruction process. We then performed a least-squares fit of the basis vectors to the pixels that remained, and used those coefficients to construct a background frame from the full unmasked principal component basis vectors (Long et al. 2018; and independently in Hunziker et al. 2018).

We held back 25% of the sky frames in order to cross-validate the procedure, finding that a background model with the first six principal components reproduced background counts with an rms error of 13 counts (≈0.1% of the average background level) when the fit excluded a mask representative of those for science frames. Although the final starlight subtraction is equally capable of subtracting the sky background, we found separating these steps into two stages improved our ability to create aligned cutouts and simplified later stages of the pipeline. The error in the background reconstruction of 13 counts amounts to 7 × 10⁻⁸ in contrast units for a star of Sirius' brightness, which we judge to be negligible.

In half-second exposures with the coronagraph splitting the incoming beam in two, Sirius still saturated the PSF core completely, along with some of the diffraction rings.

In order to model any companion accurately, we needed to infer the missing flux by fitting a model to the unsaturated portions of the Sirius PSF in each of the two gvAPP-180 halves. Shorter-integration Sirius data taken for calibration were also saturated, so we used a radial profile from stacked o Gruis observations with the same filter and coronagraphconfiguration to serve as our unsaturated reference. The o Gruis observations were taken the same night (2015 November 29 UT) under similar conditions in terms of PWV and with the MagAO system correcting 300 modes. The reported rms residual wave front error was approximately 100 nm, for a Strehl ratio S ≈ 97% before noncommon-path errors. In an ideal system, there would be no flux difference between the two gvAPP-180 halves, but Otten et al. (2017) noted that one half of the split PSF was brighter than the other by several percent, and we also observed this in our data. As the PSFs are split by their left- or right-circular polarization, this discrepancy indicates a small amount of circular polarization is introduced by the instrument itself. As we observed the effect on multiple targets, we concluded the origin is instrumental, rather than any astrophysical polarization signal.

We divided each frame into two cutouts centered on the complementary PSFs, which we aligned to subpixel precision. We established alignment to a known center of rotation by simulating a gvAPP-180 PSF without wave front errors, using said PSF to align an unsaturated calibration target PSF (in this case, the o Gruis template), and using this second, more realistic, PSF to align the Sirius frame cutouts. The result was two data cubes with dark-hole regions on opposite sides of the core, aligned to subpixel precision. We computed a radial profile of the template PSF, normalized to unit flux, and each frame's science PSFs. At each frame, we fit the scale factor that matched the radial profiles in the unsaturated region best, giving us the total Sirius flux. A comparison of the profiles in the top and bottom PSFs is given in Figure 2. The effect of the "lost" flux is illustrated by Figure 3.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Once in possession of two aligned data cubes, the amplitude difference between complementary gvAPP-180 PSFs was divided out using the per-frame scale factors.

The brightness of Sirius required short integration times, which in turn led to a large number of frames containing distinct realizations of the system PSF. The KLIP algorithm of Soummer et al. (2012) leverages the difference between the number of realizations n_obs and the number of random variables (i.e., pixels) p, under the assumption that n_obs ≪ p, to reduce the size of the matrix to be eigen decomposed to obtain the Karhunen–Lòeve basis. With n_obs ∼ 10,000 and p ∼ 7000 for our analysis, this assumption no longer held. In fact, when n_obs ≃ p, the cost of forming the covariance matrix outweighed any time savings compared to a direct computation of the singular value decomposition (SVD) of the data matrix (Long & Males 2021).

A common technique, adopted by Vigan et al. (2015) and many others, is to bin the resulting data temporally to reduce the dimensions of the problem. However, Males et al. (2021) found that we should expect speckle lifetimes of roughly 20–100 ms on bright stars, well under even our 500 ms Sirius exposures. As such, we expected that capturing the information from every frame within an interval rather than temporally binning will result in improved starlight-subtraction performance. This is, in fact, what we found for locations near the host star, as discussed in Section 3.6.

Initial attempts to apply the SVD downdate algorithm described in Long & Males (2021) to accelerate the reduction of the Sirius vAPP data proved insufficient. KLIP itself is vulnerable to self-subtraction in ADI data, as a planet PSF rotating through the field has enough overlap with itself from frame to frame that the resulting starlight-subtraction process is able to remove or severely attenuate it. This is usually mitigated by use of a PSF reference star or an angular exclusion criterion, removing frames temporally adjacent to the one under analysis. However, applying an angular exclusion criterion undid some of the speed gains of the SVD downdate algorithm, as the inner SVD had to grow by a row and a column for each frame to be removed. In such a high frame rate data set, when probing locations of interest near the IWA large numbers of frames had to be excluded.

For this reason, we developed the PCA-Temporal (PCAT) algorithm detailed in Section 3.5. By operating on pixel time series rather than frames, it is straightforward to mask a planet's path from the reference sample. This in turn means a single decomposition is required for all pixels (and therefore frames) to starlight subtract.

After normalizing the PSF cubes with the profile fitting procedure in Section 3.2, we conducted "fake planet" signal injection–recovery tests with the o Gruis template PSF as a realistic proxy for the unsaturated Sirius PSF. The paired PSFs of the gvAPP-180 required injecting the appropriate companion PSF template into each half of the data, as companions could rotate out of one dark hole and into another over the course of an observation.

To perform the injection–recovery test, a final focal plane location (ρ, θ) was selected (where ρ is the separation from the center of the stellar PSF and θ is the position angle (PA) in degrees E of N). For each frame and each half, the template PSF was scaled by amplitude A and translated to (ρ, θ + ϕ_i ) in the focal plane (where ϕ is the negative of the derotation angle that places north up and east left for frame i) before adding it to the original image. To mitigate the influence of the injected planet signals on each other and the potential biasing of the estimated contrast, we only injected one companion at a time to measure the contrast limit in a location.

The starlight in each frame was subtracted following the algorithm in Section 3.5, and the residuals were derotated to place north up and east left, and combined as the pixel-wise median of the stack of frames.

The signal and noise samples were computed from the final image by simple aperture photometry in λ/D diameter apertures, excluding apertures on either side of the measurement location (i.e., known injected planet). The signal-to-noise ratio (S/N) was then computed following the corrected equation for small-sample statistics from Mawet et al. (2014).

Each focal plane pixel can be treated as a time series that can be represented in a basis of eigen time series. This technique is employed by the Temporal Reference Analysis of Planets (TRAP; Samland et al. 2021), which also includes a constant term and a model light curve in the design matrix to be fit. However, the unique constraints of the paired vAPP PSF required development of a bespoke pipeline. By excluding pixels at the same separation as our injected companion from the basis determination, we avoid self-subtraction from overlapping companion flux. In essence, we have transposed the problem defined by Soummer et al. (2012) from eigenimages that describe the column space of a data matrix with one column per frame, to eigen time series that live in the row space as shown in Figure 4. In this way, we may use a single basis of eigen time series to starlight subtract every frame, eliminating even the (reduced) per-frame cost of the SVD downdate. Furthermore, by excluding an annulus of pixels from our reference set, we may reuse a single decomposition for multiple companion PA values at a given separation, as the reference sets are identical. We found that a more selective mask covering only the notional companion's track, while better for including information on diametrically opposite speckles, also led to overfitting.

Zoom In Zoom Out Reset image size

To begin, a mask is optionally applied to select pixels from an N-frames sequence of images, and each frame's p selected pixels are unwrapped into vector x _i , combined as the columns of the data matrix X ''_p×N . Let $\tilde{{\boldsymbol{x}}}$ be a vector whose entries correspond to the median value of each row in matrix X ''. The median-subtracted data are then ${{\boldsymbol{X}}}_{p\times N}^{{\prime} }={{\boldsymbol{X}}}_{p\times N}^{{\prime\prime} }-{\tilde{{\boldsymbol{x}}}}_{p\times 1}{{\bf{1}}}_{1\times N}$.

To reduce the effect of brighter pixels on the decomposition, we divide each row j of ${\boldsymbol{X}}^{\prime}$ by a whitening factor given by the standard deviation σ_j of that row. Let matrix W be a diagonal matrix ${\boldsymbol{W}}=\mathrm{diag}({\sigma }_{0}^{-1},...,{\sigma }_{p}^{-1})$. The preprocessed observation matrix X is then given by

$\begin{eqnarray*}&&{{\boldsymbol{X}}}_{p\times N}={\boldsymbol{W}}({\boldsymbol{X}}{{\prime\prime} }_{p\times N}-{\tilde{{\boldsymbol{x}}}}_{p\times 1}{{\bf{1}}}_{1\times N}).\end{eqnarray*}$

3.5.1. Special Considerations for the gvAPP-180 Data

3.5.2. Reference Time-series Selection

3.5.3. Eigen Time-series Computation

Having selected the reference pixel time-series vectors to form R , the algorithm proceeds with an SVD to obtain the left and right singular vectors as well as their associated singular values. We use SVD to mean the "economy" SVD that omits the decomposition of the null space of the matrix. The SVD allows us to decompose R into

$\begin{eqnarray*}&&{\boldsymbol{R}}={{\boldsymbol{U}}}_{(p-b)\times K}{{\boldsymbol{\Sigma }}}_{K\times K}{\left({{\boldsymbol{V}}}_{N\times K}\right)}^{T},\end{eqnarray*}$

where $K=\min ((p-b),N)$. To prevent overfitting and improve starlight subtraction, we retain the k < K column vectors in the decomposition corresponding to the greatest singular values. (If k ≪ K, this may allow the use of faster algorithms that find the partial SVD.) The decomposition is now only approximate, and given by

$\begin{eqnarray*}&&{\boldsymbol{R}}\approx \tilde{{\boldsymbol{R}}}={\tilde{{\boldsymbol{U}}}}_{(p-b)\times k}{\tilde{{\boldsymbol{\Sigma }}}}_{k\times k}{\left({\tilde{{\boldsymbol{V}}}}_{N\times k}\right)}^{T}.\end{eqnarray*}$

The question of the optimal choice of k is addressed in Section 3.6. The matrix $\tilde{{\boldsymbol{V}}}$ defined above gives us the eigen time series with which we will model each pixel's evolution.

3.5.4. Subtraction of the Estimated Starlight

Up until now, we have referred to observation column vectors ${{\boldsymbol{x}}}_{i}\in {{\mathbb{R}}}^{p}$. Now we operate on the rows of X , which we will call ${{\boldsymbol{y}}}_{j}\in {{\mathbb{R}}}^{N}$. For every pixel time series we use from the original data cube, we proceed to project it into the basis of eigen time series $\tilde{{\boldsymbol{V}}}$ and reconstruct it in terms of the first k eigen time series

$\begin{eqnarray*}&&{\tilde{{\boldsymbol{y}}}}_{j}\approx \tilde{{\boldsymbol{V}}}{(({\tilde{{\boldsymbol{V}}}}^{T})}_{k\times N}\,{{\boldsymbol{y}}}_{j}).\end{eqnarray*}$

The residual noise plus companion signal (if any) is a vector ${{\boldsymbol{s}}}_{j}\in {{\mathbb{R}}}^{N}$ given by ${{\boldsymbol{s}}}_{j}={{\boldsymbol{y}}}_{j}-{\tilde{{\boldsymbol{y}}}}_{j}$.

By repeating this procedure for all p rows of X , we obtain the starlight-subtracted signal S whose rows are s ₁, s ₂, ..., s _p

$\begin{eqnarray*}&&{{\boldsymbol{S}}}_{p\times N}={\left({{\boldsymbol{X}}}^{T}-\tilde{{\boldsymbol{V}}}({\tilde{{\boldsymbol{V}}}}^{T}{{\boldsymbol{X}}}^{T})\right)}^{T}.\end{eqnarray*}$

To reverse the whitening transformation, we use W ⁻¹ and obtain the final starlight-subtracted data D

$\begin{eqnarray*}&&{\boldsymbol{D}}={{\boldsymbol{W}}}_{p\times p}^{-1}{{\boldsymbol{S}}}_{p\times N}.\end{eqnarray*}$

The entries of the columns of D may then be mapped back into pixel locations in the original data cube to construct a cube of residuals and further postprocessed (e.g., by derotating and stacking images for ADI).

Due to the unique structure of the gvAPP-180 data, the optimal number of modes to subtract k varied with both separation and PA. Additionally, when calibrating the S/N = 5 contrast floor, we observed that injecting and recovering a signal of amplitude A producing an S/N of S ≫ 5 would predict a larger value A_S/N=5 = (5A)/S for the minimum contrast detectable at S/N = 5 than an injection–recovery test closer to the true S/N = 5 amplitude. Or, to put it another way, calibrating an accurate minimum contrast value depends on injecting and recovering a signal close to the minimum amplitude that is detectable at a good S/N, making this an iterative process. Calibrating the S/N = 5 contrast floor therefore required tuning of two hyperparameters for every location probed: the injected signal amplitude and the number of PCAT modes. The width of the annular mask excluding pixels from the reference sample was fixed at b = 16 pixels (≈2λ/D).

The hyperparameter-tuning problem is an area of ongoing research in the machine-learning literature (see, e.g., Yang & Shami 2020 for a review), with multiple techniques to explore the parameter space efficiently at our disposal. Since each parameter evaluation is comparatively expensive, we chose Bayesian optimization (Snoek et al. 2012), descended from the "kriging" technique (also known as Gaussian process regression) initially developed to predict the profitability of mining gold deposits from a small number of bore-hole samples (Krige 1951).

The optimization process attempts to learn a relationship between the number of PCAT modes, k, injected signal amplitude, A, and the resulting recovered S/N S at some location. The objective function we chose to maximize is

$\begin{eqnarray*}&&f(A,S)=-S{\mathrm{log}}_{10}A.\end{eqnarray*}$

This depends only on the recoverability of the injected signal, and does not attempt to optimize the S/N from reducing the no-injection data set. A separate forthcoming publication (J. D. Long et al. 2023, in preparation) will discuss the particulars of our approach to hyperparameter optimization and the distributed computing infrastructure developed for it. The optimization process was run for 200 iterations of selecting parameters, injecting a signal, reducing the data with said parameters, and measuring the recovered signal. Figure 5 shows the rapid convergence of the optimizer, with only marginal contrast gains after 50 iterations.

Zoom In Zoom Out Reset image size

In recognition of the difference in noise statistics between points in the speckle-dominated region close to the star, and background-limited points further out, we did compute contrast limits with varying amounts of coadding (shown in Figure 6). Retaining every frame in full temporal sampling produced better contrast limits at small separations, while performance at larger separations benefited somewhat from coadding temporally adjacent frames.

Zoom In Zoom Out Reset image size

We validated the performance of PCAT by comparison with our implementation of KLIP, with the same optimization scheme. We optimized the 10:1 coadded reduction with KLIP following the same procedure as PCAT, and show the results in Figure 7.

Zoom In Zoom Out Reset image size

To compute the final contrast surface, we collected all combinations of parameters evaluated for a particular focal plane location where the injected signal was recovered with at least an S/N ≥ 8. This cutoff was chosen arbitrarily to limit the analysis to points with relatively trustworthy estimates of the S/N = 5 contrast level A_{S/N = 5}=5A/8. The parameter combination with the smallest A_S/N=5 was saved, and those values define the surface in Figure 8. The k modes, injected amplitude A, and temporal downsampling factor n were allowed to vary, though we only evaluated fixed steps of n ∈ {1, 10, 50} due to computing resource constraints. The spatial distributions of these parameters are shown in Appendix.

Zoom In Zoom Out Reset image size

Predicting the properties of planets that one would be sensitive to depends on choosing a model for planet evolution and spectra. We adopt the latest "Bobcat" models in the "Sonora" series of Marley et al. (2021), which span a range of 0.5 < M/M_J < 84 in planet mass and 1 Myr < t < 15 Gyr in age. They model planets and substellar objects down to T_eff = 200 K with associated evolutionary tracks, high-resolution simulated spectra, and different metallicities. We perform our own synthetic photometry for this analysis using the published spectra.

The Bobcat models are published for [M/H] ∈ {−0.5, 0.0, +0.5}. Bond et al. (2017) models the evolution of Sirius A, reporting a slight metal deficiency relative to solar abundances and that an [Fe/H] of −0.13 or −0.07 reproduces their observations. We adopt the [M/H] = 0.0 Bobcat models as a conservative option, because lower metallicity models have higher [3.95] fluxes. As an example, an absolute [3.95] mag of 14 corresponds to masses of 7.5, 8.5, and 8.8 M_J in the [M/H] −0.5, 0.0, and +0.5 model suites, respectively.

4.1.1. Incorporating the Effects of Irradiation from Sirius A

Instellation from Sirius A would be a major contributor to the effective temperature of a giant planet at a small separation from the star. Males et al. (2014) showed that the mass–separation–absolute-magnitude relationship depends only weakly on mass at small separations, effectively making lower-mass planets brighter and thus more accessible to current-generation instruments. To account for this in our mass limits, we computed the equilibrium temperature

$\begin{eqnarray}&&{T}_{\mathrm{eq}}={T}_{\mathrm{Sirius}}\sqrt{\frac{{R}_{\mathrm{Sirius}}}{2a}}{\left(1-{A}_{B}\right)}^{1/4},\end{eqnarray} \tag{ 1 }$

where A_B is the Bond albedo and a is the separation between star and planet. Marley et al. (1999) found that cloud-free planets with incident flux from an A star have a mass- and temperature-dependent albedo ranging from 0.39 ≤ A_B ≤ 0.43 in their simulations. That range extends to 0.39 < A_B < 0.93 when considering different cloud models. For this analysis, we adopt A_B = 0.5 for consistency with Jupiter (Li et al. 2018).

The giant planets we might expect to find will glow with the heat of their formation and cool slowly over time to their equilibrium temperatures. We approximate the contribution to the equilibrium temperature T_eq from incident light by summing the planet luminosity

$\begin{eqnarray}&&{L}_{\mathrm{eff}}={L}_{\mathrm{evol}}+{L}_{\mathrm{eq}},\end{eqnarray} \tag{ 2 }$

where

$\begin{eqnarray}&&{L}_{\mathrm{evol}}={\sigma }_{\mathrm{sb}}4\pi {R}_{\mathrm{planet}}^{2}{T}_{\mathrm{evol}}^{4},\end{eqnarray} \tag{ 3 }$

$\begin{eqnarray}&&{L}_{\mathrm{eq}}={\sigma }_{\mathrm{sb}}4\pi {R}_{\mathrm{planet}}^{2}{T}_{\mathrm{eq}}^{4},\end{eqnarray} \tag{ 4 }$

and T_evol is the evolutionary temperature from the Bobcat isochrones. Substituting the Stefan–Boltzmann law for L_eff gives us the following relationship for T_eff of the irradiated planet

$\begin{eqnarray}&&{T}_{\mathrm{eff}}^{4}={T}_{\mathrm{evol}}^{4}+{T}_{\mathrm{eq}}^{4}.\end{eqnarray} \tag{ 5 }$

This modified relationship in the specific case of Sirius and the [3.95] filter is shown in Figure 12. We do not claim to have captured any of the subtleties of the evolution of an irradiated planet, only to have chosen an appropriate T_eff subject to what we know about the minimum equilibrium temperature for an assumed separation and albedo value. Because of the way the Bobcat models were initialized, models with very high T_eff for very low masses were not available. Therefore, the effective minimum modeled mass for our analysis is 1 M_J.

Zoom In Zoom Out Reset image size

4.1.2. Magnitude-to-mass Relationship

To obtain the relationship between the modeled masses and their fluxes in the [3.95] filter, we assumed a system age of 242 Myr based on the determination of Sirius A's age in Bond et al. (2017). For every modeled companion mass, T_eff and surface gravity were computed from the evolution tracks subject to the modification described above. We then interpolated a spectrum from the Bobcat grid, and computed a [3.95] absolute Vega magnitude for such an object.

To compute the synthetic photometry, we adopted the latest Vega model from HST CALSPEC ¹¹ (Bohlin et al. 2014) and an atmosphere with 2 mm of PWV and $\sec (z)=1.15$. The modeled atmosphere retrieved from ESO SkyCalc (Noll et al. 2012; Jones et al. 2013) represents the La Silla observatory, quite close to Las Campanas where these data were taken, and its slightly lower altitude means its extinction is slightly pessimistic compared to the actual situation.

With a [3.95] mag for each mass computed, we inverted the relationship. This relationship is shown for a selection of T_eq values in Figure 13. Where the relationship is double valued for masses around the minimum mass for deuterium burning, we made it monotonic by taking the minimum mass value for that magnitude, resulting in the apparent discontinuity around m_[3.95] = 11.5.

Zoom In Zoom Out Reset image size

4.1.3. Mass Limits

The minimum mass we would detect at S/N = 5 is shown as a function of separation in Figure 14. The overlaid scatter points and shaded region illustrate the variation in that mass limit as a function of angle at that separation, which is a consequence of the rotation of the dark-hole region, changes in conditions over the course of the observation, etc.

Zoom In Zoom Out Reset image size

To visualize the two-dimensional variation in the mass limit, when considering a point at angular coordinates (ρ, θ) in the focal plane, $a={d}_{\mathrm{Sirius}}\tan \rho$ is used as the separation for the purpose of computing a modified magnitude-to-mass relationship. The minimum detectable masses at each position are shown in Figure 15, adopting a Bond albedo of A_B = 0.5.

Zoom In Zoom Out Reset image size

While this is sufficient to exclude a planet at (ρ, θ) of that mass (or greater) and that albedo at separation a, it cannot necessarily exclude a planet of that mass but at a greater separation that happens to fall on (ρ, θ) in the focal plane when viewed from Earth. A planet at that projected separation, but far enough away to experience negligible irradiation, is equivalent to a planet with an albedo of A_B = 1.0. The mass limits under this assumption are given in Figure 16.

Zoom In Zoom Out Reset image size

Because this relationship depends on the projected separation (for the contrast limit) and the true separation (for the T_eq), it is more convenient and informative to fold these variables into the completeness analysis.

Our sensitivity is a function of the contrast floor (itself a function of position in the focal plane), companion mass, semimajor axis, and distance from the host star at the moment of observation. To capture this in our completeness analysis, we drew 10⁵ random orbital elements for a grid of masses and semimajor axis values. The eccentricities are drawn from a linearly descending prior given in Nielsen et al. (2019).

Projecting the position of the notional companion into the plane of the sky gives us the separation and PA we would observe, while the orbit determines the actual separation and hence the irradiation-affected effective temperature. The assumed age and metallicity, plus the companion mass and the projected and actual separations, then allow us to compute a [3.95] mag for the companion. In Figure 17, the fraction of the trial orbits that results in a [3.95] mag brighter than our S/N = 5 detection level gives the completeness at each mass and semimajor axis point.

Zoom In Zoom Out Reset image size

Using the S/N = 5 contrast versus separation curves of each of the studies, as well as their dates of observation, we were able to assess whether a randomly drawn orbit would have been detected in at least one of the studies. We converted the other studies' contrast limits into mass limits following the same procedure as our own, assuming radial symmetry. For each point in (mass, semimajor axis) space, we drew 10⁵ sets of orbital elements. Each was converted into a set of positions along the orbit using the observation dates of the studies.

We used the resulting true separations and projected separations to generate synthetic photometry in the bandpasses corresponding to each published contrast curve. As before, irradiation from the primary was taken into account. If any of the studies would have detected the planet at the corresponding position for their observation date, we counted it as detectable in the final combination of the completeness grid. The overall completeness from all four studies as a function of mass and semimajor axis is shown in Figure 18.

Zoom In Zoom Out Reset image size