NEW SPATIALLY RESOLVED OBSERVATIONS OF THE T Cha TRANSITION DISK AND CONSTRAINTS ON THE PREVIOUSLY CLAIMED SUBSTELLAR COMPANION

Since their discovery, transition disks have been regarded as natural laboratories for the study of protoplanetary disk evolution and perhaps planet formation. These objects' spectral energy distributions (SEDs) lack near- to mid-infrared emission, yet display a far-infrared excess (Strom et al. 1989). Initial studies attributed these SED features to a lack of warm dust at inner, AU-scale radii, suggesting that they were "in transition" from protoplanetary disks with excess throughout the infrared, to debris disks with only very weak far-infrared excess (e.g., Lin & Papaloizou 1986, 1993, p. 749; Bryden et al. 1999; Calvet et al. 2002). More recent modeling of Spitzer spectra (e.g., Calvet et al. 2005; Brown et al. 2007; Espaillat et al. 2007a, 2007b; Merin et al. 2010) also associated the mid-infrared deficits with disk cavities or gaps on AU scales. Followup submillimeter imaging (e.g., Brown et al. 2009; Hughes et al. 2009; Andrews et al. 2011) has directly confirmed the presence of these features.

Studies have identified several processes that could play a role in the formation of gaps and cavities, including photo-evaporative winds, grain growth, and dynamical interactions with companions. While photoevaporative winds would clear out only the gas and small dust in the inner disk (Clarke et al. 2001; Alexander et al. 2006), accounting for radial drift of solids can lead to dissipation of dust at small radii as well (Alexander & Armitage 2007). Furthermore, X-ray winds may drive disk depletion at a faster rate than UV winds, suggesting that X-ray photoevaporation could clear our inner disk radii more efficiently (e.g., Ercolano et al. 2008; Drake et al. 2009; Owen et al. 2010). However, the cavity sizes and mass loss rates observed in transition disks are too large compared to their X-ray luminosities to be consistent with clearing by photoevaporation (Andrews et al. 2011; Owen et al. 2011). Furthermore, measurements of outer disk masses are too large compared to results of simulations that cause disks to go through an "inner hole" phase (Alexander et al. 2006). While photoevaporation could explain some inner clearings, it alone cannot be responsible for transition disk structure.

