# Two Rings and a Marginally Resolved, 5 au Disk around LkCa 15 Identified via Near-infrared Sparse Aperture Masking Interferometry

Dori Blakely, Logan Francis, Doug Johnstone, Anthony Soulain, Peter Tuthill, Anthony Cheetham, Joel Sanchez-Bermudez, Anand Sivaramakrishnan, Ruobing Dong, Nienke van Der Marel, et al.

► To cite this version:

Dori Blakely, Logan Francis, Doug Johnstone, Anthony Soulain, Peter Tuthill, et al. Two Rings and a Marginally Resolved, 5 au Disk around LkCa 15 Identified via Near-infrared Sparse Aperture Masking Interferometry. The Astrophysical Journal, 2022, 931, 10.3847/1538-4357/ac6586. insu-03690858

HAL Id: insu-03690858

https://insu.hal.science/insu-03690858

Submitted on 8 Jun 2022

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
OPEN ACCESS

<figure><img src="image_2.png" /></figure>

CrossMark

# Two Rings and a Marginally Resolved, 5 au Disk around LkCa 15 Identified via Near-infrared Sparse Aperture Masking Interferometry

Dori Blakely<sup>1,2 &</sup>, Logan Francis<sup>1,2 &</sup>, Doug Johnstone<sup>1,2 &</sup>, Anthony Soulain<sup>3,4</sup>, Peter Tuthill<sup>3</sup>, Anthony Cheetham<sup>3 &</sup>, Joel Sanchez-Bermudez<sup>5,6 &</sup>, Anand Sivaramakrishnan<sup>7,8 &</sup>, Ruobing Dong<sup>1 &</sup>, Nienke van der Marel<sup>1,9 &</sup>, Rachel Cooper<sup>7 &</sup>, Arthur Vigan<sup>10 &</sup>, and Faustine Cantalloube<sup>10 &</sup>

*Received 2022 January 31; revised 2022 April 6; accepted 2022 April 6; published 2022 May 18*

## Abstract

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 analyses of SAM observations of young stars surrounded by dusty disks have had difficulties disentangling planet and extended disk emission. We analyze 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 ~17 and ~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 an FWHM of 5 au and a similar inclination and orientation to the outer rings. Finally, we recover no clear evidence for candidate planets. By modeling the extended disk emission, we are able to place a lower limit on the near-infrared companion contrast of at least 1000.

*Unified Astronomy Thesaurus concepts: Protoplanetary disks (1300); Exoplanet detection methods (489); Planet formation (1241)*

## 1. Introduction

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 (Andrews et al. 2011; van der Marel et al. 2016; Huang et al. 2018) 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 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_\odot$, and luminosity 1.3 $L_\odot$

(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 pre-transition disk class of objects by Espaillat et al. (2007) owing to the presence of a significant near-infrared excess over the stellar photosphere in the spectral energy distribution (SED), 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 multiepoch 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 and 19 au) using SAM on the Large Binocular Telescope. They also detected LkCa 15b in the accretion-tracing $H\alpha$ emission using spectral differential imaging with

<figure><img src="image_3.png" /></figure>

