New Spatially Resolved Imaging of the SR 21 Transition Disk and Constraints on the Small-Grain Disk Geometry

We present new 0.6 - 4 $\mu$m imaging of the SR 21 transition disk from Keck/NIRC2 and Magellan/MagAO. The protoplanetary disk around SR 21 has a large (~ 30 - 40 AU) clearing first inferred from its spectral energy distribution and later detected in sub-millimeter imaging. Both the gas and small dust grains are known to have a different morphology, with an inner truncation in CO at ~7 AU, and micron-sized dust detected within the millimeter clearing. Previous near-infrared imaging could not distinguish between an inner dust disk with a truncation at ~7 AU or one that extended to the sublimation radius. The imaging data presented here require an inner dust disk radius of a few AU, and complex structure such as a warp or spiral. We present a parametric warped disk model that can reproduce the observations. Reconciling the images with the spectral energy distribution gathered from the literature suggests grain growth to ~2 - 5 $\mu$m within the sub-millimeter clearing. The complex disk structure and possible grain growth can be connected to dynamical shaping by a giant-planet mass companion, a scenario supported by previous observational and theoretical studies.


INTRODUCTION
Protoplanetary disks -the circumstellar disks of dust and gas that remain after star formation -are understood to be the sites of planet formation. Highresolution images of protoplanetary disks in the visible, infrared, and millimeter have informed our understanding of this process by revealing dust and gas properties (e.g. Ansdell et al. 2016), constraining the timescale for disk dispersal (∼ few Myr; e.g. Haisch et al. 2001), and suggesting the dynamical influence of young planets (e.g. Pinilla et al. 2015c;Dong et al. 2018). Direct images of young and forming planets will evince the details of how planets form, constraining accretion rates (e.g . Eisner 2015;Zhu 2015), initial entropies (e.g. Spiegel & Bur-rows 2012), and early atmospheric properties (e.g. Fortney et al. 2008). Recent direct images of protoplanet candidates both in the infrared and in accretion tracers such as Hα have allowed us to estimate their masses and accretion rates (e.g. Keppler et al. 2018;Wagner et al. 2018;Sallum et al. 2015a), as well as atmospheric properties (e.g. Müller et al. 2018).
Transition disks, protoplanetary disks with large clearings in dust, are particularly good targets for protoplanet searches. Their solar-system-sized clearings were first identified through their spectral energy distributions (Strom et al. 1989) and were later confirmed in sub-millimeter imaging (e.g. Andrews et al. 2011). Transition disk holes, gaps, and low stellar accretion rates relative to full disks (e.g. Najita et al. 2007) can be explained by planets accreting material that would have otherwise fallen onto the star (e.g. Bryden et al. 1999). Other disk features observed in some of these objects -such as warps, asymmetries, and spirals seen in both the millimeter (e.g. Pérez et al. 2014;Isella et al. 2013;van der Marel et al. 2013) and in infrared scattered light (e.g. Muto et al. 2012;Marino et al. 2015) -may also suggest the gravitational influence of protoplanets. Directly detecting protoplanets in transition disk clearings would inform our understanding of disk dispersal and disk-planet interactions.
Young and forming planets are predicted to have relatively low contrast at infrared wavelengths (e.g. Eisner 2015; Zhu 2015). However, the large distances to star forming regions and transition disks (e.g. Torres et al. 2007) make imaging solar system scales in the infrared extremely difficult. The technique of non-redundant masking (NRM; e.g. Tuthill et al. 2000b) is well suited for resolving transition disk clearings at these wavelengths. NRM turns a filled-aperture into an interferometric array, providing a smaller dark fringe, 1 and superior point spread function characterization compared to a conventional telescope. NRM can thus probe tighter angular separations than filled-aperture imaging systems, even those equipped with high performance coronagraphs (e.g. Guyon et al. 2014). Masking has led to the detections of both companions and circumstellar disks at separations close to or within the diffraction limit (Ireland & Kraus 2008;Biller et al. 2012;Kraus & Ireland 2012;Sallum et al. 2015a;Huélamo et al. 2011;Willson et al. 2019).
We present infrared NRM and visible filled-aperture observations of the transition disk around SR 21. 2 SR 21 is a G3 star in Ophiuchus, at a distance of 138 pc (Gaia Collaboration et al. 2016). 3 Estimates of its stellar mass range from 1 M , based on millimeter observations and spectro-astrometry (e.g. Pontoppidan et al. 2008;Brown et al. 2009), to 2.5 M , based on stellar evolutionary models (e.g. Prato et al. 2003). It is classified as a Class II young stellar object, with an accretion rate of 10 −7.9 M yr −1 (Manara et al. 2014). Its nature as a transition disk was established by Brown et al. (2007), who required an inner dust disk from 0.22-0.39 AU and a dust disk gap from 0.39-16 AU to explain the spectral energy distribution.
Followup millimeter observations have confirmed this dust clearing and further constrained the morphology of the dust disk. A cavity extending to ∼ 28 AU was detected in 880 µm SMA data (Brown et al. 2009). Reanalysis of these observations yielded a cavity radius of ∼ 40 AU (Andrews et al. 2011). ALMA observations at 450 µm and 870 µm led to cavity radii estimates of 34 AU and 41 AU, respectively (Pinilla et al. 2015b). Furthermore, 450 µm ALMA data show a large-scale disk asymmetry that can be fit by a vortex, and residuals of the fits may suggest spiral structure (Pérez et al. 2014).
CO observations of SR21 reveal a different morphology in the gas. Spectroastrometry from VLT/CRIRES revealed a gas disk truncation at ∼ 7 AU (Pontoppidan et al. 2008) that may be caused by a companion. A gas density drop within ∼ 25 AU has more recently been derived using ALMA 13 CO and C 18 O observations at ∼ 335 GHz ). The gas disk truncation and density drop within the millimeter dust clearing are both consistent with predictions for sculpting by planetary mass companions.
Infrared observations reveal small-grain material within the millimeter clearing and also suggest substellar companions. Imaging at 8.8 µm and 11.6 µm yield characteristic mid-IR sizes of 67 and 92 mas, respectively, corresponding to 9 and 13 AU at 138 pc (Eisner et al. 2009). These data suggest the presence of a compact, warm companion -possibly circumplanetary material -at a disk gap edge at approximately 6 AU (Eisner et al. 2009). H band images also evince small grains within the millimeter clearing, but cannot constrain the location of the disk rim (Follette et al. 2013). The best fit disk position angle appears to change with stellocentric radius, suggesting a warp or spiral structure (Follette et al. 2013). Compact, asymmetric mid-IR emission, dust segregation, and disk warps may all be signs of a substellar companion in the SR 21 transition disk.
The near-infrared and visible light observations of SR 21 presented here probe tighter separations than previous direct imaging studies. They resolve the location of the companion posited in Eisner et al. (2009) at wavelengths where it is predicted to be relatively bright. In Sections 2 and 3 we describe the observations and data reduction, respectively. In Section 4 we describe our analysis of the imaging data and spectral energy distribution, and in Sections 5 and 6 we present our conclusions.
We collected these data in pupil-stabilized mode, which allows the sky to rotate relative to the detector throughout the night. This facilitates calibration by keeping quasi-static speckles fixed on the detector, while true astrophysical signals rotate with the sky (e.g. Marois et al. 2006). Since SDI provides its own reference PSF, this dataset consists of a single visit to SR 21 that we also self-calibrate using angular differential imaging (ADI; e.g. Marois et al. 2006). While we were allocated a half night for these observations, clouds limited the total observing time to ∼ 90 minutes, and the data were of low quality (FWHM ∼ 160 mas).

