High Angular Resolution Imaging of the V892 Tau Binary System: A New Circumprimary Disk Detection and Updated Orbital Constraints

We present a direct imaging study of V892 Tau, a young Herbig Ae/Be star with a close-in stellar companion and circumbinary disk. Our observations consist of images acquired via Keck 2/NIRC2 with non-redundant masking and the pyramid wavefront sensor at K$^\prime$ band (2.12$\mu$m) and L$^\prime$ band (3.78$\mu$m). Sensitivity to low-mass accreting companions and cool disk material is high at L$^\prime$ band, while complimentary observations at K$^\prime$ band probe hotter material with higher angular resolution. These multi-wavelength, multi-epoch data allow us to differentiate the secondary stellar emission from disk emission and deeply probe the structure of the circumbinary disk at small angular separations. We constrain architectural properties of the system by fitting geometric disk and companion models to the K$^\prime$ and L$^\prime$ band data. From these models, we constrain the astrometric and photometric properties of the stellar binary and update the orbit, placing the tightest estimates to date on the V892 Tau orbital parameters. We also constrain the geometric structure of the circumbinary disk, and resolve a circumprimary disk for the first time.


INTRODUCTION
Circumbinary planets have only recently been observed (e.g.Sigurdsson et al. 2003;Doyle et al. 2011), despite long-standing predictions of planets existing in stable orbits around stellar binaries (e.g.Szebehely 1980; Rabl & Dvorak 1988;Benest 1993).Out of the ∼5000 confirmed exoplanets, only ∼75 have been confirmed to be in a circumbinary (CB) orbit, an orbital path where a planet orbits both the primary and secondary stars (e.g.Akeson et al. 2013).These mature CB planets are an enigma for planet formation theory, since their semi-major axes are close-in, near the limits of instability (Doyle et al. 2011;Socia et al. 2020).These closein separations suggest that migration is one of the key mechanisms in CB planet formation, but many aspects of this process remain untested (e.g Masset et al. 2006;Pierens & Nelson 2008;Kley & Haghighipour 2014;Penzlin et al. 2021;Coleman et al. 2022).Understanding the details of the locations and the timescales on which CB planets form requires better constraints on the youngest of these systems where planets are actively forming.
High angular resolution studies of circumbinary disks enable us to characterize the initial conditions of CB planet formation.Constraining their dust distributions can advance our understanding of where planet formation may be ongoing (e.g. in reservoirs of circumbinary, circumprimary, or circumsecondary material Müller & Kley 2012;Lines et al. 2015).Such studies can also inform our understanding of migration mechanisms and how they are influenced by the spatial properties of the CB disk (Guilloteau et al. 2008;Monnier et al. 2009;Boehler et al. 2017;Kurtovic et al. 2018).Furthermore, theoretical predictions of protoplanet spectral energy distributions (SEDs) suggest that they have low contrasts at infrared (IR) wavelengths (e.g.Zhu 2015; Eisner 2015; Baraffe et al. 2003).Hence, in addition to disk characterization, searching for and characterizing actively-forming or recently-formed planets in the IR would directly constrain the initial conditions of CB planet formation.

V892 Tau
One natural laboratory for studying circumbinary planet formation is V892 Tau, a young Herbig Ae/Be star located at a distance of ∼135 pc within the Taurus-Auriga star-forming region (Gaia Collaboration 2020).Spectral type estimates of the primary vary from A0 to B9 (Hillenbrand et al. 1992;Strom & Strom 1994;Hernández et al. 2004;Herczeg & Hillenbrand 2014).

Vides
When imaged at near-IR wavelengths, V892 Tau is shown to host an almost equally bright secondary companion at angular separations ranging between 40-65 mas, with the most recent constraints on the orbital parameters being: semi-major axis a = 7.1 ± 0.1 au, period P = 7.7 ± 0.2 yr, eccentricity e = 0.27 ± 0.1 and inclination i = 59.3 ± 2.7 • (Smith et al. 2005;Monnier et al. 2008;Long et al. 2021).From CO Keplerian gas rotation, the total mass of the system is determined to be 6.0 ± 0.2 M ⊙ (Long et al. 2021).
The circumstellar environment of V892 Tau and its companion has also been well studied in both the midinfrared and the millimeter continuum.The circumbinary disk was first discovered when imaged at 10.7 µm.An elongated structure with two bright lobes was detected and fit with a 2-dimensional Gaussian model to estimate the inclination of the disk (Monnier et al. 2008).In the millimeter continuum, V892 Tau has a radially asymmetric dust ring with peak mm emission at 0.2 ′′ and enough mass to form giant planets (Long et al. 2021).Warping in the CB disk has been tentatively identified in mm and 10.7 µm imaging (Monnier et al. 2008;Long et al. 2021).This scenario is consistent with the V892 Tau binary eccentricity and mass ratio, since highly eccentric orbits are known to induce warping and tidal truncation within the circumbinary disks of near equal mass binaries (Artymowicz & Lubow 1994;Hirsh et al. 2020;Miranda & Lai 2015).In addition to characterizing the CB disk, mid-IR long baseline interferometry has tentatively suggested the presence of a resolved circumprimary disk, but was unable to distinguish a circumprimary disk from an additional dusty companion (Cahuasquí 2019).