Original content from this work may be used under the terms of the [Creative Commons Attribution 4.0 licence](https://creativecommons.org/licenses/by/4.0/). Any further distribution of this work must maintain attribution to the author(s) and the title of the work, journal citation and DOI.

<footnote n="1"><sup>1</sup> Department of Physics and Astronomy, University of Victoria, 3800 Finnerty Road, Elliot Building, Victoria, BC V8P 5C2, Canada</footnote>

<footnote n="2"><sup>2</sup> NRC Herzberg Astronomy and Astrophysics, 5071 West Saanich Road, Victoria, BC V9E 2E7, Canada</footnote>

<footnote n="3"><sup>3</sup> Sydney Institute for Astronomy, School of Physics, University of Sydney, NSW 2006, Australia</footnote>

<footnote n="4"><sup>4</sup> Univ. Grenoble Alpes, CNRS, IPAG, F-38100 Grenoble, France</footnote>

<footnote n="5"><sup>5</sup> Instituto de Astronomía, Universidad Nacional Autónoma de México, Apdo. Postal 70264, Ciudad de México, 04510, México</footnote>

<footnote n="6"><sup>6</sup> Max-Planck-Institut für Astronomie, Königstuhl 17, D-69117 Heidelberg, Germany</footnote>

<footnote n="7"><sup>7</sup> Space Telescope Science Institute, 3700 San Martin Drive, Baltimore, MD 21218, USA</footnote>

<footnote n="8"><sup>8</sup> Johns Hopkins University Department of Physics and Astronomy 3400 North Charles, Baltimore, MD 21218, USA</footnote>

<footnote n="9"><sup>9</sup> Leiden Observatory, Leiden University, 2300 RA, Leiden, The Netherlands</footnote>

<footnote n="10"><sup>10</sup> Aix Marseille Univ, CNRS, CNES, LAM, Marseille, France</footnote>
the Magellan Adaptive Optics System. More recent direct imaging work, however, using SPHERE ZIMPOL and IRDIS (Thalmann et al. 2015, 2016), 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.

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 subapertures. Each subaperture 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) = \int I(x, y)\exp(-i2\pi(ux + vy))dydx, \quad (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 owing 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.

<figure><img src="image_4.png" /><figcaption>Figure 1. LkCa 15 VLT/SPHERE-IRDIS (u, v)-plane coverage. The largest and smallest baselines, B, correspond to spatial scales (λ/B) ~70–260 mas (10–40 au).</figcaption></figure>

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: Section 2 outlines how the data were obtained and processed; Section 3 details the techniques used for image reconstruction, including the geometrical models that were fit; Section 4 outlines how the data were analyzed using these techniques and the results of applying these techniques; Section 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 Section 6.

## 2. 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 2018 January 4, using the seven-hole (21-baseline) VLT/SPHERE-IRDIS (Dohlen et al. 2008; Beuzit et al. 2019) SAM mode (Tuthill et al. 2010; Lacour et al. 2011; Cheetham et al. 2016), over the course of about 2 hr. 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 (north up, 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⁴), before closure phases and visibilities were extracted from the Fourier transforms of the images. Each complex visibility (over 21 baselines) is computed over a subsample of pixels using a weighted average. These pixel weights are retrieved from the absolute positions of the seven 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.<fnref n="11" /> A σ-clipping selection was applied to reject eventual bad frames due to bad seeing conditions or adaptive optics 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 owing 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 color variations between 2.1 and 2.3 µm across the components of LkCa 15 are not significant and therefore that the benefit of increasing (u, v)-plane coverage (see Figure 1) outweighs the additional reconstruction uncertainty obtained by considering the two wavelengths together during image reconstruction.

### 3. SAM Analysis Methods

SAM 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

$$V_{\text{binary}}(u, v) = \alpha + \beta \exp(-i2\pi(ux_0 + vy_0)), \quad (2)$$

where α is the central star brightness, β is the off-center source brightness, and x₀ and y₀ 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, underconstrained 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 Fourier transform (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.

#### 3.1. 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 (Dong et al. 2015). 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

$$I_{\text{GR}}(r, \theta) = I_0 \exp \left( - \frac{(r - r_0)^2}{2\sigma_r^2} - \frac{(\theta - \theta_0)^2}{2\sigma_\theta^2} \right), \quad (3)$$

where r₀ is the central radius of the ring, θ₀ is the position angle (defined east of north, where east is positive R.A. and north is positive decl.) of the bright Gaussian peak, σr defines the radial width of the ring, σθ defines the azimuthal width of the arc, and I₀ 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 = θ₀. 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., Dong et al. 2015). 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. 2007; Mulders et al. 2010; Alencar et al. 2018), and thus in some cases we also allow for a central, inclined Gaussian (∞e⁻(x²+(y/cos(i))²)/²) disk component to model this emission.

<footnote n="11">¹¹ See Soulain et al. (2020) for details.</footnote>
The Gaussian disk component is allowed to have its own inclination and position angle, as well as a small offset from the central star.

## 3.2. 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 that 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 $\theta$, is given by

$$P(D|\theta, M) = \frac{1}{\prod_i^n \sqrt{2\pi\sigma_i^2}} \exp\left(\sum_i^n \frac{-(y_i - y(x_i))^2}{2\sigma_i^2}\right), \quad (4)$$

where $y_i$ is an individual measurement with measured standard deviation $\sigma_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, 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 unbiased 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, $\Omega_\theta$, of the likelihood, $P(D|\theta, M)$, multiplied by the prior, $P(\theta|M)$,

$$P(D|M) = \int_{\Omega_{\theta}} P(D|\theta, M)P(\theta|M)d\theta. \quad (5)$$

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.

## 3.3. Direct Image Reconstruction

