High-contrast imaging has opened a new way of exploring and characterizing exoplanets around young stars (e.g. HR8799 bcde; Marois et al. 2008 Marois et al. , 2010 . This is achieved with large ground-based telescopes by correcting wavefronts of the incident light that are distorted by turbulence of the Earth's atmosphere. Adaptive optics (AO: Beckers 1993) provides real-time corrections to the wavefront distortions using a guide star, which results in a point spread function (PSF) close to a diffraction-limited pattern. Wavefront sensing has classically been performed at optical wavelengths, in which case the visible brightness of the guide star is very important for the wavefront sensing and correction. However, for some specific science targets such as M-type stars or extincted young stellar objects (YSOs) that are faint at optical wavelengths while bright in the near-infrared (NIR), wavefront sensing at NIR wavelengths provides better AO performance. Recently, a NIR Pyramid wavefront sensor (PWFS) was installed in the AO system at Keck-II (Bond et al. 2018 (Bond et al. , 2020 in press) as part of the Keck Planet Imager and Characterizer (KPIC: Mawet et al. 2016) . Science operations using the NIR-PWFS with Keck/NIRC2 began in early 2019.
EM* SR 21 (hereafter SR 21) is a YSO in the Ophiuchus star forming region with a distance of 138.4 ± 1.1 pc (Gaia Collaboration et al. 2018) . Herczeg & Hillenbrand (2014) found a spectral type of F7 and an age of 10 Myr assuming 121 pc for its distance. Sallum et al. (2019) inferred a spectral type of G3 for assuming the Gaia-based distance. Its SED shows strong IR excess but lacks emission at ∼ 5µm, which suggests a transitional disk with an inner gap (Brown et al. 2007 ). High-angular resolution studies have shown a large cavity (e.g. van der Marel et al. 2016) , which together with complicated features such as gaps and spirals suggests possible planet formation within the disk (Sallum et al. 2019; Muro-Arena et al. 2020 ). This YSO is faint at optical wavelengths but bright at NIR wavelengths (e.g., V =14.1, H=7.5; Cutri et al. 2003; Rebull et al. 2018) , making it a suitable testbed for high-contrast imaging with the PWFS to better explore its surrounding disk features and/or companions. In this study we present the latest observation of SR 21, which is among the earliest science targets high-contrast imaging with with the NIR-PWFS following Wang et al. (2020) , who observed PDS 70 using Keck/NIRC2 with the NIR-PWFS. Section 2 describes our observations, data reduction, and the results. Data analysis and discussion including the comparison with the archived Keck/NIRC2 data are described in Section 3. Finally our observations and results are summarized in Section 4.
We observed SR 21 on 2020 May 31 UT with Keck/NIRC2+PWFS in L band combined with the vector vortex coronagraph mask (inner working angle: IWA∼100 mas; Serabyn et al. 2017; Mawet et al. 2017 ) and angular differential imaging (ADI: Marois et al. 2006 ). The wavefront sensing was carried out at H band. The total exposure time is Figure 1 . Simulated throughput map around the vortex mask. The off-axis SR 21 data used in this study were taken with an offset of ∼ 0. 28 from the vortex mask, which affected the throughput of the off-axis PSF.
4110 sec (0.3 sec per exposure × 100 coadds × 137 frames) and the parallactic angle change is 56. • 31. We also took an off-axis image offset by ∼ 0. 28 from the vortex mask, which is comprised of 100 coadds each having an integration time of 0.0075 sec. Typical full width at half maximum (FWHM) was measured to be 8.2 pix (or 81.8 mas with a pixel scale of 9.972 mas/pix). We compared this unsaturated PSF with the simulated ideal PSF at L band, where we took into account throughput loss by the vortex mask (see Figure 1) , and we obtained a Strehl ratio of ∼ 87% (approximately equivalent to the root mean square (RMS) of the wavefront errors 200 nm). The current throughput of the NIRC2-PWFS is ∼ 9% and the maximum elevation of SR 21 was ∼ 45 • during our observations, which affected the AO performance compared with the ideal case (Bond et al. 2018 (Bond et al. , 2020 . The data were taken as a part of a shared-risk program using the NIR-PWFS (PI: Dimitri Mawet) which was aimed at science verification with redder targets. In this program we did not observe suitable reference stars for reference-star differential imaging (RDI: Ruane et al. 2019 ) and thus we decided to use other data sets (see Table 1 and Section 2.2.1) that have the most comparable luminosity at H and W 1 bands to SR 21 among targets observed at this epoch as reference PSFs for the RDI reduction.
We first carried out preprocessing, including bad pixel correction, flat fielding, sky subtraction and image registration following Xuan et al. (2018) and Ruane et al. (2019) . We then performed post-processing to subtract the stellar halo. We investigate two approaches of post-processing here: RDI to investigate the faint disk features and inner potential companions, and ADI to investigate outer potential companions. NOTE-1. Mean DIMM seeing at the summit of Mauna Kea on each date. We note that during these observations part of the seeing information was not recorded. It is possible that the NIRC2 data were taken at different seeing conditions. 2. 2MASS photometry (Cutri et al. 2003) . 3. AllWISE photometry (Cutri et al. 2013 ).
For both approaches, we used the pyKLIP package 1 (Wang et al. 2015 ) that adopted Karhunen-Loève Image Projection (KLIP: Soummer et al. 2012; Amara & Quanz 2012) to produce the most likely reference PSF for the target exposure.
As we did not take suitable data sets for the reference PSF, the RDI results can be affected by the quality of other stars in the reference library. In addition to the SR 21 data, we included WaOph 6 and WSB 52 data sets (PI: Charlotte Bond, Dimitri Mawet) in the PSF library (see Table 1 for the observing logs). We decided to use these two stars as reference stars because their W 1 magnitudes are the closest to SR 21 among all targets taken on and around 2020 May 31. These two stars are young stellar objects (YSOs) surrounded by protoplanetary disks, which have been resolved by Huang et al. (2018a,b) . The NIRC2 L -band PSF for these targets may be influenced by the disks, which could impact the RDI reduction at particularly inner separations.
As mentioned above, the reference stars (WaOph 6 and WSB 52) used in this study are not ideal for the RDI reduction of SR 21 because 1) their W 1 magnitudes are not exactly the same, 2) observational conditions may vary from target to target, and 3) they possibly include the protoplanetary disk signals at L band, which can degrade the efficiency and output of the post-processing with RDI. We first conducted frame selection among the PSF library by taking advantage of the mean square error (MSE: Wang et al. 2004; Ruane et al. 2019) . We evaluated MSE scores at separations ≤ 390 mas (∼5 λ/D) and removed the worst 20% of frames in each data set. After the frame selection the field of view (FoV) was cropped into 2. 0 × 2. 0, which is used to make the reference PSF. We performed RDI reduction including all the science data sets and the selected good reference data sets using pyKLIP. Figure 2 shows the RDI-reduced result (Karhunen-Loèvethe number of basis vector; KL=5) in total intensity, which is overlaid with contours from SPHERE/IRDIS H-band coronagraphic (left) and J-band non-coronagraphic (right) ob- The SPHERE data set presents very faint feature of the disk surface of WaOpH 6, and the SPHERE-based SR 21 spiral arm is ∼30 times brighter than the WaOph 6 disk at the same position, indicating that the contamination of the WaOph 6 disk on the NIRC2 PSF of SR 21 should be relatively very small. For WSB 52 there is no public SPHERE data in the VLT archive.
However, our observation did not confirm other features the SPHERE observations reported. The inner area of the processed image has regions with both positive and negative intensity. Part of the positive areas might correspond to the real disk features but they are likely affected by the poor RDI reduction due to the non-ideal reference stars.
We conducted basic ADI reduction using pyKLIP using all of the 137 frames. We divided the NIRC2 FoV into 9 annuli × 4 subsections to produce the most likely reference PSF at each area. We also assumed the real astrophysical signal, if any, moves 1 pix from frame to frame to avoid heavy self-subtraction caused by the ADI reduction. Figure 3 shows the ADI result (KL=20) for the larger FoV (left) and for the zoomed-in image (right) to compare with the RDI result. We did not detect any significant sources within ∼ 3 " . With a smaller number of KL we also see the spiral feature in the ADI-reduced image (Figure 4 ). In the zoomed-in image of Figure 3 , there are two marginal pointlike sources to the Southwest and South that are coincident with the spiral feature and stable among KLs between 10 and 30. However, their signal to noise ratios (SNRs) are not large enough to robustly claim a detection (∼3.5, 3.7, respectively). There is another marginal source outside the spiral feature at Southeast but this is fainter than the sources mentioned above. In this study, we regard them as the residuals of the spiral feature or uncorrected speckle noise and future observation that can explore deeper. contrast at ∼ 0. 2 − 0. 3 will help to identify the nature of these sources especially in the context of planet-disk interactions in the SR 21 system. We also injected fake sources at separations ranging 0. 1-3. 0 and at a variety of position angles to test whether the inject sources are reliably reproduced after the ADI reduction. We altered the contrast levels of the fake sources as a function of radius (7.5 × 10 −4 at r ≤ 0. 3, 1.5 × 10 −4 at 0. 3 < r ≤ 0. 7, and 5.0 × 10 −5 at r > 0. 7) because a bright source can result in large flux loss. Figure 5 shows a result of the injection test. Although part of the disk features may affect the outcome of the ADI reduction, the injection test shows that the point sources can be reliably recovered.
The disk feature detected in L band may be explained by two scenarios: scattering or self-luminosity. If the spiral feature is self-luminous it may be a result of heating a shock around the protoplanet, which would indicate planet-disk interactions (e.g. Lyra et al. 2016) .
We use our RDI result to extract the surface brightness of the spiral feature because RDI causes only modest signal loss (e.g.
15% with the PCA-based RDI reduction of the disk feature; Uyama et al. 2020b ) compared with self-subtraction of ADI, particularly at such small separations. We converted the pixel count into surface brightness by comparing the photometric result of the off-axis PSF (aperture radius=3.5×FWHM) with the Spitzer/IRAC-3.6µm value (1.23±0.01 Jy; Gutermuth et al. 2009 ). We traced the spiral feature and investigated several areas, which measures ∼60, 100, 80, 70 mJy arcsec −2 at PA=180 • , 210 • , 240 • , and 270 • , respectively 2 . Assuming these values reflect thermal emissions from the optically-thick surface of the spiral, the typical temperature of the spiral corresponds to ∼ 200 K. On the other hand, assuming L SR21 = 11L (Francis & van der Marel 2020) and 0. 23 for the typical separation of the spiral, the effective temperature at the spiral is estimated at 80 K and the derived temperature from our RDI result suggests the scattering effect or another heating mechanism (e.g. the shock scenario; Lyra et al. 2016) .
We also calibrated the SPHERE H-band Q φ image following the IRDAP pipeline. We measured that polarimetric intensity at the spiral surface is 35 -40 mJy arcsec −2 along the spine, with an error of ±1.3 mJy arcsec −2 , which is defined as standard deviation in an annular area (at radii between 175 and 225 mas) of the calibrated SPHERE U φ image that is smoothed with a Gaussian kernel (σ = 2 pix). By comparing these values with radiative transfer simulations one can infer the scattering characteristics and the possibility of thermal emission at the spiral surface, but we note that simply comparing our L -band total intensity with the SPHERE H-band polarimetric imaging result leaves uncertainty on polarization ratio. Obtaining total intensity and polarization intensity at the same wavelength is useful to characterize the scattering mechanisms at the disk surface in detail. As mentioned in Section 2, the reference stars we used in this study are not suitable for SR 21 and thus the extracted surface brightness may be biased. Future follow-up observations with more suitable reference stars will help to better investigate the whole disk feature including the spiral in L band.
We convolved the ADI-reduced image by a circular aperture with a radius of FWHM/2 and then calculated the radial noise profile, which is then compared with the unsaturated PSF of SR 21. We calibrated this radial profile by the flux loss ratio at each separation estimated from the injection test result. Figure 6 shows our 5σ contrast limit: we achieved 10 −4 at 0. 4 and 2 × 10 −5 at 1. 0. With the COND03 model (Baraffe et al. 2003) assuming 10 Myr overlaid in Figure 6 , we set a constraint on ∼ 16 M Jup , 10 M Jup , and 5 M Jup at ∼ 40 au, 55 au, and 110 au, respectively. We assumed 5.8 mag at L band, which is different from Muro-Arena et al. (2020) . Their assumption of 6.8 mag is equivalent to 0.48 Jy at 3.8µm and this value is fainter than photometric values at adjacent wavelengths (e.g. ∼ 1.1 − 1.3 Jy at 3.5-4.6µm; Gutermuth et al. 2009; Cutri et al. 2013; Schlafly et al. 2019 ). Our assumption of 5.8 mag is equivalent to the Spitzer/IRAC-3.6 µm value. Our mass limit is consistent with an updated mass limit of the archival data from Figure We also took into account extinction by the interstellar medium considering A V = 6.2 (Herczeg & Hillenbrand 2014) and wavelength dependence of extinction (A λ ∝ λ −1.75 : Draine 1989) , which corresponds to A L = 0.2.
The evolutionary model we adopted in Figure 6 does not include the accreting mechanism. Given that the gas in the SR 21 disk has not been fully depleted at the predicted location the potential protoplanet would likely be actively accreting. So we also compare our contrast limits with a 'circumplanetary-disk' model (Zhu 2015) , with which the L -band luminosity is calculated as a function of the product of planet mass (M ) and accretion rate (MṀ ) and inner radius of the circumplanetary disk (R in ). We follow Ruane et al. (2017) who investigated constraints on the potential accreting planets in the TW Hya disk. Muro-Arena et al. (2020) predicted a 1M Jup protoplanet at 44 au and PA ∼ 11 • and we use our contrast limit to set a constraint on the mass accretion rate onto this potential protoplanet. Figure 7 shows the comparison between our contrast limit at 0. 3 (contrast: 2.6×10 −4 ) and the cirucmplanetary-disk model at given inner radii (1R Jup ≤ R in ≤ 4R Jup ). Even if we assume 1 M Jup for the potential protoplanet as predicted in Muro-Arena et al. (2020), we are not able to solve the degeneracy between the mass accretion rate and the inner radius of the circumplanetary disk because no observational evidence of the circumplanetary disk in the SR 21 disk has been reported. Therefore our detection limit could set an upper limit of the accretion rate toṀ 1.2 × 10 −5 M Jup yr −1 by referring to a R in = 4R Jup case in the circumplanetary-disk model.
These estimations of the mass/accretion-rate limits do not include potential extinction of the SR 21 disk itself. This leaves uncertainty on the mass estimation but this may not be as important because the L -band detection limits are less subject to the extinction than shorter NIR wavelengths (e.g. JHK bands). 2020) did not report any disk features in the archival Keck data. Our detection limit also achieves a better contrast level particularly at inner separations than the archival data (see Figure 6 for the comparison of the detection limits). At separations ≥ 0. 5 the difference of the contrast limits is likely related to difference of the exposure time (4110 sec and 2160 sec for our observation and the previous observation, respectively). We also note that difference of field rotation (56. • 31 and 18. • 19, respectively) may also affect the detection limits at inner separations because a ratio of flux loss made by the ADI reduction can be affected by the rotation angle.
The main difference is that we used the vortex coronagraph and NIR-PWFS while the archival data were taken with the optical wavefront sensing and no coronagraph mask. We took into account the effect of the vortex mask when we measured the Strehl ratio (see Section 2.1). The observing conditions are also different -typical DIMM seeings were 0. 50 and 0. 77 on 2020 May 31 and 2016 April 23 respectively. We measured the Strehl ratio of the archival data to be ∼69%, which is smaller than ∼87% that we measured for our observation. Although seeing was slightly larger on 2016 April 23 when the archival data were taken, the majority of the additional wavefront errors due to the increased seeing should be outside the AO correction radius. The increase of the Strehl ratio by ∼20% is likely produced by the NIR-PWFS suggesting that the NIR-PWFS meaningfully improved the AO performance of Keck/NIRC2 on such a red target.
We have presented the latest Keck/NIRC2 L -band observation of SR 21. The observation was conducted with a NIR Figure 6 . 5σ contrast limits of the ADI (KL=20) result of our data (red) compared with approximate values of the archival data (blue). We adopted L =5.8 mag, which is different from Muro-Arena et al. (2020), for the SR 21 luminosity. For the plot of the previous data we checked the contrast limit presented in Figure (H band) PWFS and the Strehl ratio is measured to be ∼87% (approximately 200 nm for RMS of the wavefront errors). We used the pyKLIP package to post-process the observation utilizing RDI and ADI techniques. In the RDI-reduced image, we confirmed the spiral feature in total intensity at L band, which is consistent with the VLT/SPHERE polarimetric observation. However, we did not observe suitable reference stars for SR 21 but instead used WaOph 6 and WSB 52, which might also include disk signals. Therefore our RDI result may have been affected by these non-ideal reference stars. Our ADI reduction did not detect any convincing companion candidates but we confirmed the spiral feature at small KLs (KL ≤ 5). In addition to the spiral feature at small KLs, we see two point-like sources with SNRs of ∼ 3.5 − 3.7 at KL = 10-30 but they may be residuals of the spiral feature. We calculated a 5σ contrast limit from the ADI result, which achieved 10 −4 at 0. 4 and 2 × 10 −5 at 1. 0. Assuming COND03 model and 10 Myr, we converted the contrast limit into a mass limit -16 M Jup , 10 M Jup , and 5 M Jup at 40 au, 50 au, and 110 au, respectively. Compared with the archival data taken with optical wavefront sensing we achieved the better Strehl ratio and better contrast level, which suggests that the NIR PWFS improved the AO performance of Keck/NIRC2 and will provide more opportunities for red targets. As the spiral feature is affected by the bias of the reduction method and the reference stars used in this study, future observation of this system will provide better RDI results and enable a detailed discussion on the surface brightness distribution of the protoplanetary disk.
Facility: Keck:II (NIRC2) Software: pyKLIP (Wang et al. 2015) , IRDAP (van Holstein et al. 2017 We wish to acknowledge the critical importance of the current and recent Mauna Kea Observatories daycrew, technicians, telescope operators, computer support, and office staff employees, especially during the challenging times presented by the COVID-19 pandemic. Their expertise, ingenuity, and dedication is indispensable to the continued successful operation of these observatories. The authors wish to acknowledge the very significant cultural role and reverence that the summit of Mauna Kea has always had within the indigenous Hawaiian community. We are most fortunate to have the opportunity to conduct observations from this mountain.