Paper Outline
We present 2-4 micron high angular resolution observations of the V892 Tau system.The paper is formatted as follows: Section 2 describes the observations and data reduction.Section 3 describes the image reconstruction method and analytical model fitting to the data.In Section 4, we statistically compare the results of each model and the correlations between the models and the data; then we place constraints on the geometry of the system.We then append our results to prior astrometric data and fit an orbit to the V892 Tau stellar companion.In Section 5, we discuss the implications of the results and estimate our sensitivity to additional, planetary-mass companions in the system.We conclude in Section 6.

Non-Redundant Masking
Nearby star-forming regions are located at distances ≳ 100 pc, where the angular resolution provided by traditional direct imaging methods is insufficient for IR planet searches on orbital scales of ≲ 10 − 15 AU (e.g.Guyon et al. 2013).Resolving smaller scales at such distances requires interferometric techniques.One method is non-redundant masking (NRM), an aperture masking technique in which a conventional telescope is turned into an interferometric array by placing a mask with discrete holes in the pupil plane (e.g.Tuthill et al. 2000).Each pair of holes, also known as a baseline, has a distinct orientation and separation, such that each baseline has its own unique spatial frequency (hence the term non-redundant).
The image on the detector, or the interferogram, shows the interference fringes formed by the mask.We take the Fourier transform of the interferogram to get the complex visibilites (which have the form Ae iϕ ).The complex visibility for each baseline is located in an extended region in Fourier space due to the finite size of the holes and the wavelength coverage set by the observing bandpass.From the appropriate regions in Fourier space, we calculate two quantities: squared visibilities and closure phases.Squared visibilities are the squares of the complex visibility amplitudes, and give the power corresponding to each baseline (e.g.Jennison 1958).Closure phases are sums of phases around baselines that form a triangle (e.g.Baldwin et al. 1986).Closure phases are highly sensitive to asymmetries and cancel first order wavefront errors, leaving only intrinsic phase and higher order residual errors.Closure phases and squared visibilities can be used to understand the source brightness distribution via both model fitting and image reconstruction.
Here we apply NRM on Keck 2/NIRC2 to deeply probe the structure of the V892 Tau circumbinary disk.

Observations
We used the 9-hole NRM in Keck 2/NIRC2 in conjunction with the Pyramid Wavefront Sensor (PyWFS) to directly image V892 Tau with the L ′ filter (central λ = 3.776 µm) and K ′ filter (central λ = 2.124 µm).Observations took place from 10:06 UT until 15:37 UT on November 6th, 2020 (L ′ ) and at 4:57 UT until 10:33 UT on January 21st, 2022 (K ′ ).The median seeing for the first half-night of November 5th, 2020 was 0.76 ′′ , with a minimum of 0.49 ′′ , a maximum of 1.26 ′′ , and a standard deviation of 0.16 ′′ as measured by the Differential Image Motion Monitor.The median seeing for the first halfnight of January 20th, 2022 was 1.11 ′′ , with a minimum of 0.75 ′′ , a maximum of 1.95 ′′ , and a standard deviation of 0.3 ′′ .
At both wavelengths, we observed V892 Tau for a halfnight centered on transit.We observed in vertical angle mode, which allows baselines to rotate on the sky.This fills the Fourier plane and allows for astrophysical signals to rotate while instrumental systematics remain fixed, enabling angular differential calibration.Figure 1 shows the rotation of each baseline with parallactic angle for the duration of each half-night at both wavelengths.We obtained ∼64 degrees and ∼108 degrees of sky rotation at L ′ and K ′ , respectively.
During each night we alternated between observing the science target and point spread function (PSF) calibrators, which are used to estimate higher order wavefront errors.To choose appropriate calibrators, we optimized between matching WFS brightnesses for similar quality adaptive optics (AO) correction, brightnesses at the science wavelength for efficient integration times, and separations on the sky to maximize common atmospheric paths and minimize slew times.We chose HD 283520, HD 281928 and HD 283577, whose coordinates and fluxes relative to V892 Tau are listed in Table 1.The calibrators have similar brightnesses to V892 Tau in the science bandpasses, allowing for efficient integration times.Despite the calibrators' brighter fluxes at H band (the PyWFS bandpass), the WFS frame rate for all objects was 1054 Hz, resulting in similar quality AO correction.The close angular separations between V892 Tau and the calibrators minimize slew overheads as well as calibration errors caused by differential refraction (Ireland 2013).
We subframed the 1024 × 1024 pixel detector to 512 pixels on each side and dithered on the detector, taking 10 frames in the top-left and bottom-right corners for each pointing to enable background subtraction.We spent equal amounts of time on the science target and PSF calibrators and alternated between dither-pair sequences.The coadds and integration times were chosen to build signal to noise with the detector in a linear response regime with minimal readout overheads.We obtained ∼40 min and ∼30 min of total integration time at L ′ and K ′ band, respectively.