Many effective, well-tested algorithms exist for directly reconstructing images from sparse interferometric data (e.g., SQUEEZE, Baron et al. 2010; 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 hyperparameter 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 \lesssim 300$ mas). The total variation

hyperparameter was found using the L-curve method (Hansen 1992).

The algorithm uses the Adam optimizer (Kingma & Ba 2014) with a learning rate of 0.1 to reconstruct an image by minimizing the reduced $\chi^2$ of the output image closure phase and squared visibility measurements plus the regularization term.

# 4. Analysis of the VLT/SPHERE-IRDIS Data

## 4.1. Geometrical Model Reconstruction

Previous near-infrared observations of LkCa 15 have attributed the observed emission to scattered light from a disk surface (e.g., Thalmann et al. 2015, 2016; Currie et al. 2019), motivating our use of the arc-like polar Gaussian ring geometrical model (Section 3.1). This model and each subsequent multicomponent model were fit to 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°, 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 \pm 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 that the present SAM observations can probe and therefore should contribute a nonnegligible amount to our visibility measurements. We included this outer disk using a second polar Gaussian ring (2PG; middle panel of Figure 2), where, to avoid a degeneracy between the rings, the deprojected 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 $\pm 10^\circ$ and $5^\circ$, 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
<figure><img src="image_5.png" /><figcaption>**Figure 2.** Best-fit geometrical models based on *dynesty* fits to 2.1 and 2.3 µm VLT/SPHERE-IRDIS closure phase and squared visibility data of LkCa 15 (images are shown on the same fixed intensity scale with the normalization such that the peak polar Gaussian ring component brightness in the bottom panel is one). Top: single polar Gaussian ring model (1PG). Middle: two polar Gaussian rings model (2PG). Bottom: two polar Gaussian rings plus central, inclined, Gaussian model (2PG+IG), with the inner Gaussian component shown by its half-maximum and 1/10 maximum contours (same color scale but normalized by its peak brightness).</figcaption></figure>

et al. (2014), Oh et al. (2016), and Currie et al. (2019; see also Section 5.1). The two-component polar Gaussian ring model is significantly more consistent with the VLT/SPHERE-IRDIS visibility measurements, having $\ln(P(D|M)) = -577.7 \pm 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 SED studies (e.g., Espaillat et al. 2007, 2010; Mulders et al. 2010). A sub-au disk component also was used in the LkCa 15 radiative transfer modeling 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 \pm 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 nonoffset model was found within the uncertainties, with only a small, insignificant, spatial offset required ($\Delta$decl. $= 4_{-13}^{+9}$ mas and $\Delta$R.A. $= -4_{-4}^{+6}$ 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 \pm 0.3$.

## 4.2. Direct Image Reconstruction

The direct image reconstruction technique described in Section 3.3 was used to further investigate the robustness of the geometrical model fits. Details of the results are provided in Appendix, 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 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 undersampled measurement set, yield results similar to our
<table><caption>**Table 1**
Two Polar Gaussian Rings plus Central, Inclined, Gaussian Disk Model (2PG+IG) Parameters Estimated from Dynamic Nested Sampling Marginalized Distributions</caption><thead><tr><th></th><th>log<sub>10</sub>I<sub>0</sub><br>(pixel<sup>-1</sup>)</th><th>r<sub>0</sub><br>(mas)</th><th>i<br>(deg)</th><th>PA<br>(deg)</th><th>FWHM<sub>r</sub><br>(mas)</th><th>FWHM<sub>θ</sub><br>(deg)</th></tr></thead><tbody><tr><td>Inner arc</td><td>7.78<sup>+0.02</sup><sub>-0.02</sub></td><td>107<sup>+3</sup><sub>-3</sub></td><td>44<sup>+2</sup><sub>-2</sub></td><td>-24<sup>+1</sup><sub>-1</sub></td><td>87<sup>+6</sup><sub>-6</sub></td><td>98<sup>+6</sup><sub>-6</sub></td></tr><tr><td>Outer arc</td><td>7.34<sup>+0.04</sup><sub>-0.05</sub></td><td>285<sup>+11</sup><sub>-9</sub></td><td>-2.2<sup>+2.4</sup><sub>-1.8</sub></td><td>-5.9<sup>+1.2</sup><sub>-1.4</sub></td><td>104<sup>+9</sup><sub>-9</sub></td><td>70<sup>+8</sup><sub>-8</sub></td></tr><tr><td>Central disk</td><td>9.45<sup>+0.24</sup><sub>-0.22</sub></td><td>32<sup>+5</sup><sub>-5</sub></td><td>46<sup>+2</sup><sub>-2</sub></td><td>-34<sup>+3</sup><sub>-3</sub></td><td>1.2<sup>+0.2</sup><sub>-0.2</sub></td><td>239<sup>+20</sup><sub>-19</sub></td></tr></tbody></table>

**Notes.** Upper and lower uncertainties are the 84th and 16th percentiles, respectively. The central point-source intensity was added as a constant in the Fourier domain and was fixed at $log_{10} I_s = 12$.

<sup>a</sup> The parameters $r_0$ and $\phi_0$ define the offset of the central disk in polar coordinates, where $\phi$ is defined in the same manner as the position angle.

<figure><img src="image_6.png" /><figcaption>**Figure 3.** Left: ring locations (solid lines) and half-maximum/quarter-maximum brightness contours (dashed lines) for the best-fit two polar Gaussian rings plus central, inclined, Gaussian disk model (2PG+IG) plotted over the K-band coronagraphic image of LkCa 15 by Currie et al. (2019). The white dashed circle denotes a radius of 60 mas, approximately the inner working angle of the CHARIS instrument (Groff et al. 2013). Right: the same model plotted over the mm wavelength ALMA image of LkCa 15 by Facchini et al. (2020).</figcaption></figure>

expectations and thus provide additional confidence that the geometrical models presented here represent robust solutions.

# 5. Discussion

## 5.1. 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 an FWHM of ~5 au, and the radii of the inner and outer polar Gaussian rings in our model are ~17 and ~45 au, respectively, with an inclination of ~45°. 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. (2014, 2015, 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. (2014, 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. (2015, 2016), Oh et al. (2016), and Currie et al. (2019). This component resides at an ideal angular scale given our SAM (u, v)-plane coverage ($\lambda/B \sim 70-260$ mas). Our outer ring resides at the same location as the original ring that was found (Espaillat et al. 2007) and that 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
<figure><img src="image_7.png" /><figcaption>**Figure 4.** 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 one to three companions); PG = polar Gaussian ring (with one or two independent rings); IG = inclined Gaussian.</figcaption></figure>

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 (Dong et al. 2015). 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 evidence values, which says that there is strong evidence for one model in favor 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 provide 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 three-dimensional 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,