Ks and L Observations
We observed SR 21 at K s (2.2 µm) and L (3.8 µm) using the 6-hole NRM at Magellan / MagAO / Clio2 (e.g. Close et al. 2012;Morzinski et al. 2014) and the 9hole NRM at Keck / NIRC2 (e.g. Tuthill et al. 2000a). The MagAO observations were taken in 2013, and the Keck datasets were taken between 2017 and 2018. The effective field of view of NRM observations is λ / d hole , where d hole is the diameter of a single sub-aperture in the mask. This corresponds to 950 mas, 700 mas, and 400 mas for MagAO L , Keck L , and Keck K s , respectively.
All data were obtained in pupil-stabilized mode, with the sky rotating on the detector throughout the night. We broke our observations up into "visits", each of which was a single pointing composed of n frame frames. For the MagAO datasets, we positioned the first half of each visit on the top of the frame, and the second half on the bottom of the frame; for Keck the dither positions were in the top left and bottom right quadrants of NIRC2.
In addition to SR 21 we observed unresolved stars as PSF calibrators (see Tables 1 and 2), following a pointing pattern: ...cal 1 -target -cal 2 -target.... To limit calibration errors due to differential refraction, we chose calibrators that were close to SR 21 on the sky. We in general also chose calibrators with similar brightness in the wavefront sensing bandpass (R band), to ensure similar adaptive optics performance. However, the April 2013 and April 2017 observations used calibrators that were much brighter than SR 21 at R. This resulted in comparable calibration for the phase observables, but degraded the squared visibility calibration compared to a well-matched R-band calibrator. For the remainder of the observations, we used calibrators that were better matched to SR 21 in the visible.
We initially select our calibrators using JMMC's searchcal (Chelli et al. 2016). This software uses an interferometrically-measured relation between photometry and stellar diameter to estimate angular diameters (Delfosse & Bonneau 2004). All calibrators have small angular diameters relative to Keck's 10-meter maximum baseline. After observing, we vet all calibrators both by checking that their squared visibilities are consistent with 1, and by fitting companion models and reconstructing images. No resolved structure was found for any PSF calibrator.