Data Reduction
We use a well-tested pipeline, SAMpy (Sallum & Eisner 2017;Sallum et al. 2022), to reduce the data.We first calibrate the images by flattening and correcting for bad pixels by replacing them with the mean of the adjacent pixels.The median of one dither is then subtracted from each image in the other dither position for each dither-pair sequence for sky subtraction.The images are then cropped and Fourier transformed to obtain complex visibilities.
We crop the L ′ band images to 161 x 161 pixels and the K ′ band images to 91 x 91 pixels, then pad the images with zeros such that their sizes are 1024 x 1024 pixels before taking their Fourier transforms (FTs).We then sample the FT using all pixels that correspond to each baseline, and square the amplitude of the FT to obtain squared visibilities.The closure phases are then calculated such that the (u,v) coordinates of each closing triangle satisfy: (1) Sampling the FT at (u,v) coordinates satisfying the above equation, we calculate bispectra, which are products of the complex visibilities.We average the bispectra over multiple pixels for each baseline triangle for each frame, and then across the frames for each pointing.The phases of the bispectra are then taken to obtain the closure phases (Section 2.1).There are 36 squared visibilities and 84 closure phases calculated per pointing for the NIRC2 9-hole mask.
The squared visbilites and closure phases are next calibrated by fitting polynomial functions in time to the PSF calibrators.Polynomial orders are representative of instrumental noise variation, with a zeroth order polynomial indicating constant noise throughout the night and a high-order polynomial indicating high variability.We sample the polynomial function at the time of the science observations to estimate the instrumental systematics present in both the science closure phases and squared visibilities.To calibrate, we subtract the instrumental closure phases from the science closure phases, and divide the science squared visibilities by the instrumental squared visibilities.
We perform multiple calibrations with a variety of polynomial orders (ranging between zero and N-1, where N is the number of pointings).For the final calibrated closure phases, we adopt the order that minimizes their scatter, corresponding to a first-order polynomial at L ′ band and a second-order polynomial at K ′ band.For the squared visibilities, we assess the quality of the calibration by finding orders that minimize not only the scatter, but also the number of outliers with values > 1.A first-order polynomial satisfies these criteria at L ′ band,  and a fourth-order polynomial at K ′ band.The higher order polynomials at K ′ band suggest higher variability in the seeing, which is consistent with the seeing values reported in Section 2.2.We estimate error bars for the calibrated data by measuring the scatter of the calibrator squared visibilities and closure phases associated with each baseline and closing triangle, respectively.Rather than estimating the statistical error by measuring the standard deviation around the mean for each pointing, we estimate the calibration error by measuring the standard deviation across all pointings.This captures variability caused by changing systematics such as quasi-static speckles.We assign the estimated calibration error as the error bar for each science target squared visibility and closure phase.

Vides
This method is more appropriate than assigning statistical error bars, since calibration errors are the dominant error source in NRM observations (Ireland 2013).However, since some of the variations in the systematics are by definition removed during calibration, this method is conservative and tends to overestimate error bars.Lastly, we note that this approach means that error bars between baselines and closing triangles vary in a single pointing, but the error bars for all observables V892 Tau 5 associated with a given baseline or closing triangle are constant across all N pointings.

Image Reconstruction
After the observables are calibrated, we reconstruct images of the science target with SQUEEZE (e.g.Baron et al. 2010), an algorithm that uses Markov-Chain Monte Carlo (MCMC) methods to fit a model image to the closure phases and squared visibilities.SQUEEZE allows for simultaneous image reconstruction and model fitting, with several analytic model components that can be included in the reconstructed images.We reconstruct two images for each dataset using two different SQUEEZE models -a single point source model, and a binary model since V892 Tau has a known stellar companion.The single point source model includes a central, unresolved delta function to represent the central star, and we allow its fractional flux to vary.The binary model includes two unresolved delta functions to represent the central star and the companion, and we allow their fractional fluxes and the companion position to vary.We run SQUEEZE in parallel tempering mode to efficiently explore the image parameter space.The images have a platescale of 5 milliarcseconds (mas) per pixel and a size of 100 pixels on each side.