<table><caption>**Table 2**

Derived Disk Parameters from Oh et al. (2016) and Currie et al. (2019)</caption><thead><tr><th>Source</th><th>r<sub>0</sub><br>(au)</th><th>i<br>(deg)</th><th>PA<br>(deg)</th></tr></thead><tbody><tr><td colspan="4">Outer Ring</td></tr><tr><td>Oh et al. (2016)<sup>a</sup></td><td>59.0 ± 1.4</td><td>44 ± 1</td><td>-31 ± 2</td></tr><tr><td>Currie et al. (2019)<sup>b</sup></td><td>55</td><td>50</td><td>-30</td></tr><tr><td>This work</td><td>45<sup>+2</sup><sub>-1</sub></td><td>43 ± 2</td><td>-30 ± 1</td></tr><tr><td colspan="4">Inner Ring</td></tr><tr><td>Oh et al. (2016)<sup>a</sup></td><td>29.8 ± 2.0</td><td>44 ± 2</td><td>-18 ± 2</td></tr><tr><td>Currie et al. (2019)<sup>b</sup></td><td>20</td><td>51.5</td><td>-30</td></tr><tr><td>This work</td><td>17.0 ± 0.5</td><td>44 ± 2</td><td>-24 ± 1</td></tr></tbody></table>

**Notes.** The ellipse fits contain a sometimes considerable offset from the star not shown here.

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 three-dimensional 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 modeling techniques.

<footnote><sup>a</sup> Ellipse fit to middle of ring.</footnote>