Hα Imaging
To reduce the VisAO data, we bias and dark subtract all frames before dividing by a flat field and subframing down to a field of view of 1.6" (200 pixels). We then use pyKLIP (Wang et al. 2015), a python implementation of Karhunen-Loève Image Processing (KLIP; e.g. Soummer et al. 2012) to carry out PSF subtraction. Since accreting companions would have lower contrast at Hα compared to the continuum, while forward scattered light would have equal contrast in the two narrowband filters, we process the VisAO observations using KLIP in KLIP uses several parameters to control how the reference PSF library is built. These include n an , the number of annuli in the field of view that are treated independently; n sub , the number of subsections into which each annulus is divided, IW A, the radius inside of which pixels are discarded; φ mov , the number of degrees by which the image must have rotated to be included in the library; and n b , the number of Karhunen-Loève basis vectors subtracted from the science image. Different source morphologies (e.g. circumstellar disks versus planets) will be better preserved for different choices of these parameters. For example, point-like companions are more easily detected with aggressive parameter choices (e.g. large n an , n sub , n b and small φ mov ), while extended sources such as circumstellar disks may appear as a collection of point sources in an aggressive reduction.
We reduce the continuum-subtracted (SDI) observations using aggressive parameter choices (n an = 6, n sub = 4, φ mov = 6 • , n b > 10), and the summed observations using less aggressive parameters optimized for disk imaging (n an = 1, n sub = 1, φ mov = 10 • , n b < 10). The aggressive number of annuli was chosen so that each annulus had a width roughly equal to the stellar FWHM. For both reductions, we vary the IW A and n b ; changing IW A and n b does not change the resulting images significantly. We also mask an annulus between 27 and 42 pixels (216 -336 mas) in all raw images, since this corresponds to MagAO's control radius, where there are short lived speckles that are not easily subtracted (e.g Follette et al. 2017). We smooth the final KLIP-processed image by a Gaussian kernel with the stellar FWHM (∼ 17 pixels).
Since the conditions and PSF quality for this dataset were variable, we also used bootstrapping to estimate mean images and SNR maps for the combined and SDI KLIP reductions. For each reduction, we KLIP processed a large number (> 1000) of datacubes with the same length as the MagAO data, sampling randomly from the entire set of images. We took the mean of the set of KLIP-processed bootstrapped datacubes as a representative image for each reduction, and we used the scatter over the set to generate the SNR map. Figure 1 shows the final KLIP images, the bootstrapped mean images, and the bootstrapped SNR maps. These tests suggest that both datasets are consistent with noise. Since bootstrapping ADI observations effectively decreases the amount of parallactic angle coverage, we checked that we would still have detected injected companions in the bootstrapped images for both reductions.   Figure 1. MagAO Hα + continuum (top), and SDI (bottom) reductions. The left panels show KLIP images for each dataset, the middle panels show mean bootstrapped images for each dataset, and the right panels show SNR maps created from the bootstrapping tests. The grey shading shows regions that have been masked in the reduction: radii less than 7 pixels, and an annulus from 27-42 pixels corresponding to the AO control radius.

K s and L Imaging
We use an updated version of the data reduction pipeline described in  to process the NRM observations. For each visit, we flatten all frames and then perform dark and sky subtraction by subtracting the median of one dither position from each frame in the other dither position. We Fourier transform the calibrated images to calculate complex visibilities that contain the amplitudes and phases corresponding to each baseline in the mask. Due to the observing bandpass and the finite mask hole size, information from each baseline is encoded in an extended region in the Fourier transform. To use information from the entirety of this region, we calculate closure phases using the method described in Monnier (1999), in which we average bispectra for multiple triangles of pixels for each triangle of baselines in the mask. We then average the bispectra over the cube of images in each visit before taking the phase as the mean closure phase. Since closure phases are correlated, we project them into linearlyindependent combinations of phases called kernel phases (e.g. Martinache 2010; Sallum et al. 2015b). We also calculate squared visibilities by summing the power in all pixels corresponding to each mask baseline and averaging over the cube of images.
We calibrate the kernel phases and squared visibilities using observations of unresolved stars. We use a method presented in , in which we fit a polynomial function in time to each calibrator squared visibility and kernel phase. We sample this function at the time of each target observation to assign instrumental kernel phases / visibilities. We then subtract the instrumental kernel phases from the target kernel phases and divide the instrumental visibilities into the target visibilities. We do this using 0 th − 5 th order polynomial functions and use the order that provides the lowest kernel phase / visibility scatter to calculate the final calibrated observations.
All of the NRM datasets are limited by calibration errors; the observed scatter in the calibrated phases and visibilities is larger than the random scatter over each cube of images. To assign more realistic error bars, we assume that the errors for each dataset are uniform and we fit Gaussian distributions to each night of calibrated kernel phases, closure phases, and squared visibilities. We first subtract the mean signal for the night, to ensure that the best fit distributions have zero mean. Figure 2 compares the best fits to the observed, mean-subtracted phases and visibilities; most of the observations are wellfit by Gaussian functions. However, the noisiest datasets are not and as a result they may have poorly-estimated errors (e.g. squared visibilities for June 28, 2018 and July 24, 2018, and closure / kernel phases for April 5, 2013). While the assigned errors for these datasets should thus be treated with caution, they still significantly down-weight the noisiest nights during fitting.

ANALYSIS
Here we present K s and L images of SR 21 and use them to constrain the morphology of material within the millimeter clearing. We reconstruct these images from the observed closure phases and squared visibilities. Since image reconstruction often relies on assumptions about the source morphology, we first fit parametric models to the data. We fit simple, single companion models, which reveal asymmetries in single-epoch observations and can constrain orbital motion over multiepoch datasets. We also fit geometric and radiative transfer disk models to the data, which can account for static asymmetries and resolved (< 1) squared visibility signals (which imply extended structure). We then reconstruct images from the observations in a parametric way, including model components that have been constrained by these simpler tests. Due to their low quality and small amount of sky rotation, we do not include the Hα observations in the companion and disk fitting that follows. However, we do check that our final constraints on SR 21's morphology are consistent with the MagAO data.