Geometric Models
Due to the sparsity of NRM Fourier coverage and incomplete recovery of phase information, image reconstruction is an under-constrained problem.We thus fit models to the observables to understand the morphology of the system and to test the robustness of the reconstructed images (Sallum & Eisner 2017, see also Section 4.3).We explore geometric models that include disk and unresolved companion components, fitting them to the Fourier observables.The three classes of models that we include are: (1) companion-only, (2) disk-only, and (3) disk-plus-companion.We fit them to each wavelength independently, since as geometric models they do not apply physically-motivated constraints on the relative fluxes of each component at K ′ and L ′ bands.
In the companion-only model, we analytically take the FT of two delta functions representing the primary and secondary stars, which have separation (S), position angle (PA), and the height of the secondary representing the contrast (CC).We convert the contrast from magnitudes to a flux ratio and give the delta function representing the primary a height of 1 and the delta function representing the secondary a height equal to the flux ratio.We let the separation, position angle, and contrast of the secondary vary between 0 and 500 mas, 0 and 360 degrees, and 0 and 8 magnitudes, respectively.We sample the model FTs at the locations corresponding to the mask baselines to calculate model closure phases and squared visibilities.To fully explore parameter space, we use emcee (Foreman-Mackey et al. 2013), an MCMC fitting package, in parallel tempering mode with 20 temperatures, 100 walkers, and 10,000 steps.
To model extended emission from the disk, we follow a procedure similar to that described in Appendix B of Sallum et al. (2021).The brightness distribution of the disk is defined as: where (3) and and x and y are locations in image space (increasing up and to the right).Here, θ is the position angle of the disk major axis, which is measured east of north and allowed to vary from 0 • to 180 • .The skew amplitude of the disk is given by A s and ranges from 0 to 1, and the peak skew position is given by ϕ s , which is allowed to vary between 0 • and 360 • measured east of north.The minor-to-major axis ratio (a ratio ) and the full-width-half-maximum (FWHM) of the Gaussian disk along the major axis are given by: We let the FWHM vary from 0 to 500 mas and allow for a hole that occupies a fraction f h of the FWHM.A delta function with a fractional flux b represents the central star.Both f h and b are free parameters between 0 and 1.We use emcee to explore parameter space using the same parallel tempering settings as the companion-only model.
The disk-plus-companion model is a combination of the Gaussian disk and the two delta functions.We vary the fractional fluxes occupied by the disk and companion.In the general disk-plus-companion model, the disk is allowed to be interior or exterior to the companion.

Vides
However, we also explore a model that forces the disk to be exterior to the companion which is discussed in more detail in Section 4.2.We use the same emcee parallel tempering settings as the previous two models.

Contrast Curve Generation
We generate contrast curves from companion-only models to place mass constraints on additional, undetected companions.We generate a grid of evenly spaced companion models ranging 0 to 500 mas in separation, 0 • to 360 • in position angle, and 0 to 8 magnitudes in contrast.We then fit the companion models to the residuals between the data and the best-fit disk-pluscompanion model, calculating a χ 2 for each model.To obtain a contrast curve, we average over the grid in position angle and then calculate χ 2 intervals between the null model (no companion) and the companion-only models at each separation and contrast.At a given separation, we adopt the contrast with a χ 2 interval of 25 as the 5σ contrast (e.g.Sallum et al. 2019).This method gives similar results as fitting the grid of companion models to a PSF calibrator that has undergone the calibration process (see Section 2.3).

Fourier Observables and Reconstructed Images
Figure 2 shows the Fourier observables -closure phases (left) and squared visibilities (right) -at both wavelengths.At L ′ band, the squared visibilities fall off rapidly as a function of baseline length independent of position angle, which is characteristic of a centrosymmetric, extended morphology.At K ′ band, some baselines fall off more sharply than others, indicating a lower degree of centro-symmetry.The closure phases at both bands have high values compared to the calibrators, which suggests asymmetry at close-in angular separations.
Figure 3 and Figure 4 show the reconstructed images that include the SQUEEZE single point source model and binary model, respectively.These images also suggest that the K ′ band morphology is less centrosymmetric than that at L ′ band.Approximately ∼ 36% of the non-stellar flux is enclosed in the compact signal located southwest of the star at K ′ band, compared to ∼ 6.4% in the L ′ band signal located southeast of the star.In Figure 4, the compact signals are removed from the images since they are captured by the SQUEEZE binary model.The K ′ image has a higher fraction of the total flux accounted for by the binary components, at 0.87 compared to 0.79 at L ′ .The residual extended emission in the L ′ band image and the removal of the southwest compact signal in the K ′ images are consistent with the single point source model images.
The SQUEEZE images show different morphologies between the datasets and the models used during reconstruction (single point source versus binary).The binary model places the companion at 26.9 ± 0.7 mas with a PA of 144.9 ± 2 • at L ′ band and 39.6 ± 0.03 mas with a PA 239.8 ± 1 • at K ′ band.The changes in the separation and position angle of the companion are due to orbital motion between the two epochs.We find agreement between the positions of the companion delta functions in the binary model images and the locations of the compact signals in the single point source model images, but the fractional fluxes of the central star differ between the two models.With the single point source model the fractional fluxes of the central stars are 0.71 and 0.50 at L ′ and K ′ band, respectively.Using the binary model, the fractional fluxes of the central stars are 0.51 and 0.62, with secondary fractional fluxes of 0.28 and 0.25, at L ′ and K ′ band, respectively.
As we further demonstrate in Section 4.3, these discrepancies are due to biases introducted in the fitting process.SQUEEZE has difficulty simultaneously matching the squared visibilities and closure phases when the single point source model is used (χ 2 =559.99).The binary model provides a better match to the data (χ 2 =327.85) and is also a more appropriate model given the known stellar companion.In the single point source model, achieving a better match to both observables would require arbitrary changes to the error bar scalings of the squaredvisibilities and closure phases (e.g.Sallum & Eisner 2017).Since the error bar scaling process would be motivated not by the data but by the SQUEEZE algorithm, we refrain from doing this and instead show all reconstructions using original error bars.
Although the use of the binary model is more physically motivated, we include both SQUEEZE models to demonstrate the effects of adding different components during the reconstruction process.While the single point source model allows us to more freely place circumstellar emission at any location in the images, the binary model images reveal complex structure in the disk that is not visually apparent in the single point source images (Figure 4).The K ′ band image shows a point-like feature to the southeast of the central star, and the L ′ band image shows complex asymmetries in the form of multiple arcs.To explore these features in the context of different disk plus companion scenarios, we perform simulated image reconstructions described in Section 4.3.ities and closure phases at both wavelengths are best described by the disk-plus-companion model.To as-sess the significance of the preference for the disk-pluscompanion model over the others, we compare the im-V892 Tau 9 provements in χ 2 values between models to distributions with N degrees of freedom, where N is the difference in number of parameters between two models.We find that the disk-plus-companion model is preferred at >5σ for both wavelengths.

