Two Rings and a Marginally Resolved, 5 AU, Disk Around LkCa 15 Identified Via Near Infrared Sparse Aperture Masking Interferometry

Sparse aperture masking interferometry (SAM) is a high resolution observing technique that allows for imaging at and beyond a telescope's diffraction limit. The technique is ideal for searching for stellar companions at small separations from their host star; however, previous analysis of SAM observations of young stars surrounded by dusty disks have had difficulties disentangling planet and extended disk emission. We analyse VLT/SPHERE-IRDIS SAM observations of the transition disk LkCa\,15, model the extended disk emission, probe for planets at small separations, and improve contrast limits for planets. We fit geometrical models directly to the interferometric observables and recover previously observed extended disk emission. We use dynamic nested sampling to estimate uncertainties on our model parameters and to calculate evidences to perform model comparison. We compare our extended disk emission models against point source models to robustly conclude that the system is dominated by extended emission within 50 au. We report detections of two previously observed asymmetric rings at $\sim$17 au and $\sim$45 au. The peak brightness location of each model ring is consistent with the previous observations. We also, for the first time with imaging, robustly recover an elliptical Gaussian inner disk, previously inferred via SED fitting. This inner disk has a FWHM of ~5 au and a similar inclination and orientation as the outer rings. Finally, we recover no clear evidence for candidate planets. By modelling the extended disk emission, we are able to place a lower limit on the near infrared companion contrast of at least 1000.