Single Companion Models
We fit single companion models to the kernel phases using a grid search. The models consist of a central delta function and a second delta function having a separation s, position angle PA, (measured east of north), and contrast ∆. We sample the grid in single degree intervals for PA, every 2 mas in s (out to 750 mas), and every 0.1 magnitudes (up to 10 mag) in ∆. Since any orbital motion at resolved angular separations (∼ 23 mas at K s and ∼ 40 mas at L ) would be small on the timescale of one night, we fit the June 6 -7, 2017 data together and the June 28 -29, 2018 data together. Table 3 lists the best fit model for each epoch, Figure 3 shows the observed versus best fit model kernel phases, and Figure  4 shows χ 2 slices at the best fit companion contrast for each epoch.
We calculate the χ 2 interval between the best fit companion model and the null model to estimate the significance of the best fit relative to the null model. Since χ 2 intervals can be corrupted by poor error estimation, we also use Monte Carlo methods to estimate the detection significance. We fit a single companion model to many (> 10 5 ) Gaussian noise realizations drawn from the distributions shown in Figure 2. We use these fits to calculate the false positive probability for each epoch in two ways: (1) We find the number of noise fits that have the same separation as the best fit to the data, and take the false positive probability to be the number of noise fits with equal or lower contrast than the fit to the data (e.g. Kraus et al. 2011). (2) We find the distribution of F-statistics, the best fit χ 2 divided by the null model χ 2 , and calculate the fraction of F-statistics for noise that are less than the F-statistic for the data (e.g. Protassov et al. 2002;Sallum et al. 2015b). Figure 5 shows the distributions of fits to simulated kernel phase noise as well as the distributions of F-statistics. Table 3 lists the χ 2 intervals, the corresponding significance levels, and the false positive probabilities.
Since statistically-significant companion model fits exist, we investigate whether they can be explained by an orbiting companion in the plane of the outer disk (Figure 6), for both of SR 21's stellar mass estimates (1 M and 2.5 M ). The best fit single-companion models do not have positions consistent with orbital motion in the disk plane. The June 2017 best fit, which lies southwest of the star, is ∼ 180 • offset from the April 2017 and June 2018 best fits. While companions to the northeast of the star are allowed by all datasets at 1σ, the 1σ allowed regions are inconsistent with orbital motion in the disk plane for both stellar mass estimates (Figure 6). Extended, static brightness distributions (often associated with scattered light from disk material) have been known to cause spurious companion phase signals in single-epoch datasets (e.g. Huélamo et al. 2011;Sallum et al. 2015b;Cheetham et al. 2015), and position angles that change erratically between local χ 2 minima in multi-epoch datasets (e.g. Sallum et al. 2016). The observed asymmetries in the multi-epoch dataset presented here are more consistent with a static, extended structure than with an orbiting companion.
The squared visibilities are also inconsistent with a simple single-companion model. The companion models that reproduce the kernel phases over-predict the squared visibilities; the models are under-resolved compared to the observations. Furthermore, the squared visibilities have consistent values over the entire range of parallactic angles in each dataset, suggesting a more symmetric morphology than a single companion model. Indeed, single companion fits to the squared visibilities alone cannot reproduce them. Figure 7 shows this: to match the lowest visibilities, a binary model would require low contrast. In order to match the uniformity with parallactic angle, however, the brightness distribution must be relatively centro-symmetric. Both the kernel phases and the squared visibilities point toward a relatively static, extended brightness distribution.

Geometric Disk Models
Since we do not observe orbital motion in the single companion fits, and since a simple companion model cannot reproduce the visibilities, we next fit a grid of geometric disk models to the data to estimate the size of the near infrared emission. Assuming a static disk, we fit all epochs at each wavelength simultaneously. Each model consists of a 2-dimensional Gaussian with an axes ratio (r) and major axis position angle (θ; measured east of north), as well as a central δ function accounting for some fraction f of the total image flux (e.g. ). We use a χ 2 interval to derive the parameter uncertainties. Due to their centro-symmetry, these models all predict zero kernel and closure phases. We first fix the axes ratio at 1, ignoring the major axis position angle, and we then let r and θ vary; Table 4 lists the results. The circular Gaussian models prefer similar sizes for the L and K s emission, while the elliptical Gaussian models prefer a smaller axis ratio and larger FWHM for Ks compared to L , with nearly-overlapping parameter error bars. In both cases, the central delta function flux is significantly higher at Ks than at L .