Best-fit Geometric Models
At K ′ band, we find that the reduced χ 2 values are closest to 1 for the disk-plus-companion model, ranging between 1 and 2 for the two Fourier observables.At L ′ band, the individual disk-plus-companion reduced χ 2 values for the closure phases and squared visibilities are < 1, which would imply over-fitting for perfect error bars.However, the conservative error bar calculations applied here (see Section 2.3) may bias the reduced χ 2 toward low values.We thus base the χ 2 model selection primarily on the raw, rather than reduced values.Figure 5 shows the correlation between the disk-pluscompanion model and the data for both observables and bands.Given the conservative error bars, the χ 2 model selection tests, and these correlations, we accept the disk-plus-companion model as the best fit for both wavelengths.
Figure 6 shows the best fit models at L ′ band and K ′ band.At L ′ band, the geometry of the system is best described by a circumbinary disk and companion.We find that the FWHM of the CB disk is 189.9 +16.5/-19.7 mas (∼ 25.5 AU) and the separation of the stellar companion at the time of observation is 26.0 +0.7/-0.6 mas located at 147.4 • ± 1.4 • measured east of north.The contrast of the stellar companion is 0.60 ±0.03 magnitudes relative to the host star.From the total flux of the system and the contrast, the flux of the secondary companion is calculated and converted to an absolute magnitude, M L = 6.63 ±0.03 mag.
At K ′ band, we find that the geometry of system is best represented with a circumprimary disk and companion.We detect a circumprimary disk with a FWHM of 15.7 +2.3/-2.0mas (∼ 2 AU), and a companion with a separation of 42.1 +0.70/-0.63mas (∼ 5.6 AU) located at 238.51 • +0.98/-0.83• east of north.The contrast of the secondary star relative to the primary star is 0.67 +0.02/-0.06mag, giving M K = 6.49+0.02/-0.06mag.The skew of the circumprimary disk is roughly aligned with the PA of the companion.
To test that the detection of the circumprimary disk is a physical feature of V892 Tau and rule out a local minimum in the fitting, we explore a set of models where we set the upper bound of the companion separation (S) prior to: where α = PA − θ. (8) In the above equations, f h is the fraction of the semimajor axis that is occupied by a hole and θ is the position angle of the major axis of the disk.These equations ensure that companion is always within the disk, forcing a CB disk.Table 3 lists the results of the forced circumbinary disk-plus-companion model fit to the K ′ band data.The χ 2 values indicate that this model does not adequately describe the data.The discrepancy between the data and model is especially apparent in the closure phases, where the reduced χ 2 value, χ 2 r , is 5.11 for the forced CB model, as opposed to 2.18 for the unrestricted diskplus-companion model.The χ 2 difference between the unrestricted disk-plus-companion model and the forced CB disk model is 2240.16.We compare this value to significance estimates for a χ 2 distribution with 10 degrees of freedom, the number of parameters in the CB disk-plus-companion model.The circumprimary disk scenario is preferred at >5σ.Furthermore, the forced circumbinary parameters are poorly constrained due to the existence of many local likelihood maxima, which also include non-physical scenarios.We thus find the circumprimary disk detection to be robust.

SQUEEZE Model Reconstruction
Image reconstructions capture the true source brightness distribution to varying degrees depending on the mask (u,v) coverage, the amount of sky rotation in the observations, and the particular regularization choices in the individual algorithms (e.g.Sallum & Eisner 2017).We thus check whether the best-fit geometric model reproduces the structure in the images reconstructed from the observations (Figure 3 and Figure 4).We generate model closure phases and squared visibilities by sampling the best fit geometric models with the same sky rotation and (u,v) coverage as the observations.We then add noise to the model closure phases and squared visibilities such that the scatter matches that of the observations.Figures 7 and 8 show the resulting images, which are generally consistent with the reconstructions generated from the data (Figure 3 and 4).
The fractional fluxes of the central stars in the simulated reconstructions using SQUEEZE's single point source model (0.79 at L ′ band and 0.51 at K ′ band) are similar to those for the data (0.71 at L ′ band and 0.50 at K ′ band).In the reconstruction using SQUEEZE's binary model, the fractional fluxes are also similar to the observations, with central star fractional fluxes of 0.57 at L ′ band and 0.65 at K ′ band.The fractional fluxes of the companions are 0.31 at L ′ band and 0.27 at K ′ band.This demonstrates in a controlled way that the variation in fractional fluxes for different SQUEEZE configurations is an algorithmic bias.Furthermore, the consistency between the fractional fluxes for the simulations and the observations shows that approximately the same algorithmic biases are introduced in the reconstructions generated from the data and from the geometric models.