Rather than clearing away disk material to lower the infrared emission, grain growth decreases its emissivity (Draine 2006; D'Alessio et al. 2006). Growing grains to mm/cm sizes can create SED deficits comparable to those observed in transition disks (e.g., Dullemond & Dominik 2005; Tanaka et al. 2005; Birnstiel et al. 2012). However, disk evolution simulations by Birnstiel et al. (2012) failed to generate the particle size distributions required to produce mm wavelength cavities. Additionally, soon after growing from ∼1 μm size to mm size, direct collisions between silicate particles could become destructive (e.g., Windmark et al. 2012), and the resulting smaller particles could then produce emission to fill in the SED deficit. This suggests that, while grain growth must impact disk evolution at some level, this process alone cannot shape transition disk cavities.

Dynamical interactions with companions are the best explanation to date for forming disk gaps. Models have demonstrated that stellar mass companions can open cavities in disks (e.g., Artymowicz & Lubow 1994; Pichardo et al. 2008), and some observed disk gaps, such as those in CoKu Tau 4 (Ireland & Kraus 2008), HD 98800 (Furlan et al. 2007), and Hen 3-600 (Uchida et al. 2004), have been associated with stellar mass binaries. However, high resolution imaging has ruled out companions with masses higher than ∼20–30 M_Jup for approximately half of the known transition disks (Kraus et al. 2011; Evans et al. 2012). This leaves the exciting possibility that planetary mass companions are clearing out cavities and gaps, accreting material that would have otherwise fallen onto the star (e.g., Najita et al. 2007). Simulations have shown that tidal interactions with a planetary mass companion can indeed open gaps in disks (e.g., Lin & Papaloizou 1986; Bryden et al. 1999; Crida et al. 2006).

Here we discuss one transition disk thought to be shaped by a substellar mass companion, T Chamaeleontis (T Cha). T Cha is a G8 type, 1.5 M_☉ star, first categorized as a weak-line T Tauri star due to its low Hα equivalent width of <10 Å (Alcala et al. 1993). This suggested it had entered the final stages of accretion. The classification conflicted with T Cha's infrared excess, thought to result from its outer disk; this spectral feature would place it as a classical T Tauri star. Later observations showed that it is in fact a transition disk object, perhaps in the intermediate stages between a protoplanetary disk and a disk-free planetary system. Brown et al. (2007) found that T Cha's SED could be reproduced by a disk with a gap between 0.2 and 15 AU, and SED modeling by Olofsson et al. (2011) supported this, with a best fit gap extending from 0.17 to 7.5 AU.

Imaging observations suggested the presence of a companion of L' contrast 5.1 mag at a separation of 62 mas—6.7 AU at a distance of 108 pc—(Huelamo et al. 2011), within T Cha's disk gap. However, Olofsson et al. (2013) showed that the detected signal could be modeled almost equally well by an asymmetry caused by forward scattering from the upper layers of the outer disk. Observing orbital motion of the companion candidate would distinguish between these two scenarios. To this end, we acquired new observations of T Cha with the Magellan AO (MagAO) system. We also present a re-analysis of the original discovery data from VLT/NaCo, as well as a new analysis of previously unpublished NaCo data from the Very Large Telescope (VLT) archive.

Non-redundant masking (NRM) transforms a filled aperture into a sparse interferometric array using a pupil-plane mask. While blocking the majority of the light, this provides much better knowledge of the point-spread function than a conventional telescope. A resulting image (called an interferogram) then shows the interference fringes formed by the mask, and subsequent image reconstruction or model fitting relies on quantities calculated from its Fourier transform. Since the mask is non-redundant (no two baselines have the same length and orientation), each baseline has unique (u, v) coordinates. The symmetry of the Fourier transform means that identical information for a baseline can be found in two points in (u, v) space—at (u, v) and (− u, −v). The finite size of the mask holes, as well as the width of the filter bandpass, causes this information to spread out. This means that the Fourier transform will have several distinct "splodges," two coming from each mask baseline.

Using the Fourier transform at the locations of these splodges, we find the complex visibility for each baseline, which has the form Ae^iϕ, where A is the amplitude and ϕ the phase. Since atmospheric and instrumental effects corrupt the complex visibilities, we calculate two other quantities, squared visibilities and closure phases. Squared visibilities measure the power in the Fourier transform as a function of baseline length, and closure phases are sums of phases around three baselines that form a triangle. Closure phases eliminate atmospheric and instrumental phase offsets that corrupt measurements taken using single baselines. These obey the relation:

$$\begin{equation} \Phi _{cp} = \phi \left(u_1,v_1\right) + \phi \left(u_2,v_2\right) + \phi \left(u_3,v_3\right), \end{equation} \tag{ 1 }$$

where u_i and v_i are the sampling coordinates of the ith baseline in the Fourier plane. An N-hole mask will provide ${N \choose 2}$ baselines and visibilities, and ${N \choose 3}$ closure phases, ${N-1 \choose 2}$ of which are independent.

3.1. New Magellan/MagAO/Clio2 Data

We observed T Cha and two unresolved calibrators, HD 101251 and HD 102260 using the 6.5 m Clay telescope, MagAO adaptive optics system (Close et al. 2013; Morzinski et al. 2014), and Clio2 science camera (Freed et al. 2004; Sivanadam et al. 2006) on 2013 April 5. A six-hole non-redundant mask was mounted in a pupil plane filter wheel in Clio2. Table 1 lists the parameters of the Clio2 mask. We used two calibrators to lessen the probability that detected signals were being injected by a resolved or binary calibrator. Our exposure times for T Cha, HD 101251, and HD 102260 were 1.0 s, 0.9 s, and 1.3 s, respectively. We took four to six 50-frame data cubes (called "visits" in Table 2) before switching objects according to the pattern target-cal1-target-cal2. We acquired six visits for T Cha, two visits for HD 101251, and three visits for HD 102260. These resulted in 28.3 minutes, 9 minutes, and 15.16 minutes of total integration, respectively. We observed in L', λ_c = 3.78 μm, and the total change in sky rotation angle was 47°, shown in Figure 1.

Zoom In Zoom Out Reset image size

Table 1. Mask Parameters

Mask	Nholes	Nb	Ncp	Nkp	Baseline Range	Throughput
					(m)	(%)
NaCo	7	21	35	15	1.77–6.45	16%
MagAO/Clio2	6	15	20	10	1.68–5.02	11%

Download table as:  ASCIITypeset image

Table 2. Summary of Observations

Target	Right Ascension	Declination	tint	Nframesa	Nvisitsb	Total Time	Seeing	τ0
	(hh mm ss.sss)	(dd mm ss.sss)	(s)			(minutes)	(arcsec)	(ms)
L' Observations
2010 Mar 17: VLT/NaCo
T Cha	11 57 13.550	−79 21 31.537	0.4	800	9	48	0.6	8
HD 102260	11 45 13.822	−78 36 58.633	0.4	800	10	53.3
2011 Mar 14: VLT/NaCo
T Cha	11 57 13.550	−79 21 31.537	0.4	800	7	37.3	2.0	1.5
HD 102260	11 45 13.822	−78 36 58.633	0.4	800	9	48
2012 Mar 8: VLT/NaCo
T Cha	11 57 13.550	−79 21 31.537	0.4	800	5	26.7	1.0	6
HD 101251	11 37 49.220	−79 14 31.604	0.4	800	2	10.7
HD 102260	11 45 13.822	−78 36 58.633	0.4	800	2	10.7
2013 Mar 25: VLT/NaCo
T Cha	11 57 13.550	−79 21 31.537	0.4	800	2	10.7	0.75	7
T Cha	11 57 13.550	−79 21 31.537	0.3	1057	7	37.0
T Cha	11 57 13.550	−79 21 31.537	0.2	1405	1	4.7
HD 102260	11 45 13.822	−78 36 58.633	0.4	800	1	5.3
HD 102260	11 45 13.822	−78 36 58.633	0.3	1057	4	21.1
HD 102260	11 45 13.822	−78 36 58.633	0.5	1057	1	8.8
HD 101251	11 37 49.220	−79 14 31.604	0.4	800	1	5.3
HD 101251	11 37 49.220	−79 14 31.604	0.3	1058	2	10.6
HD 101251	11 37 49.220	−79 14 31.604	0.15	1687	1	4.2
2013 Apr 5: Magellan/MagAO/Clio2
T Cha	11 57 13.550	−79 21 31.537	1.0	300	5	28.3	0.6	N/A
T Cha	11 57 13.550	−79 21 31.537	1.0	200	1	28.3
HD 101251	11 37 49.220	−79 14 31.604	0.9	300	2	9
HD 102260	11 45 13.822	−78 36 58.633	1.3	300	1	15.16
HD 102260	11 45 13.822	−78 36 58.633	1.3	200	2	15.16
Ks Observations
2010 Jul 1: VLT/NaCo
T Cha	11 57 13.550	−79 21 31.537	0.5	800	3	20	1	4
HD 101251	11 37 49.220	−79 14 31.604	0.5	800	1	6.7
HD 102260	11 45 13.822	−78 36 58.633	0.5	800	1	6.7
2011 Mar 15: VLT/NaCo
T Cha	11 57 13.550	−79 21 31.537	2.0	200	6	40	0.75	3
HD 102260	11 45 13.822	−78 36 58.633	2.0	200	6	40
2013 Mar 26: VLT/NaCo
T Cha	11 57 13.550	−79 21 31.537	1.0	280	2	9.3	1	4
T Cha	11 57 13.550	−79 21 31.537	0.8	427	4	22.8
T Cha	11 57 13.550	−79 21 31.537	0.7	497	3	17.4
T Cha	11 57 13.550	−79 21 31.537	0.5	807	2	13.5
HD 102260	11 45 13.822	−78 36 58.633	1.0	350	1	5.8
HD 102260	11 45 13.822	−78 36 58.633	0.5	707	2	11.8
HD 102260	11 45 13.822	−78 36 58.633	0.4	707	1	4.7
HD 102260	11 45 13.822	−78 36 58.633	0.3	1057	1	5.3
HD 101251	11 37 49.220	−79 14 31.604	0.2	1407	3	14.1
HD 101251	11 37 49.220	−79 14 31.604	0.11	2782	1	5.1
H Observations
2013 Mar 27: VLT/NaCo
T Cha	11 57 13.550	−79 21 31.537	0.4	605	1	4.03	0.75	3.5
T Cha	11 57 13.550	−79 21 31.537	0.4	847	1	5.64
T Cha	11 57 13.550	−79 21 31.537	0.5	707	4	23.6
HD 102260	11 45 13.822	−78 36 58.633	0.15	2107	3	15.8
HD 101251	11 37 49.220	−79 14 31.604	0.11	3157	2	11.6
HD 101251	11 37 49.220	−79 14 31.604	0.11	1353	1	2.5

Notes. ^aNumber of frames in each visit ^bEach visit consists of all images taken before switching between target and calibrator.

Download table as:  ASCIITypeset image

3.2. Previously Published 2010 VLT/NaCo Data

3.3. Unpublished VLT/NaCo Data

We have developed a suite of software in Python to perform basic data reduction, calibration, and visibility and closure phase calculations. We first flat-field and bad-pixel correct all images. For a given set of two dithers, we perform sky subtraction for one position by taking the median of all images in the other dither position and subtracting the median from each image. We then apply a super-Gaussian window function to reduce the noise associated with the low-signal edges of the interferogram and spatially filter the data. A super-Gaussian has the form exp (− kx⁴); we choose k such that the half width at half max is λ/d_sub, where d_sub is the mask sub-aperture diameter (e.g., Bernat 2012). After windowing, we Fourier transform the data.

Next, we calculate squared visibilities. We first simulate an interferogram and resulting Fourier transform using the locations and sizes of the holes in the mask, along with the observation's wavelength. This provides us with the pixel locations of the splodges (see Section 2) in the Fourier transform. We then square the Fourier transform of the data, and sum all pixels within the splodges corresponding to each baseline, normalizing by the total power in the interferogram.

The typical uncertainty due to random errors for the visibilities ranges between 0.02 and 0.12 for all L' and 2013 Ks data sets. For comparison, a binary with separation 62 mas and contrast ΔL' = 5.1 mag produces a change in visibility between the shortest and longest baseline of approximately 0.035 for both the MagAO/Clio2 and NaCo masks. While some of the followup data sets' random visibility errors are smaller than this signal, visibilities' dependence on the AO system (e.g., Lacour et al. 2011; Kraus & Ireland 2012) renders them harder to calibrate. Differences in AO performance over the night, or between target and calibrator observations can introduce additional error. To investigate this, we divide the visibilities for each set of two adjacent calibrator scans; the calibrated visibilities for a point source should be equal to 1 for all baselines. We then take the scatter in these calibrated visibilities as an estimate for the systematic uncertainties in the target visibilities. This results in calibrated visibilities with scatters ranging between 0.04 and 0.11. Due to these large uncertainties, as in previous NRM studies (e.g., Huelamo et al. 2011; Kraus & Ireland 2012) we restrict our binary fitting to phases only, rather than including the squared visibilities.

For each triangle, we form the bispectrum by multiplying the complex values in the Fourier transform at three (u, v) coordinates. We then average over the individual frames, and take the phase of the average bispectrum to be our closure phase. Of the ${N \choose 3}$ closure phases, only ${N-1 \choose 2}$ are independent. For this reason, we perform fits on kernel phases, linearly independent combinations of closure phases (see Martinache 2010; Kraus & Ireland 2012; and Ireland 2013). We find our kernel phases in a way similar to Martinache (2010). Appendix C gives a detailed description of our projection method.

Figure 2 shows histograms of the kernel phases for individual images in the Ks and L' data sets, and can be taken to represent a comparison of the snapshot kernel phase errors for the different data sets. For a given dither, we calculate the ${N-1 \choose 2}$ mean kernel phases. We then calculate the kernel phases for every individual image in the dither and subtract the mean. We use the same total integration time, 39 s, to calculate the mean kernel phase. This process yields the distributions shown in Figure 2. The snapshot errors in the 2013 MagAO data (σ = 073) are much smaller than the VLT/NaCo data sets (σ = 2.85–438 for L' kernel phases.) We also compare the MagAO/Clio2 and NaCo data with different amounts of time-averaging. As long as these averages are performed within a single dithering sequence, the chosen interval does not change the noise levels significantly. Due to their high scatter and the resulting unreliable fits we show kernel phase histograms for the 2010 and 2011 Ks as well as the 2013 H band data in Appendix A.