Protoplanetary disks are the expected sites of planet formation and are thus prime targets for young planet searches. Many protoplanetary disks, especially transition disks, have observable substructure in the millimeter (mm) dust continuum and molecular lines (Huang et al. 2018;Andrews et al. 2011;van der Marel et al. 2016) that is similar to that expected from planet-disk interaction theory (Pinilla et al. 2012). Observations of these structured systems at high angular resolution in the near infrared allow us additional characterization 1 Released on XX, XX, XXXX of the inner dust disk and potentially the discovery of planets in formation. Beyond just the discovery of forming planets, it is expected that these data will increase our understanding as to the accretion of material onto the planet in connection with the transport of material through the disk gap (e.g. Lubow & D'Angelo 2006).
LkCa 15 is a nearby (160 pc), young, T Tauri star with spectral type K2, mass 1.32 M , and luminosity 1.3 L (Francis & van der Marel 2020) surrounded by a protoplanetary disk with a central dust depleted cavity and three wide-orbit, narrow dust rings observed in the mm at ∼47, ∼69, and ∼100 au (Facchini et al. 2020). LkCa 15 was classified as a prototype for the pretransition disk class of objects by Espaillat et al. (2007a) due to the presence of a significant near-infrared excess over the stellar photosphere in the spectral energy distribution, suggesting a compact and optically thick inner disk. An unresolved inner disk may also be present at mm wavelengths, as deep ALMA observations reveal significant but low-level emission within the cavity located interior to the ∼47 au ring (Facchini et al. 2020).
Previous studies of LkCa 15 utilizing near infrared sparse aperture masking interferometry (SAM) claimed to have observed planets in the process of formation embedded in its disk. Kraus & Ireland (2011) reported a multi-epoch detection of a single blue compact source (deprojected orbital radius ∼16-20 au) surrounded by redder material using SAM on Keck II. Sallum et al. (2015) reported the detection of three point sources (LkCa 15 b,c,d with deprojected orbital radii between ∼15-19 au) using SAM on the Large Binocular Telescope. They also detected LkCa 15b in the accretiontracing Hα emission using SDI with the Magellan Adaptive Optics System. More recent direct imaging work, however, using SPHERE ZIMPOL and IRDIS (Thalmann et al. 2015(Thalmann et al. , 2016 and Subaru/SCExAO coupled with CHARIS, and Keck/NIRC2 data have disputed some of these results (Currie et al. 2019), by showing that the SAM detections of LkCa 15 b,c,d can likely be attributed to disk emission from a ring with an inner radius of ∼20 au.
Sparse aperture masking interferometry (SAM; Baldwin et al. 1986;Haniff et al. 1987;Tuthill et al. 2000;Greenbaum et al. 2015) is a high resolution optical interferometric observing technique that uses a mask (at the pupil plane) to transform the full telescope aperture into a set of smaller sub-apertures. Each sub-aperture pair acts as a two station interferometer, sampling a discrete component of the Fourier transform of a target's brightness distribution. SAM can resolve objects with tighter separations than can be resolved with a filled aperture (Lloyd et al. 2006;Sallum et al. 2019), making it a good technique for observing planets at close separations from their stars. Due to the sparse nature of SAM data (tens of data points per exposure), accurate image reconstruction is difficult and thus inner disk signals can be confused with planet detections (Cieza et al. 2013;Kraus et al. 2013).
The van Cittert-Zernike theorem (van Cittert 1934;Zernike 1938) states that the complex visibility to go from the image plane to the pupil plane is given by the two dimensional Fourier transform V (u, v) = I(x, y) exp(−i2π(ux + vy))dydx. (1) where u and v are the coordinates in the Fourier domain (measured in units of the observing wavelength), x and y are the coordinates in the image plane, and I(x, y) is the target brightness distribution. It is not currently possible to obtain well-calibrated phases on long baselines from ground based observations due to calibration challenges. Thus, more robust quantities are often derived from the individual complex visibility measurements, with an accompanying slight loss in information content. The first robust SAM observable is the squared visibility, or the modulus squared of the complex visibility. The second observable is the closure phase. The closure phase is a source-dependent quantity made by summing phases around triangles formed by holes in the mask. The closure phase is useful because the phase errors introduced by the telescope can be removed (Monnier 2007); phase shifts produced by noise sources, such as the atmosphere or electronics, will cause equal but opposite phase shifts for connected baselines in each closure triangle. As previously noted, a difficulty of working with SAM observations is that the data can be sparse, only consisting of a limited number of squared visibilities and closure phase measurements. This makes reconstructing an image of arbitrary complexity difficult, if not impossible. The lack of an easily obtainable, reliable image reconstruction becomes especially challenging when prior knowledge about a target is limited. Given sufficient prior knowledge about the types of objects being observed, however, a powerful approach for analyzing sparse interferometric data can be to use geometrical models. Geometrical models are constructed using analytic descriptions to model the major components in the scene using only a limited number of parameters.
Here we present an analysis of new near infrared SAM observations obtained by VLT/SPHERE-IRDIS (Cheetham et al. 2016) of the transition disk LkCa 15. We use Bayesian model fitting techniques to determine the structure of the inner tens of au around LkCa 15, testing the limits of using these high resolution techniques for observing the inner regions of protoplanetary disks. We also discuss the application of these methods to observing planets within transition disks. The structure of the paper is as follows: § 2 outlines how the data were obtained and processed; § 3 details the techniques used for image reconstruction, including the geometrical models that were fit; § 4 outlines how the data were analysed using these techniques and the results of applying these techniques; § 5 discusses the image reconstructions and performs a companion analysis, where point source models are considered in combination with the more complex geometrical models. We close with a summary in § 6.

OBSERVATIONS
LkCa 15 and calibrator star HD 284581 were observed at 2.11 (2.1) and 2.25 (2.3) µm (K12 dual-band imaging filters, Vigan et al. 2010) on January 4, 2018, using the 7 hole (21 baseline) VLT/SPHERE-IRDIS (Beuzit et al. 2019;Dohlen et al. 2008) SAM mode (Lacour et al. 2011;Tuthill et al. 2010;Cheetham et al. 2016), over the course of about two hours. The spatial scales covered by our observations (see Figure 1), defined by the longest and shortest baselines, are ∼70-260 mas or 10-40 au at the distance of LkCa 15. Note that the standard diffraction limit of a telescope with an aperture the size of the longest baseline would be 83 mas. Due to the interferometric nature of our data, we anticipate sensitivity to scales as small as ∼35 mas or ∼5 au. The LkCa 15 data set comprises four exposures, with each consisting of four integrations. Three exposures with the same number of integrations and frames per exposure were acquired for the calibrator HD 284581.
The IRDIS data were reduced using the SPHERE Data Reduction and Handling (DRH) automated pipeline (Pavlov et al. 2008). We applied basic corrections for bad pixels, dark current, and flat field, as well as distortion correction. Absolute orientation and pixel scales of the images were calculated using the parallactic angle and the absolute calibration provided by the SPHERE consortium. The current respective estimates of the pixel scale and true north for IRDIS are 12.255±0.009 mas/pixel and TN=-1.75±0.08 • (Maire et al. 2016). Our observations were obtained using the pupil-stabilized mode. The absolute orientation (northup, east-left) is therefore retrieved within the DRH using both TN and pupil offset (-135.99±0.11 • ).
AMICAL (Soulain et al. 2020) was used to extract the closure phase, squared visibility, and uncertainties from each integration, as well as to calibrate the data. The Fourier extracting method applied within AMICAL required additional reduction steps to allow an optimal extraction. The processed data (post-DRH) were background subtracted, cropped, and cleaned of residual bad pixels and cosmic rays. The cleaned cubes were then windowed with a super-Gaussian function of the form, exp(−ar 4 ), before closure phases and visibilities were extracted from the Fourier transforms of the images. Each complex visibility (over 21 baselines) is computed over a sub-sample of pixels using a weighted average. These pixel weights are retrieved from the absolute positions of the 7 apertures (in meters from the pupil plane) and the pixel size of IRDIS. The appropriate combination of Fourier transform for each aperture pair takes into account the fraction of pixel not centered on the real pixel position and increases the accuracy of the method 1 . A σ-clipping selection was applied to reject eventual bad frames due to bad seeing conditions or AO instability. Calibration of the closure phases (visibilities) was performed on each wavelength individually by subtracting (dividing) a weighted sum of the corresponding measurements taken on the nearest in time integration of the calibrator star. The extracted observables were finally saved as standard interferometric file format oifits (Duvert et al. 2017), taking into account the astrometric orientation obtained from the DRH (i.e., rotation of the u, v coordinates). The observables extracted from one of the 2.1 µm integrations were discarded due to a large variance from the mean compared to the other extracted data points.
In order to optimally utilize the available data, in our analysis we assume that the colour variations between 2.1 and 2.3 µm across the components of LkCa 15 are not significant and therefore that the benefit of increasing uv-plane coverage (see Figure 1) outweighs the additional reconstruction uncertainty obtained by considering the two wavelengths together during image reconstruction.

SAM ANALYSIS METHODS
Sparse aperture masking interferometry is a technique well suited for observing multiple point sources, such as a binary star system or a star orbited by a bright planet. Such sources have a distinct signature in the Fourier domain that is simple enough to often only require tens of samples to robustly observe. For two point sources (with one at the phase center) the complex visibilities at the observed coordinates in the u, v plane are given by the four parameter formula where α is the central star brightness, β is the off-center source brightness, and x 0 and y 0 are the coordinates of the off-center point source in the image domain.
Interpreting SAM data can be significantly more challenging when complex structures beyond point sources are being observed. Image reconstruction from sparse data is an ill-posed, under-constrained inversion problem, that can be addressed with the use of priors. One method for using strong prior knowledge about a target is to utilize image domain geometrical models. Instead of having to determine complex analytic Fourier transforms for any arbitrary geometry, the fast Fourier transform can be used. Here, the complex visibilities of geometrical models are calculated using a fast-Fouriertransform (FFT) and bi-linear interpolation similar to what is done in GALARIO (Tazzari et al. 2018). This technique is implemented using Tensorflow-GPU (Abadi et al. 2015) for a significant speed-up compared with other available methods.

Geometrical Disk Models
In this work we primarily use a simple geometrical model to describe a disk component with a lopsided brightness distribution: the polar Gaussian ring (sometimes referred to in this paper as an arc). This model has been frequently used for model fitting to mm interferometer observations of disks and rings with azimuthal asymmetries (e.g. Pinilla et al. 2015;Hashimoto et al. 2021), but does not capture the full three dimensional geometry of a thick disk . To account for the dominant emission from the central star, a bright central point source is also included analytically in the Fourier domain.
The intensity distribution of the polar Gaussian ring, in the r, θ plane of the disk, is given by a five-parameter formula where r 0 is the central radius of the ring, θ 0 is the position angle (defined East of North, where East is positive RA and North is positive DEC) of the bright Gaussian peak, σ r defines the radial width of the ring, σ θ defines the azimuthal width of the arc, and I 0 is the peak brightness in the ring. Requiring the origin of the ring to be fixed at the central star allows for only two additional fitting parameters: the observed disk inclination, i, and position angle, PA. Under the assumption that the brightness asymmetry is due to scattering from the inclined disk, however, we equate the position angle with the direction of the peak brightness in the ring such that PA= θ 0 . Finally, we note that the scattered light emission from a thick disk may not follow exactly a circular arc in the plane of the disk (e.g., . Thus, for some geometrical fits we allow the origin of the ring to deviate from the central star position. To account for the possibility of multiple rings in LkCa 15, additional polar Gaussian ring components can be included in the fit. Furthermore, it is likely that a (marginally) resolved inner disk exists around the central star (Espaillat et al. 2007a;Mulders et al. 2010;Alencar et al. 2018) and thus in some cases we also allow for a central, inclined Gaussian (∝ e −(x 2 +(y/cos(i)) 2 )/2 ) disk component to model this emission. The Gaussian disk component is allowed to have its own inclination and position angle, as well as a small offset from the central star.

Bayesian Model Fitting
We use Bayesian model fitting of the LkCa 15 visibilities which allows explicit priors to be included and model comparison to be conducted. We assume our set of measurements are independent Gaussian random variables, so the likelihood, or the probability of obtaining a set of n measurements given a model M with parameters θ, is given by where y i is an individual measurement with measured standard deviation σ i and y(x i ) is the prediction generated by M . For the entirety of our analysis we assume flat priors on all model parameters.
To fit models to the data by sampling the posterior probability, dynamic nested sampling (Higson et al. 2018) with dynesty (Speagle 2020) is used. Nested sampling (Skilling 2004(Skilling , 2006) is a powerful method because, unlike with MCMC based methods, an initial guess for model parameters is not used, so the entire allowed parameter space is efficiently mapped in an un-biased manner. Nested sampling also directly calculates the (model) evidence, or marginal likelihood, which can be used for model comparison. The evidence, P (D|M ), is defined as the integral over the entire set of model parameters, Ω θ , of the likelihood, P (D|θ, M ), multiplied by the prior, P (θ|M ), Our derived model parameters are chosen to be the median of the marginalized posterior distributions and our uncertainties are the 16th and 84th percentiles of the posteriors. The images produced by each model are 256×256 pixels with a pixel size of 4.7 mas and a central star brightness fixed at log 10 I s = 12, in the normalized units used throughout this analysis. We use dynamic nested sampling with flat priors, 1500 initial live points, and batches of 500 live points added until dynesty's default convergence criteria are satisfied.

Direct Image Reconstruction
Many effective, well tested algorithms exist for directly reconstructing images from sparse interferometric data [e.g. SQUEEZE (Baron et al. 2010) and MiRA (Thiébaut 2008)]. However, for this work a simple stochastic gradient descent based algorithm was written in Python using Tensorflow (Abadi et al. 2015). This algorithm takes advantage of Tensorflow's GPU compatibility, and uses auto-differentiation (Rall 1986) to calculate derivatives. The code consists of a single 'locally-connected' layer (Gregor & LeCun 2010) that takes an initial 'prior' image as input and multiplies each pixel value by a trainable parameter that is initialized to one. For our analysis, compactness and total variation regularizers were used (see, Piétu et al. 2006). The compactness hyper-parameter value was chosen by inspecting the gradients to determine the optimal value to concentrate the flux within the region of highest sensitivity for our data (R 300 mas). The total variation hyper-parameter was found using the L-curve method (Hansen 1992).
The algorithm uses the Adam optimizer (Kingma & Ba 2017) with a learning rate of 0.1 to reconstruct an image by minimizing the reduced χ 2 of the output image closure phase and squared visibility measurements plus the regularization term.  the combined 2.1 and 2.3 µm data for LkCa 15 taken by VLT/SPHERE-IRDIS by comparing calculated model squared visibilities and closure phases to the data. This analysis was performed with dynamic nested sampling using dynesty (Speagle 2020) to find the most likely model parameters, to approximate the posterior probability distributions, and to calculate the Bayesian evidence for model comparison.
To help avoid issues with potential degeneracy between the radius and inclination of our ring models, the inclination of LkCa 15 is constrained to be between 40 and 60 degrees, as LkCa 15 is known to have a moderately inclined disk on the sky (e.g. Thalmann et al. 2014;Oh et al. 2016). All other parameters are constrained only to be physically realistic. In addition to the geometrical model parameters, two normalization parameters are included in order to scale the calculated squared visibilities and to correct for the fact that we fit a single model with data obtained at two wavelengths.
As a first test, a single polar Gaussian ring (1PG) was fit to the data. The best fit model (top panel in Figure 2) from the dynamic nested sampling, has a similar geometry and orientation to the near infrared LkCa 15 disk emission features (Currie et al. 2019) previously observed via coronagraphy, and probes angular scales well sampled by the SAM observations, ∼100 mas. The single polar Gaussian ring model has a log evidence, ln(P (D|M )) = −884.2 ± 0.2, indicative of a rather poor fit when compared to other tested models (see below) and suggesting that the model is overly simple to accurately reproduce the LkCa 15 observations. An outer disk, beyond the arc seen in the single polar Gaussian ring model, has been observed previously at similar wavelengths to our data (e.g. Thalmann et al. 2014;Oh et al. 2016;Currie et al. 2019). The bright inner edge of this outer disk component resides at an angular scale which the present SAM observations can probe and therefore should contribute a non-negligible amount to our visibility measurements. We included this outer disk using a second polar Gaussian ring (2PG; middle panel in Figure 2), where to avoid a degeneracy between the rings the de-projected central radius of the outer ring is constrained to be larger than 200 mas. Additionally, the position angle and inclination of the outer ring are constrained to be offset from the inner ring by at most ± 10 and 5 degrees, respectively. For this model, the parameters estimated using dynesty give an inner arc comparable to what is seen in the single polar Gaussian ring image, and a reasonably well constrained faint outer arc, similar to what is seen by Thalmann et al. (2014), Oh et al. (2016), and Currie et al. (2019) (see also Section 5.1). The two component polar Gaus-sian ring model is significantly more consistent with the VLT/SPHERE-IRDIS visibility measurements, having ln(P (D|M )) = −577.7 ± 0.2.
The final component added to the modeling was a marginally resolved central, inclined (elliptical) Gaussian disk, to account for emission on scales less than ∼50 mas (< 8 au). The VLT/SPHERE-IRDIS observations provide a unique opportunity to uncover structure at this angular scale. There exists independent observational evidence for an inner disk (∼1 au) in LkCa 15 from spectroscopic observations (Alencar et al. 2018) and one has been hypothesized in spectral energy distribution (SED) studies (e.g. Espaillat et al. 2007a;Mulders et al. 2010;Espaillat et al. 2010). A sub-au disk component also was used in the LkCa 15 radiative transfer modelling by both Thalmann et al. (2014) and Currie et al. (2019). The best fit model image for the two polar Gaussian rings plus a central, inclined, Gaussian disk (2PG+IG) is shown in the bottom panel of Figure 2. By including the central disk an additional large improvement in the quality of the fit is obtained, with an increase in the log evidence to ln(P (D|M )) = −227.7 ± 0.3. The model parameters derived for this three component model by the dynamic nested sampling are shown in Table 1.
For completeness, the 2PG+IG model was also tested while allowing for a single offset from the origin for both of the polar Gaussian components. A nearly identical solution to the non-offset model was found within the uncertainties, with only a small, insignificant, spatial offset required (∆ DEC = 4 +9 −13 mas and ∆ RA = −4 +6 −4 mas), indicating that there is no requirement for these extra parameters. We also explored a small-scale stellar binary model as an alternative to an inner Gaussian disk but obtained a much poorer fit to the data, ln(P (D|M )) = −386.6 ± 0.3.

Direct Image Reconstruction
The direct image reconstruction technique described in § 3.3 was used to further investigate the robustness of the geometrical model fits. Details of the results are provided in Appendix A and here we present an overview of the main results.
From starting with an axisymmetric smooth distribution of emission, after 250 iterations the direct reconstruction technique recovers a feature very similar to the inner arc found by all our geometrical models, 1PG, 2PG, and 2PG+IG. In addition, with this direct reconstruction, marginally extended bright emission remains near the central star, similar to the central Gaussian disk component found in the 2PG+IG geometrical model. Finally, the outer arc found by the 2PG and 2PG+IG a The parameters r 0 and φ 0 define the offset of the central disk in polar coordinates where φ is defined in the same manner as the position angle.
geometrical model fits, at angular distances poorly measured by the SAM observations, is not recovered.
To test the robustness of the direct image result, we followed a similar approach to that used by Sanchez-Bermudez et al. (2021) and directly reconstructed an image using simulated data from the best fit 2PG+IG model instead of the observations themselves. For this test we obtain a very similar image to that discussed in the previous paragraph. Furthermore, an image was directly reconstructed from simulated data of a multiple point source model, consisting of LkCa 15 and three candidate planets. The reconstructed image clearly retains each of the three distinct point sources. In combination, the direct image reconstructions, although utilizing a significantly under-sampled measurement set, yield results similar to our expectations and thus provide additional confidence that the geometrical models presented here represent robust solutions.

Comparing Disk Image Reconstructions
The best fit properties for the three component geometrical model (2PG+IG) including an inner and an outer polar Gaussian ring plus a central, inclined, Gaussian disk are provided in Table 1. The inner Gaussian (central disk) component in our model has a FWHM ∼5 au, the radii of the inner and outer polar Gaussian rings in our model are ∼17 au and ∼45 au, respectively, with an inclination of ∼45 degrees. On scales of the inner and outer rings, our model is mostly consistent with what has been previously quantified at similar wavelengths using direct or coronagraphic imaging by Thalmann et al. (2014Thalmann et al. ( , 2015Thalmann et al. ( , 2016, Oh et al. (2016) and Currie et al. (2019). The compact central Gaussian disk of our three component model is near the smallest angular scales that we expect to be able to observe with our SAM measurements. This central disk corresponds to the innermost component in the models by Thalmann et al. (2014Thalmann et al. ( , 2015 and Currie et al. (2019), which has been inferred from SED fitting (e.g. Mulders et al. 2010). Our inner ring corresponds to the inner ring observed by Thalmann et al. (2015Thalmann et al. ( , 2016, Oh et al. (2016) and Currie et al. (2019). This component resides at an ideal angular scale given our SAM uv-plane coverage (λ/B ∼ 70-260 mas). Our outer ring resides at the same location as the original ring that was found (Espaillat et al. 2007b) and which has been identified by all of the above studies. This ring is near the outer edge of the spatial scales to which our SAM observations are sensitive.
In Figure 3 we plot the 50% brightness contours for each component of our best fit model over the K-band coronagraphic image by Currie et al. (2019) to visually emphasize the agreement discussed above. The inner arc very closely matches the radial location and extent seen in the K-band coronagraphic image. The outer arc, however, does not match exactly the emission seen in the K-band image. Although the position angle of the peak brightness is in good agreement, the radial extent of the model is shifted somewhat outward. This small discrepancy is likely attributable to our decreased sensitivity at larger separations and (relatively) lower contrast. Additionally, the polar Gaussian model used here is extremely simple and therefore does not take into account expected complexity in the scattered light emission structure . In Figure 3 we also plot our model contours over an ALMA image of LkCa 15 by Facchini et al. (2020), which demonstrates that the rings we find are interior to the rings seen with ALMA. Our best fit model was built up iteratively, starting from a single arc, with the evidence improving with each added component, as shown in Figure 4. These log evidence values can be quantitatively compared using the Bayes factor (Kass & Raftery 1995;Jenkins & Peacock 2011), the ratio of evidences, which says that there is strong evidence for one model in favour of another when the difference between the log evidences is greater than 5 (Trotta 2008). The large increases seen in the evidence with each included model component provides a strong indication that such components contribute significantly to the disk signal present in our visibility measurements. This is an important result because it shows that the SAM observations have leverage over an extended range in angular scales, since both the outer arc and compact central disk are near the edge of the scales probed by VLT/SPHERE-IRDIS SAM. Furthermore, this result reveals the need to model each of these components before hunting for faint forming planets.
Our procedure uses a minimum number of free parameters, and thus the derived geometrical values are not necessarily a true reflection of the inherently 3dimensional structure of the LkCa 15 disk components. This consideration needs to be taken into account when comparing our derived parameters with previous determinations from direct imaging studies (Table 2). Due to the (relative) simplicity of our models, it is satisfying to see the strong agreement in position angle and inclination while not surprising that the radius measures have somewhat larger deviations. Unlike the fits in the previous papers, our rings do not have appreciable offsets from the star, and we have only included a single polar Gaussian for each ring, which is unable to capture the full azimuthal extent of the disk.
To derive more accurate physical properties for the LkCa 15 disk it will be necessary to consider its threedimensional structure. The best approach to accomplish this task would be to use radiative transfer models as the input parameters to the visibility measurement fitting. The difficulty with this approach is the computational cost of radiative transfer models. In this work we were able to evaluate millions of models in minutes, whereas individual three-dimensional radiative transfer simulations, in the mostly optically thick regime our observations occupy, can take on the order of minutes each. This motivates the development of accelerated, or approximate, radiative transfer modelling techniques.

Importance of Including Extended Disk
Components when Searching for Planets

Planets versus Rings
To test whether high contrast, small separation, extended structure can be robustly distinguished from point like emission the same reconstruction method as utilized to fit for rings ( § 4.1) was applied to models consisting of multiple point sources. We therefore fit for one (1PT), two (2PT), and three (3PT) planets orbiting LkCa 15, under the assumption of no extended disk emission. In Figure 4 we present the evidence for each of our multiple point source models while in Table 3 we tabulate the best fit parameters. Although the evidence increases with each additional included point source, even the eleven parameter, three companion, model (3PT) has poorer evidence compared to the eight parameter single polar Gaussian ring model (1PG) described above. This clearly demonstrates that the largest contribution to the visibility measurements is a feature consistent with smooth extended arc-like emission, rather than point-like sources.
We also investigated where in the image plane the best fit point sources are located when no extended structure is included in the model. Figure 5 shows the locations of the point source fits plotted over the (zoomed in) three component geometrical disk model (2PG+IG). As expected, the point source locations trace the peak brightness of the inner disk arc. Similarly, the bottom right panel in Figure 5 shows the claimed planet detections from Kraus & Ireland (2011) (black star) and (Sallum  . The evidence (ln(P (D|M ))) for each model discussed in this paper as a function of the number of free parameters for the model. Recall that the evidence takes into account the number of free parameters and thus can be fairly compared across these models. PT = point source companion (with 1 to 3 companions), PG = polar Gaussian ring (with 1 or 2 independent rings), IG = inclined Gaussian. et al. 2015) (blue stars). Our point source fits are consistent with the previously determined candidate planet positions, with the exception that we do not find any point sources coincident with LkCa 15b. This is not surprising as there is only one single epoch H α detection of LkCa 15 b (Sallum et al. 2015), and more recent work (Mendigutía 2019) has argued that LkCa 15's H α emission is more consistent with a disk than a planet. We also recover contrasts of a few hundred (see Table 3), which are mostly consistent with the contrasts reported by Sallum et al. (2015). We emphasize, however, that our point source model fits are significantly poorer than the extended structure model fits and included here as a demonstration of the potential pitfalls in interpreting sparse data sets using models incompatible with the underlying source structure.
Along with the goodness of fit criterion, the timedomain can be used to separate planets from disk emission. Planets should orbit at Keplerian rates around the central star. In contrast, the lopsided brightness observed from an inclined disk should remain fixed in location. To illustrate this difference, blue arrows in the bottom right panel of Figure 5 denote the expected orbital locations of the previously proposed planets at the time of our observations, following the prescription in Currie et al. (2019). As the Keplerian orbits are long, significantly more time between observations would be required to capture observable changes in the viewing geometry. For systems with disks, however, recovering Keplerian motion is likely to be a key indicator of robust planet detection.

Rings and Planets
To explore whether any point sources that were previously obscured by the bright disk emission are present when the disk emission is modelled, a point source component was added to the three component geometrical model (2PG+IG+PT), including inner and outer rings and a central disk. The best fit parameters were again extracted from the marginalized likelihood distributions obtained using dynamic nested sampling. To within the measured uncertainties, identical ring geometries are found as when the planetary point source is not included. Furthermore, the model including a planetary source has a somewhat higher evidence (Figure 4).
We note that this improvement in the fit is real, in that there is strong evidence for the planet model against the no planet model according to the computed natural logarithm of the Bayes factor of ∼16 > 5 (Trotta 2008). Nevertheless, because our simple geometric models are unlikely to capture all of the expected real structure of the system, we can not rule out the possibility that such enhanced disk features are masquerading as planets in these fits. Thus, the (relatively) small increase in evidence should be treated with extreme caution. To further probe the residuals remaining in our data we also tested the effect of including alternative disk components, such as an additional uniform Gaussian ring, and observed similar small increases in the evidence. We note, however, that these various tests provided no consistency in the additional structure identified and thus we conclude that we have reached the limits of our geometrical model analysis.
In Table 3 we show the extracted companion parameters from our best fit three component disk model plus planet (2PG+IG+PT). We note that the recovered companion brightness contrast is large, greater than 1000, and that the uncertainty in the contrast is ∼20%. Furthermore, the marginalized distributions for the companion are multi-modal, so the quoted radius and position angle are ambiguous. Thus, in the table we present also the parameters extracted from the most prominent modes. For illustrative purposes, in Figure 6 the best fit companion contours are plotted over our 2PG+IG model (left panel) along with a corner plot of the planet parameters (right panel).
From the modeling, it is clear that the rings and inner Gaussian account for the bulk of the asymmetric brightness recorded by the visibility measurements. This is in contrast to the models that only include point sources, where in order to fit the observations, the point source contrasts have to be much smaller, ∼250 for the brightest peak, and the uncertainties on the contrasts much lower (see Table 3). The lack of a specific best-fit location for the candidate planet within a three component model (2PG+IG+PT), coupled with the large required brightness contrast and only marginal improvement in the evidence (Figure 4), further emphasizes the need to treat the companion result with extreme caution. While the modeling does not rule out the presence of planets, it does limit their potential near infrared brightness. A broader, multi-band, and multi-epoch, study as well as further analysis of potential sources of disk asymmetry, beyond the polar Gaussian model, will be required to test against potential faint companion models and determine the robustness of these potential planet detections.

CONCLUSIONS
We have analyzed VLT/SPHERE-IRDIS sparse aperture masking interferometry (SAM) data of the transition disk LkCa 15 and fit geometrical models to our 2.1 and 2.3 µm data simultaneously, uncovering two previously identified arc components along with a marginally a Posterior is multi-modal, see Figure 6. b Extracted parameters from the most prominent modes found in the companion parameters of the 2PG+IG+PT fit, with r 1 , PA 1 and Contrast 1 including the parameters of the two modes at PA ∼135 deg. resolved, previously inferred, compact disk surrounding LkCa 15. The robust 32 mas FWHM (5 au) compact disk detection is an excellent demonstration of the power of SAM for imaging at and beyond the diffraction limit of a telescope. Furthermore, when using SAM, we demonstrate the importance of properly fitting extended structure before searching for candidate planets within transition disks. Our investigation has increased the lower limit contrast for candidate planets in the LkCa 15 system by greater than a factor of three. We emphasize that for our model fitting approach, the advantage of geometrical models over radiative transfer models is that we can evaluate millions of models in minutes, whereas each radiative transfer model can take minutes to compute. Using radiative transfer informed models could be a next step to analyzing more realistic sources of asymmetry in transition disks, potentially using an emulator or accelerating current options using modern machine learning frameworks such as JAX (Bradbury et al. 2018).
The main findings of our LkCa 15 investigation are as follows: • Near infrared emission from the LkCa 15 transition disk system is dominated by a series of arcs, spanning the spatial scales 5 -50 au recoverable via VLT SAM ( Figure 2 and Table 1). These features, the larger of which were previously identified by near infrared direct imaging (Thalmann et al. 2014(Thalmann et al. , 2015(Thalmann et al. , 2016Oh et al. 2016;Currie et al. 2019), lie inside a similar series of rings on scales 50 -100 au revealed in the mm using ALMA by Facchini et al. (2020) (Figure 3); • A marginally resolved compact central Gaussian disk with FWHM ∼5 au and slightly offset from LkCa 15 is robustly detected by the VLT SAM observations. This component, which has previously been inferred (Espaillat et al. 2007a;Mulders et al. 2010;Espaillat et al. 2010;Alencar et al. 2018) as belonging to LkCa 15, demonstrates the power of SAM at small angular scales; • When companions are fit to the VLT SAM observations without first considering the extended disk structure, candidate planets are uncovered (Figure 5) with similar properties to those claimed previously (Kraus & Ireland 2011;Sallum et al. 2015). However, the companion fits are significantly poorer (lower evidence, see Figure 4) than the extended disk structure results leading us to prefer the latter models; • When companions are fit to the VLT SAM observations after including the extended disk structure, no clear evidence for candidate planets is recovered ( Figure 6). The near infrared contrast threshold for planets within 50 au of LkCa 15 is found to be at least 1000 (Table 3). Figure 7. Median reconstructed images using our custom direct image reconstruction method (Section 3.3). Top Left: Image reconstructed using the SAM measurements of LkCa 15, with the white contour denoting the half maximum of the inner arc component of the 2PG + IG geometric model (see also Figure 2). Top Right: Image reconstructed using simulated measurements from the best fit 2PG + IG geometric model. Note the similarity with the top left panel at the arc location and in central asymmetry. Bottom Left: Image reconstructed using simulated measurements from the best fit 2PG + IG geometric model omitting the inner Gaussian component. Note the similarity in the arc location with the top two panels and the dissimilarity in the central structure. Bottom Right: Image reconstructed using LkCa 15 bcd point source model simulated measurements, applying contrasts from Sallum et al. (2015) and predicted locations (denoted by red crosses) as shown in Figure 5. The images in each panel are re-scaled such that the brightest pixel outside 45 mas from the central source is unity.