Vides
We find that when we reconstruct the image from the best-fit disk plus companion observables using a single point source SQUEEZE model, the model reconstruction visually matches the data at both bandpasses.However, this is not the case at both bands when we reconstruct images for the disk plus companion model observables using SQUEEZE's binary model.At K ′ band, we find that the best-fit circumprimary disk model reproduces the data, including the point-like blob to the The best fit geometric model at L ′ band.The FWHM of the circumbinary disk is 189.9 +16.5/-19.7 mas and the separation of the stellar companion is 26.84 +0.75/-0.64mas at 146.6 +1.35/-1.32degrees measured east-of-north.The image has been sub-framed and smoothed with a Gaussian to make the δ functions visible.Right: The best fit geometric model at K ′ band.The circumprimary disk has FWHM of 15.7 +2.3/-2.0mas, and a separation of 42.1 +0.70/-0.63mas located at 238.51 +0.98/-0.83degrees measured east-of-north.The image has been sub-framed and smoothed with a Gaussian to make the δ functions visible.Changes in the separations and position angles of the companions between the two epochs show the orbital motion of the stellar companion.The locations of the companion in both panels are similar to those in Figure 4; however their appearances are stretched here due to differences in the images' fields of view.

L Band
K Band We find that the orbit is consistent with the astrometry estimated by the SQUEEZE images as well.
We update the orbit of the V892 Tau stellar companion using Orbitize!(e.g.Blunt et al. 2020).We fit our astrometry and the archival data shown in Table 4, which lists the separations and position angles, measured east of north, of each data point included in the orbit fit.We exclude the first data point in Smith et al. (2005), following a similar method as Long et al. (2021).This data point was ambiguously measured; it was unclear which component was the primary or secondary star, making the uncertainty in the position angle 180 • .
We use parallel-tempered MCMC methods to fit an orbit to the data, with 10 temperatures, 100 walkers, and 10,000 steps.We use the Orbitize!default uniform priors; varying the semi-major axis from 0.001 to 10 7 AU, the eccentricity from 0 to 1, the inclination from 0 to 2π, the argument of periastron from 0 to 2π, the position angle of the ascending node from 0 to 2π (measured E of N), and the periastron passage from 0 to 1.The distance and total mass of the system are Gaussian priors that are centered on the parallax (7.4 ± 0.08 mas; L Band 11-05-2020 022) L Band 11-05-2020 Figure 10.Contrast curves for V892 Tau at L ′ band (left) and K ′ band (right) using NRM and the PyWFS.We find that the PyWFS contrast is 0.5-1 magnitudes better than observations with NIRC2 NRM behind the Shack-Hartmann wavefront sensor.We use these contrast curves to place limits on undetected companions and find that we are sensitive to ≳20 MJ planets at L ′ band and ≳50 MJ planets at K ′ band (assuming an age of ∼ 1 Myr).Gaia Collaboration 2020) and total mass of the system (6.0 ± 0.2 M ⊙ ; Long et al. 2021), respectively.Parallax uncertainties for binaries with separations less than a few arcseconds and G ≲ 18 have been underestimated by ≲ 30% (El-Badry et al. 2021).Thus, the parallax error is likely slightly underestimated due to the presence of the binary and disk.However, since the astrometric error bars dominate the orbital error budget, even a 30% inflation of the parallax error would not significantly change the orbit fit results.
The left panel of Figure 9 shows 50 random orbits in the posterior distribution with the color bar representing a single orbit.The right panel shows the separation and position angle of V892 Tau's stellar companion across all the epochs.The semi-major axis of the stellar companion is 6.8 +0.04/-0.03AU (49.0 ± 2 mas) with a period of 7.2 +0.07/-0.05years.The eccentricity of the orbit is 0.26 ±0.04 and its inclination is 58.4 • ± 3 • .From these estimates, we independently measure the dynamical mass of the system to be 6.1 +0.2/-0.1 M ⊙ , which is in agreement with measurements made by Long et al. (2021).We compare our values to previous estimates in Table 5.

Additional Companions
Figure 10 shows the 1-5σ contrast curves at L ′ band (left) and K ′ band (right).We convert the contrast limits to absolute magnitudes and compare them to magnitudes predicted by hot start models with an assumed age of 1 Myr (approximately the same age as V892 Tau e.g.Baraffe et al. 2003).We estimate that we are sensitive to ≳20 M J brown dwarf companions at L ′ band and ≳50 M J companions at K ′ band, placing a rough upper limit on planetary masses in the system.We also convert the L ′ band contrast curve to planet mass times accretion rate (e.g.Eisner 2015).The L ′ band data are sensitive to a planet mass times accretion rate of ∼4 x 10 −5 M 2 J /yr, corresponding to a rapidly-accreting Jupiter analog or less-rapidly-accreting higher mass planet.