Zoom In Zoom Out Reset image size

Random sources of noise associated with both AO performance and observing conditions cause snapshot kernel phase errors. For exposures much longer than the inverse of the AO system bandwidth, closure phase errors scale in the following way (Ireland 2013):

$$\begin{equation} \sigma _{cp} \sim \sigma _{\phi } \sqrt{\frac{1}{f_c T}}, \end{equation} \tag{ 2 }$$

where σ_ϕ is the phase noise on each sub-aperture in the closing triangle, f_c the cutoff frequency below which piston noise is white, and T the exposure time. The kernel phase errors will then be the projection of the closure phase errors. If we assume σ_ϕ and f_c are equal for all data sets, the ratio of two observations' closure phase errors scales with the ratio of their exposure times. Using 0.2–0.5 s for the NaCo data, and 1.0 s for the MagAO/Clio2 data, the ratio of the MagAO/Clio2 to NaCo closure phase errors should range from 0.44 to 0.71. The ratio of the MagAO/Clio2 and NaCo L' observed errors ranges between 0.27 and 0.18 depending on the NaCo L' data set. Hence, exposure time alone cannot explain the discrepancy between the MagAO/Clio2 and NaCo L' snapshot errors.

Other random sources of noise include photon, background, and read noise, which lead to closure phase errors (following Ireland 2013):

$$\begin{equation} \sigma _{cp} \sim \frac{N_h}{N_p V}\sqrt{1.5\left(N_p+N_b+n_p \sigma ^2_{ro}\right)}, \end{equation} \tag{ 3 }$$

where N_h is the number of holes in the mask, V the fringe visibility, N_p the total number of photons, N_b the number of background photons, n_p the number of pixels, and σ_ro the read noise. If we assume we are in the photon-noise regime, the closure phase error simplifies to

$$\begin{equation} \sigma _{cp} \sim \frac{N_h}{V}\sqrt{\frac{1.5}{N_p}}. \end{equation} \tag{ 4 }$$

The number of photons is N_p = F_target × A_tel × f_mask × T, where F_target is the photon flux from the source, A_tel the telescope collecting area, f_mask the fraction of light allowed through by the mask, and T the exposure time. The relevant values for the MagAO/Clio2 observations are A_tel ∼ (6.5 m)², f_mask = 0.11, T = 1s, and N_h = 6. For the NaCo observations, A_tel ∼ (8.2 m)², f_mask = 0.16, T = 0.2 − 0.4 s, and N_h = 7. Assuming that the fringe visibilities for the two data sets are approximately equal, taking the ratio of the MagAO/Clio2 to NaCo closure phase errors gives 0.58–0.92. Thus, the differences in telescope and observing parameters cannot fully explain the lower MagAO/Clio2 snapshot errors. Better AO performance by MagAO, which leads to lower values for σ_ϕ and higher values of f_c, could be one cause of this discrepancy. Additionally, the smaller holes (as evinced by the lower throughput in Table 1) of the Clio2 mask mean that the 2013 observations are less redundant than the NaCo observations. This could also reduce the snapshot errors for the MagAO/Clio2 data.