<footnote><sup>b</sup> Radiative transfer model (radii are inner edge of each disk component).</footnote>
<table><caption>**Table 3**
Companion Parameters Estimated from Dynamic Nested Sampling Marginalized Likelihood Distributions</caption><thead><tr><th>Model</th><th>r<sub>0</sub> (mas)</th><th>PA<sub>0</sub> (deg)</th><th>Contrast<sub>0</sub></th><th>r<sub>1</sub> (mas)</th><th>PA<sub>1</sub> (deg)</th><th>Contrast<sub>1</sub></th><th>r<sub>2</sub> (mas)</th><th>PA<sub>2</sub> (deg)</th><th>Contrast<sub>2</sub></th></tr></thead><tbody><tr><td>1PT</td><td>76.0<sup>+0.6</sup><sub>-0.6</sub></td><td>340<sup>+1</sup><sub>-1</sub></td><td>248<sup>+8</sup><sub>-7</sub></td><td></td><td></td><td></td><td></td><td></td><td></td></tr><tr><td>2PT</td><td>74.3<sup>+1</sup><sub>-1</sub></td><td>324<sup>+1</sup><sub>-1</sub></td><td>288<sup>+14</sup><sub>-12</sub></td><td>80.2<sup>+0.9</sup><sub>-0.8</sub></td><td>357<sup>+1</sup><sub>-2</sub></td><td>294<sup>+14</sup><sub>-13</sub></td><td></td><td></td><td></td></tr><tr><td>3PT</td><td>74.0<sup>+1.2</sup><sub>-1.2</sub></td><td>315<sup>+3</sup><sub>-3</sub></td><td>418<sup>+42</sup><sub>-33</sub></td><td>77.5<sup>+0.9</sup><sub>-0.9</sub></td><td>348<sup>+2</sup><sub>-2</sub></td><td>254<sup>+13</sup><sub>-10</sub></td><td>90.1<sup>+2.2</sup><sub>-2.5</sub></td><td>21<sup>+3</sup><sub>-2</sub></td><td>583<sup>+71</sup><sub>-58</sub></td></tr><tr><td>2PG+IG+PT<sup>a</sup></td><td>144<sup>+10</sup><sub>-14</sub></td><td>78<sup>+142</sup><sub>-2</sub></td><td>1292<sup>+330</sup><sub>-229</sub></td><td></td><td></td><td></td><td></td><td></td><td></td></tr><tr><td>2PG+IG+PT<sup>b</sup></td><td>143<sup>+7</sup><sub>-6</sub></td><td>76.3<sup>+1.4</sup><sub>-1.3</sub></td><td>1214<sup>+316</sup><sub>-215</sub></td><td>94<sup>+8</sup><sub>-5</sub>/150<sup>+5</sup><sub>-6</sub></td><td>134<sup>+3</sup><sub>-3</sub>/135<sup>+1</sup><sub>-1</sub></td><td>1421<sup>+410</sup><sub>-295</sub>/1290<sup>+352</sup><sub>-212</sub></td><td>166<sup>+4</sup><sub>-4</sub></td><td>220<sup>+1</sup><sub>-1</sub></td><td>1323<sup>+447</sup><sub>-227</sub></td></tr></tbody></table>

**Notes.** Upper and lower uncertainties are the 84th and 16th percentiles, respectively.

<sup>a</sup> Posterior is multimodal; see Figure 6.

<sup>b</sup> Extracted parameters from the most prominent modes found in the companion parameters of the 2PG+IG+PT fit, with r<sub>1</sub>, PA<sub>1</sub>, and Contrast<sub>1</sub> including the parameters of the two modes at PA ~ 135°.

<figure><img src="image_8.png" /><figcaption>**Figure 5.** Best-fit locations for one (1PT; top left), two (2PT; top right), and three (3PT; bottom left) companions when modeling only for point sources, plotted over the best-fit 2PG + IG model, from *dynesty* fit of 2.1 and 2.3 µm closure phase and squared visibility data of LkCa 15. The bottom right panel shows the 2014 L' detections of LkCa 15 b, c, d from Sallum et al. (2015) as blue stars, the 2009 L' planet detection from Kraus & Ireland (2011) as a black star, and the approximate predicted LkCa 15 b, c, d locations at the time of our observations (blue arrows).</figcaption></figure>

## 5.2. Importance of Including Extended Disk Components When Searching for Planets

### 5.2.1. 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

(Section 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
<figure><img src="image_9.png" /><figcaption>**Figure 6.** Left: the joint posterior density (copper contours) of the R.A. and decl. position of a companion when it was fit for, along with the two polar Gaussian ring plus inclined Gaussian model (background image + inclined Gaussian contours). Note that there are four potential locations for this planet found by the fitting. Right: corner plot showing only the companion components from the model fit. The four most likely localized spots are clearly identifiable in the PA vs. radius subplot.</figcaption></figure>

<figure><img src="image_10.png" /></figure>

11-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 of Figure 5 shows the claimed planet detections from Kraus & Ireland (2011; black star) and Sallum 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 et al. 2018) 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 time domain 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.