System Geometry
The geometry of the V892 Tau circumbinary disk measured here is consistent with the literature and we detect a new component of the disk structure with the discovery of the circumprimary.The diameter of the CB disk is in agreement with previous geometric constraints at ∼26 AU (e.g.Liu et al. 2007;Monnier et al. 2008).From the best-fit L ′ band model, we find that the position angle and orientation of the disk are similar to prior estimates made by Monnier et al. (2008) and Long et al. (2021).Previous inclination estimates (i disk = 54.6 • ) indicate that the northwest side of the disk is closest to the observer Long et al. (2021).With this orientation, the inside of the disk rim on the southeast side of the star would be most visible to the observer.This viewing angle effect, possibly combined with a puffed-up disk rim due to heating by the near-equal mass binary, may contribute to the asymmetry to the southeast of the star in the reconstructed images and geometric models.

Vides
From the L' reconstructed images, we find tentative evidence of warping in the circumbinary disk.In Section 4.3 (Figure 8), we show that the Fourier observables from the geometric model cannot reproduce the complex asymmetry in the L ′ image reconstructed from the data (Figure 4).Near-equal mass binaries with highly eccentric orbits have been shown to cause warped circumbinary disks (Artymowicz & Lubow 1994;Hirsh et al. 2020;Miranda & Lai 2015).This scenario is consistent with the V892 Tau eccentricity measurements in prior studies and this work (Table 5).
We further inform the architecture of the V892 Tau system with the first detection of a circumprimary disk with a diameter of ∼ 2 AU.Observations from Cahuasquí (2019) tentatively suggested a circumprimary disk, with differential phases in N band Mid-Infrared Interferometric Instrument (MIDI) data preferring either a circumprimary disk or an additional dusty companion.The K ′ band imaging presented here provides high enough angular resolution to firmly distinguish between these two physical interpretations.The rough alignment between the position angle of the companion and the skew angle of the circumprimary disk could suggest heating of the disk by the companion.Follow-up high angular resolution observations could identify whether the stellar companion causes the observed asymmetry by constraining the time evolution of the disk skew angle and the companion position angle.This work cannot place mass constraints on the circumprimary disk, but it nonetheless represents another potential reservoir of material for planet formation around V892 Tau.

Orbit
From the best-fit geometric models (discussed in Section 4.2), we measure the separation of the stellar companion to be 3.5 ± 0.1 AU at L ′ band.At K ′ band, we find the separation of the of the stellar companion is 5.6 AU ± 0.1 AU.The astrometric measurements are more precise at K ′ band than at L ′ band, since the angular resolution is higher due to the shorter wavelength.V892 Tau is also relatively bright at K ′ band (3.23 Jy at K ′ band and 1.75 Jy at L ′ band) and the K ′ band sky background is low.The high angular resolution, high signal-to-noise ratio, and AO correction with the Py-WFS significantly reduces the size of error bars on the observables at both wavelengths, which are then propagated statistically through the best-fit astrometry and photometry.The eccentricity and inclination estimates are within 1σ of previously-published constraints, while the semi-major axis and period have a discrepancy of 2σ compared to those studies (Long et al. 2021).We also find that our independent mass measurement of the system (6.1 +0.2/-0.1 M ⊙ ) is consistent with Long et al. (2021).

Benchmarking NRM with the PyWFS
Since V892 Tau is bright at H band (the wavefront sensing bandpass), we see an improvement in contrast when compared to observations with the Shack-Hartmann wavefront sensor (SH WFS) at both science wavelengths, allowing us to make high-precision astrometric measurements.We compare the PyWFS contrast limits in Figure 10 to contrast limits from Sallum & Skemer (2019) for a star with a similar R band (the SH WFS bandpass) magnitude to V892 Tau that was observed with Keck2/NIRC2 NRM behind the SH WFS.Like V892 Tau, this star is faint at R band, but bright at H band.We find that the contrast is ∼0.5-1 magnitudes better with PyWFS than with the SH WFS.The boost in contrast demonstrates that the PyWFS is beneficial for observing red young stars.These observations are the first benchmark of NRM with Keck's Py-WFS.Future observations of stars with fainter H band magnitudes will further quantify its performance in the lower-Strehl regime.

CONCLUSION
We presented multi-wavelength, multi-epoch Keck data of the V892 Tau circumbinary disk with NRM and the PyWFS.The data allow us to differentiate the secondary stellar emission from disk emission and deeply probe the structure of the disk at small angular separations.We fit geometric models to the L ′ and K ′ band data, and find that the morphologies of both images are best described by disk-plus-companion models.At L ′ V892 Tau 15 band, the circumbinary disk properties are consistent with results from prior studies.At K ′ band, we make the first robust detection of a circumprimary disk.From the properties of the stellar binary, we update the orbit using our own data and archival data.This work places the tightest constraints on the orbital parameters of the V892 Tau stellar companion and the geometric structure of the circumbinary disk.Future observations of the V892 Tau system may include additional monitoring of the circumprimary disk to determine whether its skew is caused by heating from the stellar companion.
These are the first published observations using NRM and the PyWFS in conjunction, providing a valuable contrast benchmark.We place mass constraints on undetected companions and compare the achieved contrast with the PyWFS to contrast limits achieved with the Shack-Hartmann WFS, finding a ∼0.5-1 magnitude boost in performance with the PyWFS.The exquisite AO correction (and thus achievable contrast) offered by the PyWFS enabled the precise astrometric measurements for V892 Tau, which improve on previous estimates by a factor of 10, and the high-angular-resolution detection of the circumprimary disk.These contrast benchmarks and the high-precision detections in the V892 Tau system demonstrate that future NRM surveys can take advantage of the PyWFS for observations of similarly young, red stars.