Due to atmospheric and instrumental systematics, the mean kernel phases themselves (subtracted off in Figure 2) can vary substantially throughout the observations. To take these effects into account, we subtract our unresolved calibrator kernel phases from our target kernel phases. We do this in several ways.

To apply a simple nearest-neighbor calibration, we first calculate the time between a given target scan and all calibrator scans (Δt). We then average all calibrator scans, weighting by Δt⁻¹⁰, and subtract the weighted-average calibrator from the target scan. To use information from all of the calibrator scans, rather than limiting ourselves to only the nearest-neighbor, we calculate an average calibrator, weighting the scans by Δt⁻¹. Finally, we apply a more optimized weighting, similar to LOCI (Lafreniere et al. 2007) techniques in direct imaging data reduction and the calibration strategy adopted in Kraus & Ireland (2012) and Ireland (2013).

For the LOCI-like calibration scheme, we find the linear combination of calibrator scans that minimizes the sum of the target's squared kernel phases. This is equivalent to minimizing the χ² for the null model. While this calibration scheme provides the lowest scatter and thus highest signal to noise, it can also subtract signal from the measurements. For this reason, it is often applied iteratively, minimizing the χ² for the null model initially and then minimizing the χ² for the best-fit model until the best-fit converges. We also LOCI calibrate without iteration, minimizing the χ² of the Δt⁻¹ model. Both LOCI schemes remove signal, and do not always give consistent results. Therefore, we focus on the fits to simpler, neighbor-like calibrations in the subsequent sections. The behavior of the LOCI calibration method will be detailed in a future paper.

The histograms for the L' and 2013 Ks data sets are shown in Figure 3. While the MagAO/Clio2 data have lower snapshot kernel phase errors than the NaCo L' observations (see Figure 2), the systematics in the three data sets are such that the scatter in the mean, calibrated kernel phases are quite similar. We speculate that the greater temporal spacing between target and calibrator observations in the MagAO/Clio2 observations could cause this. The average time between target and calibrator scans in the these data is 8.8 minutes. For the NaCo observations, the mean time between target and calibrator scans for the L' data ranges between 4.2 and 6.9 minutes.

Zoom In Zoom Out Reset image size

5.1. Consistency Checks