Aligned Models
SR 21 is known to have a small-grain disk that extends within the millimeter clearing (e.g. Eisner et al. 2009;Follette et al. 2013). To investigate whether a smallgrain circumstellar disk aligned with the millimeter disk can reproduce the observations, we fit a grid of radiative   transfer models to the imaging data and to a spectral energy distribution gathered from the literature. We include the same photometry as Follette et al. (2013), from the Spitzer Infrared Array Camera (IRAC) and Multiband Imaging Photometer (MIPS), AKARI, Infrared Astronomical Satellite (IRAS), the Submillimeter Array (SMA), and the James Clerk Maxwell Telescope Submil-limetre Common-User Bolometer Array (SCUBA). We also add Herschel PACS and SPIRE photometry from 70 µm to 500 µm (Ribas et al. 2017).
We use the radiative transfer software, RADMC-3D (Dullemond 2012) to generate synthetic disk images and spectra self-consistently. We use the python library, pdspy  the necessary input files for RADMC-3D (e.g. dust density grids and radiation sources), to call RADMC-3D to do thermal calculations and generate a large grid of model images and spectra, and to read the RADMC-3D output files. For each disk, we input a standard density profile for a flared disk: where gives the disk scale height. The values r and z are the radius and height in cylindrical coordinates, h 0 and ρ 0 represent the scale height and density at some reference radius r 0 (taken to be 1 AU), and α and β are the density and scale height power law indices. The disk surface density is given by: where γ = α − β, and Σ 0 is the surface density at r 0 . Σ 0 and ρ 0 can be found by assuming a total disk mass and integrating the disk density over all space. For each model we include two disk components: (1) a relatively thin large-grain disk, and (2) a puffier small-grain disk. We fix the total disk dust mass to 5 × 10 −5 M (e.g. Andrews et al. 2011). Each disk component has the following parameters: fractional mass f m , inner radius r i , outer radius r o , minimum grain size a min , maximum grain size a max , grain size distribution index p, as well as γ, β, and h 0 . For the large-grain disk, we fix r i,l , r o,l , γ l , β l , and h 0,l to values consistent  Figure 2. The hollow circles show the best fits to the data for each epoch. Each inset panel shows the cumulative histogram of best fit contrasts for noise fits with the same separation as the fit to the data. The vertical line shows the best fit contrast for the data. We calculate the false positive probability for each epoch by finding the fraction of noise fits with equal or lower contrast than the best fit contrast for the data. Bottom: Histograms show the distributions of F-statistics for fits to Gaussian noise realizations drawn from the distributions shown in Figure 2. The vertical lines show the best fit F-statistic for the single companion models listed in Table 3. For each epoch, the fraction of fits to noise that have F-statistics lower than the best fit companion yields a false positive probability.
with the literature (Andrews et al. 2011;Brown et al. 2009;Pérez et al. 2014;van der Marel et al. 2016 We allow all the disk geometry parameters for the small grain disk to vary, and we also vary the minimum and maximum grain sizes for both disk components, requiring that the minimum grain size for the large grain disk is not smaller than the minimum grain size in the small grain disk (see Table 5). We include non-hydrostatic values for the small-grain disk flaring index, β s , to allow for the possibility of a puffed-up inner rim and/or a disk wind without adding additional model components. We fix the disk inclination to i = 15 • , and test two position angle values consistent with CO obser-vations (P A = 14 • and 194 • ; Pontoppidan et al. 2008). 4 Following Andrews et al. (2011), we fix the fractional mass of the small-grain disk to 0.15, and for both disk components we fix the dust grain size distribution index, p, to 3.5. We adopt a similar dust composition as previous studies -65% silicates and 35% graphite (e.g.

Weingartner & Draine 2001).
Since the imaging and SED are different datasets with different error bars, following Sheehan & Eisner (2017), we create a goodness-of-fit metric (X 2 ) that combines the kernel phase, squared visibility, and SED χ 2 values with relative weights: (4) We normalize the χ 2 arrays by their minimum values, to ensure that we can easily explore weights that both up-and down-weight each dataset relative to the others. This arbitrary choice does not affect the final best fit, since the same relative weighting could be found without normalizing any of the χ 2 values by a scalar. We explore a wide range of weights (∼ 0.005 − 10) for each dataset. Figure 8 shows the best fit disk model and its parameters for w SED = 0.5, w KP = 5, and w V 2 = 1. The aligned disk model can roughly reproduce the SED and the squared visibilities, but cannot match the ker-  nel phases. The second panels from the left in Figure 8 show this; there is a slope offset between the 1:1 modelvs-data relation and the plotted points, especially at L . All asymmetries in this set of models are along the wrong position angle compared to the data. The models that most closely match the observations have the following characteristics: (1) an inner disk truncation at ∼ 4 − 7 AU, (2) a high flaring index (β ∼ 2), and thus a large disk scale height at the truncation radius, and (3) minimum grain sizes larger than ∼ 2 − 5 µm. The disk models with inner radii at the sublimation radius over-predict the squared visibilities; the hot dust close to the star creates too much unresolved flux. Models with many small grains and/or optically-thick disk rims can match the imaging data, but over-predict the 10 µm silicate feature in the SED. Disks with large scale heights and flaring indices scatter enough light to match the angular size of the K s and L imaging, as well as the magnitude (but not direction) of the asymmetry, without over-predicting the 10µm SED. However, as shown in Figure 8, the large scale height can lead to an under-prediction of the short wavelength spectral energy distribution. These tests suggest that more complex disk models are required to reproduce the data.

Warped Disk Models
Since the aligned disk models cannot match the imaging, here we explore a set of warped disk models. We allow the inner regions of the small-grain disk to have both a different inclination and position angle from the millimeter disk. We follow a procedure similar to that presented in Casassus et al. (2018), where the small grain disk has an initial inclination and position angle at its inner radius, a linearly changing orientation through intermediate radii, and the same orientation as the millimeter disk at larger radii. We force the small grain disk to have the same orientation as the millimeter disk by a stellocentric radius of 7 AU, in order to be consistent with the gas position angle determined through spectro-astrometry (Pontoppidan et al. 2008).
We first fit a coarse grid of warped disk models to roughly estimate the best initial inclination and position angle. We fix the small-grain inner radius and density structure to be the same as that shown in Figure 8. We vary the initial position angle and inclination from 0 • to 350 • and 0 • to 80 • , respectively, both with grid spacings of 10 • . We force the disk orientation to change linearly to i = 15 • and either P A = 14 • or P A = 194 • between 5 AU and 7 AU. The best fit orientation has an initial position angle of 280 • and initial inclination of 20 • , and a final position angle of P A = 194 • .
We next fit a larger grid of warped disk models informed by the initial fits (see Table 6). We vary the initial position angle and inclination around the best fit values from the coarse grid, and set the disk orientation to change linearly starting at a fraction f t of the disk inner radius (r i,s ). At a disk radius of r 1 = 7 AU, the final position angle and inclination are P A 1 = 194 • and i 1 = 15 • . We once again vary β, γ, and h 0 , and now include a range of values for the total small-grain disk mass. We introduced this free parameter to allow for different optical depths for disks with the same geometries and dust properties. We select the best fit warped disk models using the same weighting scheme as the aligned disk models. Figure 9 shows representative images, squared visibilities and kernel phases, and the spectral energy distribution. The warped disk models can better match the observables, in particular the kernel phases. Like the aligned models, the warped models prefer large grains to reproduce both the imaging data and the small 10 µm feature in the spectral energy distribution. They prefer geometries with large scale heights, even when the small grain disk mass is allowed to vary. To check whether the warped disks also prefer an inner rim at a radius of a few AU, we ran a small grid of warped models with inner radii of 0.07 AU. Like the aligned models, warped models that extend to the sublimation radius over-predict the squared visibilities.
We also checked that the warped disk model is consistent with the MagAO Hα + continuum images, by KLIP processing a datacube of the model convolved with a PSF having FWHM = 17 pixels. The signal in the final KLIP-processed image is lower than the noise in the observed Hα + continuum image. Directly injecting the PSF-convolved disk model into the Hα + continuum data did not significantly alter the KLIP-processed image.

Parametric Imaging
To look for complex structure beyond simple disk models, we reconstruct images with SQUEEZE (Baron et al. 2010), which uses Markov-Chain Monte Carlo methods to sample the image posterior distribution. We run SQUEEZE in parallel-tempering mode, where chains with different acceptance probabilities can exchange information, in order to effectively explore the parameter space. SQUEEZE can include several regularizers and model components; we explore a variety of these options for the image reconstructions. SQUEEZE can reconstruct images while simultaneously fitting parametric model components such as a central delta function or a resolved disk. We reconstruct images including both of these components -a central delta function representing the star, and a central resolved, circularly-symmetric disk to account for the drop in the observed squared visibilities. These components have three free parameters: (1) the fractional image flux in the delta function, (2) the outer radius of the circularly-symmetric disk, and (3)  image flux in the circularly-symmetric disk. SQUEEZE is free to put 0% of the flux in the resolved disk component; including this component does not force the image to be symmetric, but it can constrain the overall size of the structure in the reconstructed image. We tested three separate regularizations: (1) maximum entropy (Frieden 1972), which favors images with minimal configurational information; (2) total variation (Rudin et al. 1992), which prefers images where most gradients are zero (favoring smooth, piecewise images); and (3) the l − 0 norm, (e.g. Högbom 1974), which favors sparsity. For each epoch and regularization type, we explore a grid of regularization hyperparameters, and choose the final value using the "L-curve" approach (e.g. Thiébaut & Young 2017). We plot the regularizer value (f r ) against the reconstructed image reduced χ 2 . This resembles an "L" where images having relatively constant χ 2 with f r are under-regularized, and images having rapidly-changing χ 2 with f r are over-regularized. The transition between under-and over-regularization (the elbow in the curve) gives the optimal hyperparameter value.
Applying these optimal hyperparameter values showed that the choice of model components affected the reconstruction more significantly than the regularization choice. We thus show unregularized Bayesian images as the final reconstructed images in Figure 10. For all L epochs, 65% − 70% of the flux is contained in a central delta function component and 23% − 27% in a symmetric disk component with a diameter of ∼ 150 mas. The remaining flux is distributed asymmetrically, with the brightest peak to the northeast and a fainter peak to the southwest. At K s , roughly the same amount of flux is contained in the compact component, and slightly less (∼ 17%) in a more compact (∼ 70 mas) disk component. The K s images display an asymmetry to the east and northeast, but lack the pronounced asymmetry to the southwest.
We also explored image reconstructions with delta functions alone and resolved disks alone. In these cases, SQUEEZE's ability to reproduce the kernel phases, and the observed asymmetric structure in the images, stayed the same. The algorithm could still reproduce the squared visibilities when a resolved disk (with no delta function) was included; in this case the preferred disk sizes were ∼ 55 mas at L and ∼ 30 mas at Ks. However, SQUEEZE did not converge to an image that matched the squared visibilities when the only model component included was a central delta function. These tests sug-     4.2, 4.3, 4.4, 4.5, 4.6, 4.7, 4.8, 4.9 5.0, 5.1, 5.2, 5.3, 5.4, 5.5, 5.6, 5.7, 5.8, 5.9, 6  gest that the brightness distribution contains some level of resolved, symmetric structure.
Lastly, we check whether any of the radiative transfer disk models can lead to reconstructed images similar to those presented in Figure 10. For each epoch, we sample the radiative transfer models shown in Figures 8  and 9 with the same (u, v) coverage and sky rotation as the real data. We reconstruct images before and after adding noise to the simulated observations so that their distributions match those shown in Figure 2. We fix the size of the symmetric disk component to be the same as the best fit size for the real reconstructed images, and allow the amount of flux in the delta function and symmetric disk components to vary. Figures 11 -13 show the results. The aligned disk model cannot reproduce the reconstructed images: the peak position angle does not coincide with the peak position angle in the data reconstructions, and it is also constant between Ks and L . The warped disk model images roughly match the observations; they have an extended arc to the northeast that, with noise added, can have a peak position angle that changes slightly from epoch to epoch (middle versus bottom row in Figure 12). This model can also reproduce the peak position angle shifts from Ks to L , and provides a better match to the unresolved flux, although it slightly under-predicts the amount of flux in the extended symmetric component. Among the simple disk models explored here, the models that could better reproduce the SQUEEZE symmetric disk component either had a large 10 µm excess in the spectral energy distribution, or were unable to reproduce the asymmetric structure in the reconstructed images.

The SR 21 Small Grain Disk
Previous observations of SR 21 require a puffy smallgrain disk component, in addition to an optically thick small-grain disk, to reproduce the radial polarized intensity profiles at H band (Follette et al. 2013). The data presented here also require a puffy disk component -all of the disk models that can reproduce the imaging and the spectral energy distribution require a large flaring index (β ∼ 2) and have typical scale heights of ∼ 0.05 − 0.2 AU at 1 AU for the small grains. The disk model shown in Figure 9 provides a good match to the radial brightness profile presented in Follette et al. (2013); it has roughly an r −3 profile, with a discontinuity at 36 AU that is significantly smaller than the errors on the H band radial brightness profile. Although we use a simpler, two-component model, our imaging data support the hypothesis that SR 21 has a small-grain "atmosphere" that may represent a disk wind or a remnant envelope.
Previous scattered light observations could not distinguish between a continuous small-grain disk down to separations of ∼ 0.07 AU (the sublimation radius), or a small-grain disk wall at ∼ 7 AU (Follette et al. 2013), where a gas disk truncation has been reported (Pontoppidan et al. 2008). The higher resolution imaging data presented here suggest a small-grain disk truncation at ∼ 4 − 7 AU. This density distribution may be caused by the gravitational influence of a giant planet interior to the dust disk truncation, a scenario that has previously been modeled for SR 21 (e.g. Pinilla et al. 2015a). This scenario would also be consistent with the observed gas disk truncation at ∼ 7 AU (Pontoppidan et al. 2008) as well as the gas density drop within the millimeter clearing (e.g. van der Marel et al. 2016). The large grain sizes required to match the imaging and SED also agree with this scenario, as grain growth can take place in pressure maxima induced by companions.
Reproducing the imaging data presented here requires an inner disk warp or spiral structure in SR 21. The radiative transfer modeling and image reconstruction simulations in Figures 9 and 12 demonstrate this; matching the data required a significant position angle change between the inner and outer regions of the disk. Posi-tion angle changes with semi-major axis have previously been inferred from SR 21's isophotes at larger radii (Follette et al. 2013), and were thought to indicate spiral structure. While a warped inner disk provides a reasonably good fit to the observations, similar results could likely be achieved in radiative transfer simulations of spiral structure. In this case the Ks band reconstructions, which probe tighter angular separations than the L observations, would still exhibit an offset in peak position angle. Disk structures such as warps and spirals have recently been reported in studies of other transition disks (e.g. Follette et al. 2017;Casassus et al. 2018;Monnier et al. 2019), and features like this can be caused by giant planet companions (e.g. Dong et al. 2015).

Testing the Companion Scenario
Previous mid-infrared observations of SR 21 suggested the presence of a companion near the edge of the inner disk clearing (Eisner et al. 2009). This was interpreted as circumplanetary material around a forming sub-stellar object. While the dominant features in the imaging data presented here are well explained by a small-grain disk component, the disk morphology suggests the influence of a massive companion. Furthermore, the disk models explored here cannot reproduce SR 21's near-infrared excess (see Figures 8 and 9) while matching the imaging and without over-predicting the 10 µm silicate feature in the SED. While this may be caused by the assumed dust composition or the simplicity of the disk model, circumplanetary material could account for the excess flux at these wavelengths.
We generate a contrast curve from the Hα SDI imaging data to check whether we can rule out this scenario (see Figure 14). To estimate the contrast we first split the image into several 17-pixel (1 FWHM) wide annuli. For each annulus, we low-pass filter by a Gaussian with the same FWHM as the stellar PSF to smooth out high frequency noise. We then calculate the mean (x) and standard deviation (s) of each annulus, as well as the number of resolution elements (n) within the annulus. Following Mawet et al. (2014), we assume a Student-t distribution with n − 1 degrees of freedom and standard deviation s, and find the flux that would yield a Gaussian 5σ false positive probability (∼ 3 × 10 −7 ). We then multiply this value by 1 + 1/n and add x to calculate the 5σ contrast for a given separation / annulus.
We calibrate this contrast curve by injecting fake planets into the raw Hα images. We inject single companions  Figure 9 with the same noise properties as the data (right). Images are shown north up east left.
at separations equal to each separation in the contrast curve, at 8 different position angles spaced by 45 • . We process each planet injection with identical KLIP parameters as the final SDI image, and measure the recovered planet flux. For each fake planet, the throughput is the ratio of the recovered to input planet fluxes. We take the mean throughput for the set of 8 fake planets as the throughput for a given annulus / separation. We then divide each flux value in the contrast curve by its throughput to calibrate. Figure 14 (left y-axis) shows the resulting contrast curve. We convert this achievable contrast to a detectable planet mass times accretion rate to place limits on any Hα-bright accreting companions. We assume A v ∼ 6.3 for extinction to SR 21 (Andrews et al. 2011;Follette et al. 2013):, and similar extinction toward the star as any hypothetical companions. We also assume the same Hα line luminosity to accretion luminosity scaling relations as those found for young stars (e.g. Rigliaco et al. 2012), and a planet radius of 1.6 R J . The right y-axis of Figure 14 shows the detectable planet masses times accretion rates as a function of separation. We can only rule out accretion rates of ∼ 10 −2 − 10 −3 M 2 J yr −1 . We also check whether the K s and L imaging data can constrain this scenario by exploring a small range of disk plus companion models. We inject companions into the radiative transfer disk model shown in Figure 9, varying companion semi-major axis, and K s and L contrast with respect to the star. We force the companion position to change from epoch to epoch according to a Keplerian orbit in the plane of the outer disk. We thus also let the true anomaly vary.
A wide range of companion contrasts and semi-major axes are indistinguishable from the null (disk-only) The left y-axis shows contrast in magnitudes, and right yaxis converts that contrast to a detectable planet mass times accretion rate, assuming similar extinction to any companion as that for SR 21, and a planet radius of 1.6 RJ.
model at 1σ using a χ 2 interval. At L , contrasts ranging from 1.5 to 6 magnitudes for semi-major axes of 1 -5 AU are consistent with the disk-only model. At K s , companions with contrasts of ∼ 3 − 6 magnitudes are indistinguishable from the null model for separations larger than ∼ 2 − 3 AU, and contrasts of ∼ 0.5 − 3 magnitudes at smaller separations. For comparison, the companion model described in Eisner et al. (2009), which had a temperature of 730 K and a radius of 40 R , had expected contrasts of 1.5 and 3.3 at L and K s . Thus we cannot rule out the companion parameters presented in Eisner et al. (2009). Many disk plus companion models provide significantly better fits to the Keck data than the disk-only model, with the best fit having contrasts of 3.5 and 6.0 at L and K s , respectively. Companions consistent with these contrasts are capable of filling in SR 21's near infrared excess from ∼ 2 − 5 µm. Given the small number of disk and disk plus companion models explored, and the relatively low achievable L contrast within 5 AU, confirming or ruling out the companion scenario may require higher resolution data and more rigorous modeling efforts.

CONCLUSIONS
We presented new, spatially resolved observations of the SR 21 transition disk in the visible and at K s and L . The multi-epoch data are more consistent with for-ward scattering from a circumstellar disk than with an orbiting companion. The reconstructed images reveal complex disk structure beyond a simple inclined rim, such as a warp or spiral. The imaging also suggests an inner disk truncation at a few AU, and requires a puffy small grain component, in agreement with previous scattered light studies at shorter wavelengths. The disk models explored here require relatively large grains ( 2 − 5 µm) to match both the imaging data and the spectral energy distribution.
The radiative transfer modeling and the reconstructed images support the hypothesis that SR 21 may be shaped by a giant planet companion. This agrees with previous mid-infrared observations that suggested the presence of a ∼ 700 K companion, possibly corresponding to circumplanetary material around a forming substellar object. Our data allow for the presence of similar companion contrasts as those reported in Eisner et al. (2009), and disk plus companion models with higher companion contrasts can provide a better fit to the data than disk models alone. Future observations at multiple wavelengths and higher spatial resolution will evince SR 21's disk structure in greater detail, and will place tighter constraints on substellar companions in this system.
S.S. is supported by an NSF Astronomy and Astrophysics Postdoctoral Fellowship under award AST-1701489. J.A.E. acknowledges support from NSF award number 1745406. K.M.M.'s and L.M.C.'s work is supported by the NASA Exoplanets Research Program (XRP) by cooperative agreement NNX16AD44G