Figure 1 .
Figure 1.Scattered points show uv-coverage (the rotation of each baseline with parallactic angle) at L ′ band (left) and K ′ band (right).We obtained ∼64 • and ∼108 • of sky rotation at L ′ and K ′ , respectively.

Figure 2 .
Figure 2. V892 Tau Fourier observables.Left panels show closure phases for L ′ band (top) and K ′ band (bottom).The red points with error bars show the closure phases calculated from the data as described in Section 2.3.Hollow black circles show the closure phases calculated from the best-fit geometric model as described in Section 3.2.Right panels show the squared visibilities at L ′ band (top) and K ′ band (bottom).The red points with error bars show the squared visibilites calculated from the data as described in Section 2.3.Hollow black circles show the squared visibilities calculated from the the best fit geometric model as described in Section 3.2.bution:IM Lup

Figure 3 .
Figure 3. SQUEEZE images reconstructed from the L ′ band (left) and K ′ band (right) observables.The central stars are represented with δ functions occupying fractional fluxes of 0.71 and 0.50 at L ′ band and K ′ band, respectively.These fractional fluxes differ significantly from those estimated by the geometric models.This is not a physical feature of V892 Tau, but rather a bias in the SQUEEZE algorithm that we quantify using simulations (Sections 4.1 and 4.3).Table 3 lists the best-fit parameters and goodnessof-fit metrics for each of the three model scenarios

Figure 4 .
Figure 4. SQUEEZE images reconstructed with a binary model from the L ′ band (left) and K ′ band (right) observables.The central stars and companions are represented with δ functions equal to their fractional fluxes and the images are normalized appropriately.Changes in the separations and position angles of the δ functions show the orbital motion of the stellar companion between the two epochs.The fractional fluxes of the central stars are 0.51 and 0.62 at L ′ band and K ′ band, respectively.The fractional fluxes of the secondary components are 0.28 and 0.25 at L ′ band and K ′ band, respectively.

Figure 5 .
Figure 5. Top left: Scattered points with error bars in the left panels show the correlation between the closure phase data and the best-fit model closure phases at L ′ (top) and K ′ (bottom).Scattered points with error bars in the right panels show the correlation between the squared visibility data and the best-fit model squared visibilities at L ′ (top) and K ′ (bottom).In all panels the orange solid line represents a 1:1 correlation.

Figure 6 .
Figure 6.Left: The best fit geometric model at L ′ band.The FWHM of the circumbinary disk is 189.9 +16.5/-19.7 mas and the separation of the stellar companion is 26.84 +0.75/-0.64mas at 146.6 +1.35/-1.32degrees measured east-of-north.The image has been sub-framed and smoothed with a Gaussian to make the δ functions visible.Right: The best fit geometric model at K ′ band.The circumprimary disk has FWHM of 15.7 +2.3/-2.0mas, and a separation of 42.1 +0.70/-0.63mas located at 238.51 +0.98/-0.83degrees measured east-of-north.The image has been sub-framed and smoothed with a Gaussian to make the δ functions visible.Changes in the separations and position angles of the companions between the two epochs show the orbital motion of the stellar companion.The locations of the companion in both panels are similar to those in Figure 4; however their appearances are stretched here due to differences in the images' fields of view.

Figure 7 .Figure 8 .Figure 9 .
Figure 7. SQUEEZE images reconstructed from the disk-plus-companion model closure phases and squared visibilities at L ′ band (left) and K ′ band (right) using a single point source model.The fractional fluxes of the central stars are 0.79 at L ′ band and 0.51 at K ′ band.This approximately matches the fractional fluxes in the SQUEEZE reconstructions of the data.southeast of the star.We rigorously test this by reconstructing the image from the best-fit binary geometric model, and find that the point-like blob in Figure8is removed from the image.At L ′ band, the simulated reconstruction lacks the complex structure discussed in

Table 1 .
Properties of V892 Tau and PSF calibrators

Table 2 .
Observational setup for V892 Tau and PSF calibrators at both wavelengths Target

Table 3 .
The maximum likelihood models for each scenario.Degrees of freedom is indicated by DOF.Closure phases and squared visibilites are denoted by CP and V 2 , respectively.Chi-squared and reduced chi-squared values are indicated by χ 2 and χ 2 r , respectively.

Table 4 .
V892 Tau astrometry, showing archival constraints and the L ′ and K ′ band astrometry presented here.We omit the first row from orbital fitting due to a 180 • ambiguity in the position angle.