To search for companions in our data, we fit binary models to our kernel phases. Given the angular resolution of our observations, a binary can be approximated as two delta functions, the Fourier transform of which is an analytic function. Equation (5) gives the complex visibility for a binary with a companion separation s, position angle P.A. (measured E of N), and brightness ratio b (in units of the primary's brightness):

$$\begin{equation} V(u,v) = \frac{1}{\sqrt{2 \pi }}(1 + b e^{i\cdot s(u\cdot \sin ({\rm P.A.}) + v\cdot \cos ({\rm P.A.}))}). \end{equation} \tag{ 5 }$$

The phase measured for a binary by a baseline with coordinates (u, v), is then the angle of Equation (5):

$$\begin{equation} \phi (u,v) = \tan ^{-1}\left(\frac{{\rm Im}\, (V(u,v))}{{\rm Re}\, (V(u,v))}\right). \end{equation} \tag{ 6 }$$

To create a model set of kernel phases, we calculate the closure phases for each triangle using the locations and sizes of the mask sub-apertures, and sky rotation angles included in our observations. We then project the closure phases into kernel phases in the same way as we have done for the data. Due to the symmetry of the Fourier transform, each closure phase can correspond to one of two closing triangles. To keep the sign of our closure phases (and thus kernel phases) consistent between data and model, we sample the same closing triangles in both. Additionally, where necessary, we use observations of a known binary to calibrate the orientation of our detector on the sky. For the NaCo data, as in Huelamo et al. (2011), we use observations of the binary θ Ori C, taken in 2010 April.

We perform χ² fitting of the binary models to each individual data set using both a grid and nested sampling (Sivia & Skilling 2006). Our grid spans a range of parameters in binary position angle (P.A.), separation (s), and contrast in magnitudes (Δ). Δ can be related to b, the brightness ratio, by the following:

$$\begin{equation} \Delta = -2.5 \log _{10}(b). \end{equation} \tag{ 7 }$$

Table 3 lists the ranges and spacings for each parameter. For each set of parameters, we calculated model kernel phases and a χ² statistic. We used the best grid fit as an input for our nested sampling algorithm.

Nested sampling involves first filling the parameter space with a large number (in our case, 100) of points, and calculating a likelihood (exp (− χ²/2)) for each point. We then replace the lowest likelihood point with a random member of the remaining 99, and evolve it using a Markov chain to a higher likelihood region of the parameter space. We repeat this process until the 100 points satisfy a convergence criterion, in our case, the scatter in the ensemble must be a small fraction (∼0.1%) of the mean. We give one of the 100 nested sampling points the best grid fit as an initial value, and assign random values to the remaining 99. This is not necessary for fits in which there is one clear minimum, but it can help in preventing the nested sampling fit from converging to local minima.

We report parameter errors calculated from a χ² interval. After finding the minimum χ² using nested sampling, we then scale all χ² values for a grid of parameters so that the reduced χ² of the best-fit model is equal to 1. With the scaled set of χ² values, we find all grid points within Δχ² of 3.53 (Press et al. 1992) to place a 1σ error on the fit parameters.

Bootstrapping often gives the most conservative estimates of parameter errors. However, as noted in Press et al. (1992), the results of data fitting in Fourier space rely heavily on all grid points being present, and while the (u, v) coverage in the NaCo and MagAO data sets is good, it is by no means complete. For this reason, bootstrapped data sets do not fairly represent the noise in the data. We confirmed this by bootstrapping Gaussian noise sampled at the same (u, v) points as each data set; a given noise realization's best fit contrast ratio was much higher than the bootstrapped distribution would suggest.

For each data set, we fit simulated kernel phases to quantify our type I (false-positive) and type II (false-negative) errors. The contrast ratio and separation parameters in a binary model act much like the amplitude and frequency of a sine wave; they take on non-zero values when fit to noise. In order to determine the separations and contrast ratios that could be caused by the noise levels in these data sets, we simulated Gaussian kernel phases from distributions fit to the data. These are shown in Figure 3 for L' and 2013 Ks data, and in Appendix A for 2010 and 2011 Ks and H band data. We use (u, v) coverage and sky rotation identical to each observation. We fit binary models to 1000 of these noise realizations and create a probability distribution from the best fits. Since the best fit position angle for the noise simulations is uniformly distributed, we create two-dimensional confidence intervals from the best fit separations and contrasts. We then compare the best fit from our data to these confidence intervals. For example, if our best fit lies just outside the contour enclosing 95% of the best fits, there is a 5% chance that the fit is drawn from the distribution and thus a 5% chance that the fit resulted from noise alone.

A second way of quantifying type I errors is to calculate the F statistic, the best fit χ² divided by the null model χ², for each noise realization. Comparing the F statistic for a data set's best fit to a distribution of these simulated noise F statistics yields the probability that the best fit incorrectly rejects the null hypothesis. For example, if 95% of the simulated noise realizations have lower F statistics than a data set's best fit, there is a 95% probability that the data contain no real signal. This procedure is outlined in detail in Protassov & van Dyk (2002).

To estimate our type II errors—the likelihood that the best fit to our data would be a false negative—we simulated observations with noise plus the signal from a companion. We generate 1000 Gaussian noise realizations and fit a binary model to the simulated data. We take the fraction of fits where the input signal was not recovered to be our type II error. We also calculate F statistics from the noise + signal realizations, for comparison with both the noise F statistics and the data's best fit F statistics.

Figure 4 shows the noise simulations for all L' and 2013 Ks data sets with a Δt⁻¹ calibration. The scattered points show the results of 1000 fits to Gaussian noise comparable to the scatter in each data set, and the points with error bars show the best fit to each data set. The color scale shows a probability distribution interpolated from the results of the simulations, and the contour lines indicate 1σ, 2σ, and 3σ confidence intervals.

Zoom In Zoom Out Reset image size

Figure 5 shows the false alarm probabilities for the same data sets in Figure 4. In each panel, the black line is a histogram of all of the F statistics for 1000 Gaussian noise realizations. The red vertical lines indicate the F statistic for each data set's best fit; the intersection of these lines with the black histograms gives the probability that the best fit to the data resulted from noise. The green cumulative histogram shows the distribution of F statistics for the 1000 noise + signal realizations. The intersection of the red line with this distribution gives the fraction of noise + signal realizations with fits that look less significant (higher F statistics) than our best fit. For Figures 4 and 5, the noise properties of the Δt⁻¹⁰ calibration yielded comparable results (see Table 4).

Zoom In Zoom Out Reset image size

Table 4. Binary Fit Results

Data Set	P.A.	Separation	ΔL'	P(FA)a	P(FA)b	$\Delta _{\chi ^2,Hu}$	σHu, allowed
	(deg)	(mas)		(%)	(%)
2010 Mar 17: VLT/NaCo L'
Δt−1	83 $\mathbf {\pm ^{7}_{9}}$	88 $\mathbf {\pm ^{22}_{58}}$	5.5 $\mathbf {\pm ^{0.5}_{2.5}}$	11	< 0.1	2.63	1
Δt−10	82 $\pm ^{5}_{7}$	78 $\pm ^{22}_{48}$	5.3 $\pm ^{0.6}_{2.3}$	4	<0.1	2.56	1
2011 Mar 14: VLT/NaCo L'
Δt−1	92 $\mathbf {\pm ^{11}_{14}}$	87 $\mathbf {\pm ^{33}_{57}}$	5.3 $\mathbf {\pm ^{0.7}_{2.3}}$	25	18	2.54	1
Δt−10	92 $\pm ^{117}_{21}$	94 $\pm ^{356}_{64}$	5.4 $\pm ^{0.8}_{2.4}$	18	36	1.25	1
2012 Mar 8: VLT/NaCo L'
Δt−1	−40 $\mathbf {\pm ^{5}_{5}}$	32 $\mathbf {\pm ^{28}_{2}}$	3.4 $\mathbf {\pm ^{1.6}_{3.2}}$	1.3	< 0.1	23.10	> 4
Δt−10	−38 $\pm ^{5}_{5}$	29 $\pm ^{54}_{6}$	0.4 $\pm ^{5.3}_{0.3}$	4.3	<0.1	20.43	4
2013 Mar 25: VLT/NaCo L'
Δt−1	83 $\mathbf {\pm ^{3}_{1}}$	55 $\mathbf {\pm ^{25}_{25}}$	5.2 $\mathbf {\pm ^{0.7}_{1.9}}$	0.9	< 0.1	25.82	> 4
Δt−10	82 $\pm ^{4}_{4}$	25 $\pm ^{50}_{0.01}$	2.8 $\pm ^{3.1}_{0.2}$	3.4	<0.1	17.35	4
2013 Mar 26: VLT/NaCo Ks
Δt−1	74 $\mathbf {\pm ^{4}_{0.3}}$	42 $\mathbf {\pm ^{8}_{22}}$	5.2 $\mathbf {\pm ^{0.3}_{1.8}}$	1.3	< 0.1	21.47	4
Δt−10	77 $\pm ^{1}_{3}$	51 $\pm ^{9}_{11}$	5.4 $\pm ^{0.2}_{0.3}$	0.4	<0.1	10.64	2
2013 Apr 5: Magellan/MagAO/Clio2 L'
Δt−1	112 $\mathbf {\pm ^{176}_{99}}$	337 $\mathbf {\pm ^{153}_{147}}$	5.8 $\mathbf {\pm ^{0.6}_{0.4}}$	32	2	15.97	4
Δt−10	−131 $\pm ^{3}_{5}$	315 $\pm ^{25}_{26}$	5.8 $\pm ^{0.5}_{0.3}$	32	4.5	12.85	3

Notes. ^aUsing distribution of noise simulation best fits. ^bUsing distribution of noise simulation F statistics.

Download table as:  ASCIITypeset image

While using a χ² statistic to assign significance depends on the data error bars, which could be difficult to estimate properly, these simulations take into account only the scatter in the kernel phases. For this reason, in Section 8, we use the results of the noise simulations to estimate the significance of the best fit binary results.

Table 4 shows our results in chronological order. The first three columns list the results of binary fits to the kernel phases. Following the binary fits, the next columns list two false-alarm probabilities, the first calculated using the distribution of best fits to noise realizations and the second using the distribution of best fit F statistics (see Section 7). The last two columns list the Δχ² corresponding to the Huelamo et al. (2011) binary for that data set, alongside the corresponding confidence interval at which it is allowed. The bold values are binary fit results for the Δt⁻¹ calibrated data. We include the best fits from this calibration strategy in Section 9.

8.1. 2010 VLT/NaCo L' Data: Initial Detection

For these data, we first fit our neighbor-calibrated (weighting by Δt⁻¹⁰) kernel phases, since the calibration strategy for the closure phases in Huelamo et al. (2011) is most similar to this approach (see also Lacour et al. 2011). We find a best fit binary with position angle of $82\pm ^{5}_{7}$°, separation of $78\pm ^{22}_{48}$ mas, and $\Delta L^{\prime } = 5.3\pm ^{0.6}_{2.3}$ mag. This is comparable to the Huelamo et al. (2011) best fit—position angle of 78° ± 1°, separation of 62 ± 7 mas (6.7 AU at 108 pc), and ΔL' of 5.1 ± 0.2 mag. Our fit parameter errors are substantially larger than those quoted in Huelamo et al. (2011), which were derived using χ² intervals from a binary fit to closure phases (see also Lacour et al. 2011). Closure phases have correlated errors, which could bias the parameter errors derived using a simple χ² fit. Parameter errors also depend on whether the data error bars have been scaled such that the reduced χ² is equal to 1, which is appropriate if you assume the model to be correct, and that underestimated error bars are causing a high reduced χ². Our large separation and contrast error bars are consistent with the severe degeneracy between these two parameters. Figure 6 illustrates this degeneracy, with kernel phases plotted for separations 0.5 ×, 1 ×, and 1.5 × the separation of the Huelamo et al. (2011) binary, varying the contrast ratio. The three models are nearly indistinguishable. For this reason, we believe our fit errors to be a more realistic representation of the parameter constraints.

Zoom In Zoom Out Reset image size

After the nearest-neighbor fit, we then fit the Δt⁻¹ calibrated kernel phases, finding a best fit with position angle of $83\pm ^{7}_{9}$°, separation of $88\pm ^{22}_{58}$ mas, and $\Delta L^{\prime } = 5.5\pm ^{0.5}_{2.5}$ mag. The Δt⁻¹ calibration method resulted in kernel phases with lower scatter than the Δt⁻¹⁰ method—063 versus 066. Using a χ² interval, the Huelamo et al. (2011) binary model is within 1σ of the best fit for both calibration methods. Appendix B shows our Δt⁻¹ calibrated kernel phases with both the Huelamo et al. (2011) model and the best fit model from this work. The two models are nearly indistinguishable.

Figure 4 shows that there is a small probability that the binary detection resulted from a random noise fluctuation. The best fit falls on a contour that encloses 89% of the simulation results; there is an 11% chance that the 2010 fit was the result of noise. The 2010 L' best fit is then significant at roughly the 2σ level. The distribution of F statistics, shown in Figure 5, indicates that the false alarm probability for this data set is <0.1%, giving the fit ∼3σ significance.

8.2. VLT/NaCo 2011 L' Data

8.3. VLT/NaCo 2012 L' Data

8.4. VLT/NaCo 2013 L' Data

8.5. VLT/NaCo 2013 Ks Data

8.6. Magellan/MagAO/Clio2 2013 Data

The NaCo 2012 and MagAO/Clio2 2013 L' data sets have best fit binaries that are inconsistent with both the Huelamo et al. (2011) model and our best fit to the NaCo 2013 L' data. However, simulations show that there is a non-zero chance that these fits resulted from noise—1.3% for the NaCo and 32% for the MagAO/Clio2 observations. We simulated noise realizations (see Section 7) to estimate our type II errors for both of these data sets, using the NaCo 2013 L' signal as the input binary model. The chance that we would have missed the binary signal in these data is 15.8% for NaCo 2012 and 49.9% for MagAO/Clio2 2013.

Keeping these type II errors in mind, in the subsections concerning orbital motion and forward scattering we first assume that the tentative NaCo 2012 detection is reliable but that the MagAO/Clio2 non-detection is due to noise. We then discuss the results assuming that noise fluctuations led to a false detection in the NaCo 2012 data, while a signal compatible with the other NaCo data sets was actually present beneath the noise. We also consider the possibility that asymmetries caused by planet–disk interactions could have caused the observed kernel phases. In the last subsections, we discuss whether a chance alignment or systematic error could masquerade as a companion in the data sets where we detect a significant signal.

9.1. Companion Orbital Motion

Detecting orbital motion of a binary signal in multi-epoch data sets would confirm the claimed companion from Huelamo et al. (2011). For example, Kraus & Ireland (2012) detected a planet candidate in the LkCa15 transition disk. Multi-epoch observations of this object revealed orbital motion of the companion at the level of ∼4° a year (Ireland & Kraus 2014), strengthening the case for the protoplanet explanation of the phase signal.

We compare the predicted position of an orbiting companion to our fit results. Our best fit to the archival discovery data has a separation of 88 mas, which at T Cha's distance of 108 ± 9 pc (Torres et al. 2008) corresponds to 9.5 AU. We can use the orientation of the outer disk as well as the best-fit separation to predict the position of the planet in followup data sets.

We first take a circular orbit at the same inclination (i = 58°) and position angle (−70°) as the outer disk as determined by Olofsson et al. (2013). We choose the 2010 companion location as the initial orbital position, and then predict the projected separations for the followup data sets. Each panel in Figure 7 shows the orbital prediction for a single observation over a χ² slice at its best fit contrast ratio. While the predicted position is within 1σ of the NaCo 2011 L' best fit, all other data sets, even those with lower false alarm probabilities such as 2013 NaCo L' and Ks, rule out the presence of a companion on this orbit. Thus, this scenario cannot explain the observations.

Zoom In Zoom Out Reset image size

If we assume the NaCo 2012 L' signal is a real detection, then we need to match its fitted position along with the fits from other epochs. This would require at least one full orbit to be completed between 2010 and 2013. However, the observed projected separations (s) at each epoch place a lower bound on the apocenter distance:

$$\begin{equation} a(1+e)\ge s. \end{equation} \tag{ 8 }$$

Using this constraint and the fact that T Cha has a mass of 1.5 M_☉, we can check whether any orbits with a period of ∼3 yr could produce the observed separations. Even as e approaches 1, an orbit around a 1.5 M_☉ star cannot have a period of less than 8.5 yr. This rules out orbital motion as the cause of the different position angle in the 2012 L' data.

If we assume that the 2012 L' NaCo data missed the signal found in the other data sets, we can ask what orbits would cause the best fits to the remaining observations, which show motion compared to the Huelamo et al. (2011) model. For a grid of orbits, we calculated the predicted positions for the times of the 2010, 2011, and 2013 observations. We compared these positions to the 1σ confidence intervals in position angle and separation from the binary fitting. We take the results of a simultaneous fit to the 2013 NaCo and MagAO data as our constraint on the 2013 position. We find that inclined orbits, some of which cross into the outer disk, can produce the observations. Orbits with the same inclination as the outer disk require eccentricities higher than 0.9 to reproduce the observations. Figure 8 shows two example orbits over the 1σ confidence regions at the best fit contrast for the three epochs. While a low-eccentricity orbit in the plane of the disk is inconsistent with the observations, an orbit that is highly eccentric or substantially misaligned with the disk is compatible with the data. Since one would expect a young planet that is still accreting from the disk to be on a low-eccentricity, aligned orbit, these models do not seem physically likely.

Zoom In Zoom Out Reset image size

9.2. Forward Scattering from the Disk

9.3. Optically Thin Disk Asymmetries

9.4. Chance Alignment

We investigate the probability that the chance alignment of a foreground or background object would cause a companion signal in the 2010 NaCo data set, and that its proper motion would be mistaken for orbital motion between 2010 and 2013. Since T Cha has a proper motion of −39.61 mas yr⁻¹ in right ascension and −9.87 mas yr⁻¹ in declination, stationary objects in the foreground or background could appear to move like orbiting companions.

The number of chance alignments in a field of view with area A_FOV is

$$\begin{equation} n_{{\rm align}} = A_{{\rm FOV}} \Sigma, \end{equation} \tag{ 9 }$$

where Σ is the surface density of stars along the line of sight. We use the extent of T Cha's disk gap to define our field of view. From Olofsson et al. (2013), the gap extends from 0.17 AU to 12 AU, subtending ∼011 at 108 pc. To estimate the surface density of stars along the line of sight, we queried Two Micron All Sky Survey for all stars within 1° of T Cha, with ΔKs lower than 7 (Ks of ∼7 − 14 mag). This gives a stellar surface density of 11,544 stars deg⁻². The probability that a chance alignment with a foreground or background object would cause a companion signal in the NaCo 2010 data is then ∼7.0 × 10⁻⁶.

9.5. Systematic Errors

We presented multi-epoch observations of the T Cha transition disk, taken using VLT/NaCo and Magellan/MagAO/Clio2, in L', Ks, and H bands. Out of the nine data sets, three are too noisy to detect signals at the level of the companion candidate published in Huelamo et al. (2011); these are the 2010 and 2011 NaCo Ks and 2013 NaCo H band data.

We find companion parameters comparable to those published in Huelamo et al. (2011) in our re-reduction of the initial discovery data. Furthermore, we find comparable binary parameters to the companion candidate's in the 2011 NaCo L' data set, and tentative evidence for radial motion by the time of the 2013 NaCo and MagAO/Clio2 observations. Assuming that noise in the 2012 NaCo and 2013 MagAO/Clio2 data led to non-detections of the signal (simulations show such false non-detections to be fairly probably in these data sets), highly eccentric or misaligned orbits could result in a signal consistent with the observations.

Scattered light from T Cha's outer disk could provide another explanation for the observed kernel phases, although a binary model gives a slightly better fit to the data. Preliminary image reconstructions also suggest an asymmetric structure, perhaps consistent with disk–planet interactions, as the source of the observed signals. Lastly, the ratio of the NaCo 2013 Ks and L' best fit separations argues for the possibility that challenges associated with AO correcting a dim, southern source, could cause a systematic error that would masquerade as a close in companion. The detection of a secondary minimum in the MagAO/Clio2 data at the same position as the NaCo detections is encouraging, but follow-up observations with higher signal to noise are required to rule out the possibility of systematic errors.

This work was supported by NSF AAG grant $\#$1211329. J.R.M and K.M.M. were supported under contract with the California Institute of Technology (Caltech) funded by NASA through the Sagan Fellowship Program. This material is based upon work supported by the National Science Foundation Graduate Research Fellowship under Grant No. DGE-1143953. Any opinion, findings, and conclusions or recommendations expressed in this material are those of the authors(s) and do not necessarily reflect the views of the National Science Foundation.

Figure 9 shows the sky rotation coverage for the 2010 and 2011 NaCo Ks and 2013 NaCo H band data sets. The comparison of the snapshot errors (see Section 4) is displayed in Figure 10. The calibrated kernel phase histograms, with the Gaussian distributions used to generate the noise simulations discussed in Section 7, are shown in Figure 11.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

A.1. 2010 VLT/NaCo Ks Data: Published Non-detection

Table 5 lists the results of binary fits to these data as well as the 2011 NaCo Ks and 2013 NaCo H band data sets. The Δt⁻¹ calibration strategy produced a best fit position angle of $122\pm ^{66}_{137}$°, separation of 268$\pm ^{412}_{48}$ mas, and a contrast of $\Delta Ks = 4.8\pm ^{1.3}_{0.5}$ mag. Appendix B shows these data alongside this work's best fit and the Huelamo et al. (2011) binary model.

Table 5. Binary Fit Results

Data Set	P.A.	Separation	ΔL'	P(FA)a	P(FA)b	P(miss)	$\Delta _{\chi ^2,Hu}$	σHu, allowed
	(deg)	(mas)		(%)	(%)	(%)
2010 Jul 1: VLT/NaCo Ks
Δt−1	122 $\mathbf {\pm ^{66}_{137}}$	268 $\mathbf {\pm ^{412}_{48}}$	4.8 $\mathbf {\pm ^{1.3}_{0.5}}$	71.0	71.5	93.6	10.40	3
Δt−10	124 $\pm ^{165}_{155}$	259 $\pm ^{441}_{89}$	4.8 $\pm ^{1.8}_{0.6}$	36.0	92.6	95.5	8.19	3
2011 Mar 15: VLT/NaCo Ks
Δt−1	22 $\mathbf {\pm ^{0.03}_{165}}$	358 $\mathbf {\pm ^{342}_{8}}$	4.2 $\mathbf {\pm ^{0.6}_{0.4}}$	68.0	9.5	83.5	9.44	3
Δt−10	22 $\pm ^{1}_{2}$	359 $\pm ^{1}_{9}$	4.1 $\pm ^{0.6}_{0.3}$	31.7	5.7	82.9	9.31	3
2013 Mar 27: VLT/NaCo H
Δt−1	−27 $\mathbf {\pm ^{5}_{7}}$	31 $\mathbf {\pm ^{9}_{11}}$	3.8 $\mathbf {\pm ^{0.8}_{0.8}}$	1	< 0.1	84.1	32.01	> 4
Δt−10	−28 $\pm ^{5}_{7}$	29 $\pm ^{18}_{9}$	3.6 $\pm ^{0.9}_{0.6}$	0.8	<0.1	98.1	30.10	>4

Notes. ^aUsing distribution of noise simulation best fits. ^bUsing distribution of noise simulation F statistics.

Download table as:  ASCIITypeset image

For these data, there is a very high probability that the best fit resulted from noise, and that a companion signal would not have been recovered. The best fit to these data lies on a contour that encloses 29% of the noise simulation results (see Figure 12). There is thus a 71% chance that the best fit resulted from noise. The F statistics agree with this result; 71.5% of the F statistics are lower than that measured for the best fit, F = 0.843 (see Figure 13). Both of these suggest that the fit to this data set is consistent with noise.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

For the 2010 and 2011 NaCo Ks and 2013 NaCo H band data sets, we estimate our type II errors following the procedure outlined in Section 7. As the input signal, we add companions with the same separation (62 mas) and contrast (Δ = 5.1) as the Huelamo et al. (2011) binary. To account for biases due to sky rotation coverage, we vary the input position angle of the companions. Of the 1000 noise + signal realizations, 936 resulted in erroneous best fits. This gives a 93.6% chance that, with the noise present in the 2010 Ks data, we would not have recovered the companion signal had it been present. We thus cannot assign a non-detection to this data set with a high degree of confidence.

A.2. VLT/NaCo 2011 Ks Data

A.3. VLT/NaCo 2013 H Data

Figures 14–22 show plots of kernel phase versus scan index, a proxy for sky rotation angle, for all observational epochs. In each plot, the gray points show the data and the red solid line marks this work's best fit. Figures 14–19 show datasets with low enough noise levels to include in our Discussion. Figure 19, which displays the MagAO/Clio2 data, shows the NaCo 2013 L' best fit in green, since the two datasets were taken close in time to one another. All other figures in this section show the Huelamo et al. (2011) model in green. Figures 20–22 show datasets with high enough scatter to wash out signals of interest, and thus were not include in our discussion.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Our kernel phase projection is similar to that in Martinache (2010). We begin with the matrix A, which describes the ways in which N_a apertures (ϕ) are combined to yield N_p phases (Φ):

$$\begin{equation} \Phi = {\bf A} \cdot \phi. \end{equation} \tag{ C1 }$$

This equation can be modified for observations of a source with intrinsic signal, assuming that the source phase simply adds to the instrumental phase:

$$\begin{equation} \Phi = {\bf A} \cdot \phi + \Phi _0. \end{equation} \tag{ C2 }$$

In order to eliminate the instrumental phase, we are searching for a matrix, K, such that:

$$\begin{equation} {\bf K} \cdot {\bf A} = {\bf 0}. \end{equation} \tag{ C3 }$$

We can use singular value decomposition to find K. We decompose A^T in the following way:

$$\begin{equation} {\bf A}^T = {\bf U} \cdot {\bf W} \cdot {\bf V}^T, \end{equation} \tag{ C4 }$$

where U is an N_a × N_p column-orthogonal matrix, W is an N_p × N_p diagonal matrix with either positive or zero elements, and V is an N_p × N_p orthogonal matrix. The columns of V corresponding to zero W-values are filled into the rows of K. We then build K from a matrix, T, which describes how to combine phases into closure phases (Φ_cp):

$$\begin{equation} \Phi _{cp} = {\bf T} \cdot {\bf A} \cdot \phi + {\bf T} \cdot \Phi _0 = {\bf T} \cdot \Phi _0 \end{equation} \tag{ C5 }$$

We find B, such that:

$$\begin{equation} {\bf B} \cdot {\bf T} = {\bf K} \end{equation} \tag{ C6 }$$

and

$$\begin{equation} \Phi _{k} = {\bf B} \cdot {\bf T} \cdot {\bf A} \cdot \phi + {\bf B} \cdot {\bf T} \cdot \Phi _0 = {\bf K} \cdot \Phi _0. \end{equation} \tag{ C7 }$$

In the following equations, for any matrix M, M$_{{\rm right}}^{-1}$ and M$_{{\rm left}}^{-1}$ represent the right and left generalized inverses, respectively. These are used to invert non-square matrices; a right inverse is required when a matrix has full row rank but does not have full column rank, and a left inverse when a matrix has full column rank but does not have full row rank:

$$\begin{equation} {\bf M}_{{\rm right}}^{-1} = {\bf M}^T \cdot ({\bf M} \cdot {\bf M}^T)^{-1} \end{equation} \tag{ C8 }$$

$$\begin{equation} {\bf M}_{{\rm left}}^{-1} = ({\bf M}^T \cdot {\bf M})^{-1} \cdot {\bf M}^T. \end{equation} \tag{ C9 }$$

K has the dimensions N_k × N_p where N_k is the number of (linearly independent) kernel phases and N_p the number of Fourier phases (N_k < N_p). K has a right inverse:

$$\begin{equation} {\bf B} \cdot {\bf T} \cdot {\bf K}_{{\rm right}}^{-1} = {\bf K} \cdot {\bf K}_{{\rm right}}^{-1} = {\bf I}. \end{equation} \tag{ C10 }$$

Since B has the dimensions N_k × N_cp, (N_k < N_cp), B has a right inverse, which, according to C10, is the following:

$$\begin{equation} {\bf B}_{{\rm right}}^{-1} = {\bf T} \cdot {\bf K}_{{\rm right}}^{-1} \end{equation} \tag{ C11 }$$

and

$$\begin{equation} {\bf B} \cdot {\bf B}_{{\rm right}}^{-1} = {\bf I}. \end{equation} \tag{ C12 }$$

Furthermore, the left inverse of B$_{{\rm right}}^{-1}$ can be calculated to find the B in Equation (C12).

So:

$$\begin{equation} {\bf B} = \left({\bf B}_{{\rm right}}^{-1}\right)_{{\rm left}}^{-1} = \left({\bf T} \cdot {\bf K}_{{\rm right}}^{-1}\right)_{{\rm left}}^{-1}. \end{equation} \tag{ C13 }$$

This satisfies both Equations (C6) and (C7). We find the kernel phase covariance matrix, C_k, using the closure phase covariance matrix C_cp in the following way:

$$\begin{equation} {\bf C}_k = {\bf B} \cdot {\bf C}_{cp} \cdot {\bf B}^T. \end{equation} \tag{ C14 }$$

We can take our kernel phase variances to be the diagonal entries of C_k.