### 5.2.2. 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 modeled, 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 cannot 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 multimodal, 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, multiband, and multiepoch 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.

## 6. Conclusions

We have analyzed VLT/SPHERE-IRDIS 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 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:

1. 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, 2015; Oh et al. 2016; Thalmann 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).

2. 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. 2007; 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.

3. 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.

4. 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).

We thank the anonymous referee for the insightful comments and feedback. We thank Thomas Vandal and Benjamin Pope for their helpful comments and suggestions. We acknowledge support from the Australian Research Council (DP 180103408) that funded part of this work (AMICAL). We also thank the ESO Paranal staff for support for conducting the observations reported in this paper. D.J. is supported by NRC Canada and by an NSERC Discovery Grant. J.S.-B. acknowledges the support received from the UNAM PAPIIT project IA 101220 and from the CONACYT “Ciencia de Frontera” project 263975. This work has made use of the SPHERE Data Centre, jointly operated by OSUG/IPAG (Grenoble), PYTHEAS/LAM/CeSAM (Marseille), OCA/Lagrange (Nice), Observatoire de Paris/LESIA (Paris), and Observatoire de Lyon/CRAL, and supported by a grant from Labex OSUG@2020 (Investissements d'avenir - ANR10 LABX56). SPHERE is an instrument designed and built by a consortium consisting of IPAG (Grenoble, France), MPIA (Heidelberg, Germany), LAM (Marseille, France), LESIA (Paris, France), Laboratoire Lagrange (Nice, France), INAF - Observatorio di Padova (Italy), Observatoire de Genève (Switzerland), ETH Zurich (Switzerland), NOVA (Netherlands), ONERA (France), and ASTRON (Netherlands), in collaboration with ESO. SPHERE was funded by ESO, with additional contributions from CNRS (France), MPIA (Germany), INAF (Italy), FINES (Switzerland), and NOVA (Netherlands). SPHERE also received funding from the European Commission Sixth and Seventh Framework Programmes as part of the Optical Infrared Coordination Network for Astronomy (OPTICON) under grant number RII3-Ct-2004-001566 for FP6 (2004–2008), grant number 226604 for FP7 (2009–2012) and grant number 312430 for FP7 (2013–2016).

*Software:* Tensorflow (Abadi et al. 2015), dynasty (Speagle 2020), corner (Foreman-Mackey 2016), Matplotlib (Hunter 2007), AMICAL (Soulain et al. 2020), Astropy (Astropy Collaboration et al. 2018).

## Appendix
### Direct Image Reconstruction

The direct image reconstruction technique described in Section 3.3, along with a smooth starting image, was used to build an image of LkCa 15 (Figure 7) using the squared visibilities and closure phases. The use of this basic technique
<figure><img src="image_11.png" /><figcaption>**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 rescaled such that the brightest pixel outside 45 mas from the central source is unity.</figcaption></figure>

demonstrates that from an initial assumed brightness distribution consisting of a high-contrast, central point source surrounded by an extended symmetric Gaussian brightness distribution, the most likely reconstruction includes an extended arc to the northwest of the origin, along with marginally resolved asymmetric emission surrounding the central point source. Convergence to the arc solution required 250 iterations. The contrasts used for the initial images were randomly generated with a mean star to total Gaussian brightness contrast of roughly 2 and with the spatial extent of the Gaussian being randomly selected from a distribution with mean FWHM ~ 850 mas. This process was repeated for 10 randomly generated initial contrasts and Gaussian sizes for the input initial maps. The final reconstructed image shown in the top left panel of Figure 7 is the median of the results of 10 individual runs of 250 steps each. We note that in order to converge, it is necessary to begin with a high contrast between the central point source and the extended Gaussian.

The same direct reconstruction procedure was repeated to build images using simulated data from our best-fit 2PG + IG

model (top right panel of Figure 7), our best-fit 2PG + IG model omitting the inner compact Gaussian component (bottom left panel of Figure 7), and an LkCa 15 bcd point-source model (bottom right panel of Figure 7). These tests allow us to probe the significance of any particular recovered structure as a function of the input data set. In each case, randomly generated Gaussian noise on the level expected to be present in our data was added to the simulated data before each of the 10 resets. The similarity in arc and central asymmetric structure observed for both the reconstructed image from the SAM measurements and that from the simulated data measurements of the best-fit 2PG + IG model, as well as both of their obvious dissimilarities with the reconstructed image of the point-source simulated data measurements, supports the validity of our geometrical modeling. Additionally, the obvious dissimilarity between the asymmetric, elliptical central emission seen in the reconstructed images for the SAM data and the 2PG + IG model compared to the symmetric central emission seen in the reconstruction for the 2PG + IG model omitting the inner compact Gaussian component further validates the
inclusion of the inner Gaussian component used in our geometrical modeling.

We note that the outer ring is not at all constrained during these direct image reconstructions. Image reconstruction is an ill-posed inverse problem, with, in our case, 65,536 free parameters and only 1736 data measurements. We thus only expect to be able to recover the most prominent features in our data, including the inner arc and, potentially, the central, Gaussian disk. The outer ring is near the outer edge of the spatial scales to which our data are sensitive and appears fainter than the two inner components. We therefore expect that for it to be recovered using direct image reconstruction strong priors would have to be imposed.

## ORCID iDs

Dori Blakely https://orcid.org/0000-0001-9582-4261

Logan Francis https://orcid.org/0000-0001-8822-6327

Doug Johnstone https://orcid.org/0000-0002-6773-459X

Anthony Cheetham https://orcid.org/0000-0003-3943-4044

Joel Sanchez-Bermudez https://orcid.org/0000-0002-9723-0421

Anand Sivaramakrishnan https://orcid.org/0000-0003-1251-4124

Ruobing Dong https://orcid.org/0000-0001-9290-7846

Nienke van der Marel https://orcid.org/0000-0003-2458-9756

Rachel Cooper https://orcid.org/0000-0001-7864-308X

Arthur Vigan https://orcid.org/0000-0002-5902-7828

Faustine Cantalloube https://orcid.org/0000-0002-3968-3780

## References

Abadi, M., Agarwal, A., Barham, P., et al. 2015, TensorFlow: Large-Scale Machine Learning on Heterogeneous Systems, https://www.tensorflow.org/

Alencar, S. H. P., Bouvier, J., Donati, J. F., et al. 2018, A&A, 620, A195

Andrews, S. M., Wilner, D. J., Espaillat, C., et al. 2011, ApJ, 732, 42

Astropy Collaboration, Price-Whelan, A. M., Sipőcz, B. M., et al. 2018, AJ, 156, 123

Baldwin, J. E., Haniff, C. A., Mackay, C. D., & Warner, P. J. 1986, Natur, 320, 595

Baron, F., Monnier, J., & Kloppenborg, B. 2010, Proc. SPIE, 7734, 78

Beuzit, J. L., Vigan, A., Mouillet, D., et al. 2019, A&A, 631, A155

Bradbury, J., Frostig, R., Hawkins, P., et al. 2018, JAX: composable transformations of Python+NumPy programs, 0.2.5, http://github.com/google/jax

Cheetham, A. C., Girard, J., Lacour, S., et al. 2016, Proc. SPIE, 9907, 99072T

Cieza, L. A., Lacour, S., Schreiber, M. R., et al. 2013, ApJL, 762, L12

Currie, T., Marois, C., Cieza, L., et al. 2019, ApJL, 877, L3

Dohlen, K., Langlois, M., Saisse, M., et al. 2008, Proc. SPIE, 7014, 70143L

Dong, R., Zhu, Z., Rafikov, R. R., & Stone, J. M. 2015, ApJL, 809, L5

Dong, R., Zhu, Z., & Whitney, B. 2015, ApJ, 809, 93

Duvert, G., Young, J., & Hummel, C. A. 2017, A&A, 597, A8

Espaillat, C., Calvet, N., D'Alessio, P., et al. 2007, ApJL, 670, L135

Espaillat, C., D'Alessio, P., Hernández, J., et al. 2010, ApJ, 717, 441

Facchini, S., Benisty, M., Bae, J., et al. 2020, A&A, 639, A121

Foreman-Mackey, D. 2016, JOSS, 1, 24

Francis, L., & van der Marel, N. 2020, ApJ, 892, 111

Greenbaum, A. Z., Pueyo, L., Sivaramakrishnan, A., & Lacour, S. 2015, ApJ, 798, 68

Gregor, K., & LeCun, Y. 2010, arXiv:1006.0448

Groff, T. D., Peters, M. A., Kasdin, N. J., et al. 2013, Proc. SPIE, 8864, 88640H

Haniff, C. A., Mackay, C. D., Titterington, D. J., Sivia, D., & Baldwin, J. E. 1987, Natur, 328, 694

Hansen, P. C. 1992, SIAMR, 34, 561

Hashimoto, J., Dong, R., & Muto, T. 2021, AJ, 161, 264

Higson, E., Handley, W., Hobson, M., & Lasenby, A. 2018, Stat. Comput., 29, 891

Huang, J., Andrews, S. M., Dullemond, C. P., et al. 2018, ApJL, 869, L42

Hunter, J. D. 2007, CSE, 9, 90

Jenkins, C. R., & Peacock, J. A. 2011, MNRAS, 413, 2895

Kass, R. E., & Raftery, A. E. 1995, Journal of the American Statistical Association, 90, 773

Kingma, D. P., & Ba, J. 2014, arXiv:1412.6980

Kraus, A. L., & Ireland, M. J. 2011, ApJ, 745, 5

Kraus, S., Ireland, M. J., Sitko, M. L., et al. 2013, ApJ, 768, 80

Lacour, S., Tuthill, P., Amico, P., et al. 2011, A&A, 532, A72

Lloyd, J. P., Martinache, F., Ireland, M. J., et al. 2006, ApJL, 650, L131

Lubow, S. H., & D'Angelo, G. 2006, ApJ, 641, 526

Maire, A.-L., Langlois, M., Dohlen, K., et al. 2016, Proc. SPIE, 9908, 990834

Mendigutía, I., Oudmaijer, R. D., Schneider, P. C., et al. 2018, A&A, 618, L9

Monnier, J. D. 2007, NewAR, 51, 604

Mulders, G. D., Dominik, C., & Min, M. 2010, A&A, 512, A11

Oh, D., Hashimoto, J., Tamura, M., et al. 2016, PASJ, 68, L3

Pavlov, A., Möller-Nilsson, O., Feldt, M., et al. 2008, Proc. SPIE, 7019, 701939

Piétu, V., Dutrey, A., Guilloteau, S., Chapillon, E., & Pety, J. 2006, A&A, 460, L43

Pinilla, P., Benisty, M., & Birnstiel, T. 2012, A&A, 545, A81

Pinilla, P., van der Marel, N., Pérez, L. M., et al. 2015, A&A, 584, A16

Rall, L. B. 1986, Math. Mag., 59, 275

Sallum, S., Follette, K. B., Eisner, J. A., et al. 2015, Natur, 527, 342

Sallum, S., Skemer, A. J., Eisner, J. A., et al. 2019, ApJ, 883, 100

Sanchez-Bermudez, J., Caratti o Garatti, A., Garcia Lopez, R., et al. 2021, A&A, 654, A97

Skilling, J. 2004, in AIP Conf. Ser. 735, Bayesian Inference and Maximum Entropy Methods in Science and Engineering: 24th International Workshop on Bayesian Inference and Maximum Entropy Methods in Science and Engineering, ed. R. Fischer, R. Preuss, & U. V. Toussaint (Melville, NY: AIP), 395

Skilling, J. 2006, BayAn, 1, 833

Soulain, A., Sivaramakrishnan, A., Tuthill, P., et al. 2020, Proc. SPIE, 11446, 1144611

Speagle, J. S. 2020, MNRAS, 493, 3132

Tazzari, M., Beaujean, F., & Testi, L. 2018, MNRAS, 476, 4527

Thalmann, C., Janson, M., Garufi, A., et al. 2016, ApJL, 828, L17

Thalmann, C., Mulders, G. D., Hodapp, K., et al. 2014, A&A, 566, A51

Thalmann, C., Mulders, G. D., Janson, M., et al. 2015, ApJL, 808, L41

Thiégbaut, E. 2008, Proc. SPIE, 7013, 70131I

Trotta, R. 2008, ConPh, 49, 71

Tuthill, P., Lacour, S., Amico, P., et al. 2010, Proc. SPIE, 7735, 77351O

Tuthill, P. G., Monnier, J. D., Danchi, W. C., Wishnow, E. H., & Haniff, C. A. 2000, PASP, 112, 555

van Cittert, P. 1934, Phy, 1, 201

van der Marel, N., van Dishoeck, E. F., Bruderer, S., et al. 2016, A&A, 585, A58

Vigan, A., Moutou, C., Langlois, M., et al. 2010, MNRAS, 407, 71

Zernike, F. 1938, Phy, 5, 785