Bin Ren
email: ren@caltech.edu

Élodie Choquet

Marshall D Perrin

Dimitri Mawet

Christine H Chen

Julien Milli

John H Debes

Isabel Rebollido

Christopher C Stark

J Brendan Hagan

Dean C Hines

Maxwell A Millar-Blanchaer

Laurent Pueyo

Aki Roberge

Glenn Schneider

Eugene Serabyn

Rémi Soummer

Schuyler G Wolff

Typeset using L A T E X twocolumn style in AASTeX63 A Layered Debris Disk around M Star TWA 7 in Scattered Light

Keywords: Debris disks (363), Coronagraphic imaging (313), Planetary system formation (1257), Orbital motion (1179)

We have obtained Hubble Space Telescope (HST) coronagraphic observations of the circumstellar disk around M star TWA 7 using the STIS instrument in visible light. Together with archival observations including HST/NICMOS using the F160W filter and Very Large Telescope/SPHERE at H-band in polarized light, we investigate the system in scattered light. By studying this nearly face-on system using geometric disk models and Henyey-Greenstein phase functions, we report new discovery of a tertiary ring and a clump. We identify a layered architecture: three rings, a spiral, and an ≈150 au 2 elliptical clump. The most extended ring peaks at 28 au, and the other components are on its outskirts. Our point source detection limit calculations demonstrate the necessity of disk modeling in imaging fainter planets. Morphologically, we witness a clockwise spiral motion, and the motion pattern is consistent with both solid body and local Keplerian; we also observe underdensity regions for the secondary ring that might result from mean motion resonance or moving shadows: both call for re-observations to determine their nature. Comparing multi-instrument observations, we obtain blue STIS-NICMOS color, STIS-SPHERE radial distribution peak difference for the tertiary ring, and high SPHERE-NICMOS polarization fraction; these aspects indicate that TWA 7 could retain small dust particles. By viewing the debris disk around M star TWA 7 at a nearly face-on vantage point, our study allows for the understanding of such disks in scattered light in both system architecture and dust property.

Introduction

In comparison with the ∼20% detection rate of debris disks around nearby FGK stars (e.g., Sibthorpe et al. 2018), the relative faintness of M stars makes it challenging to detect their surrounding debris disks (e.g., Luppe et al. 2020). By far, only a handful of debris disks have been imaged around M stars in scattered light: AU Mic (Kalas et al. 2004), TWA 7 and TWA 25 [START_REF] Choquet | [END_REF], and GSC 07396-00759 (Sissa et al. 2018). With sensitive upcoming instruments, the study of the circumstellar environments around M stars is imminent.

Detection and characterization of debris disks around M stars help study the formation and evolution of these systems. On the one hand, the facts that M stars comprise more than 70% of Galactic stars (e.g., Miller & Scalo 1979;Muench et al. 2002), that more than 70% of M stars are single (e.g., Lada 2006), and that M stars likely host more planets (e.g., Howard et al. 2012), make M stars promising targets for high contrast imaging search of planets (e.g., Montet et al. 2014). On the other hand, planets can interact with debris disks, and leave observational features on the disks that can help trace the existence of planets (e.g., Ozernoy et al. 2000;Lee & Chiang 2016;Sefilian et al. 2021). Scattered light imaging maps the distribution of the smallest dust particles in debris disks. For M stars, on the one hand, stellar radiation pressure is smaller than gravitational force, which allows for the existence of the smallest dust particles (e.g., Arnold et al. 2019) produced from collisional cascade (e.g., Dohnanyi 1969). These small particles can offer a large surface area to help detect and characterize debris disks, and further aid in tracing hidden planets that perturb disk morphology. On the other hand, however, small dust can be removed by mechanisms that primarily affect M stars (e.g., stellar winds: Plavchan et al. 2005;Strubbe & Chiang 2006;Augereau & Beust 2006;Schüppler et al. 2015, coronal mass ejection: Osten et al. 2013), and these mechanisms consequently pose challenges in detecting debris disks around M stars.

Among the four M star debris disks detected in scattered light, only the view of the TWA 7 system is almost face-on (∼10 • ; e.g., [START_REF] Choquet | [END_REF]Olofsson et al. 2018), while the others are almost edge-on. Studying face-on images of M star debris disks offers the best chance in studying disk formation and evolution, since distinctive features from such processes can be distorted at high inclination (Dong et al. 2016). Specifically, by analyzing face-on images, we can directly trace features such as radial distribution, dust segregation, and azimuthal asymmetries that can relate planetdisk interaction (e.g., Lee & Chiang 2016;Chiang & Fung 2017). In this paper, we re-reduce existing scattered light observations of TWA 7 with state-of-the-art methods, and analyze them together with new HST/STIS observations to have a more comprehensive understanding of the debris disk architecture and dust distribution for this nearly face-on M star debris disk in scattered light.

We describe the observations and data reduction procedure in Section 2. In Section 3 we model the detections in all three instruments. We discuss the spatial components of the system in Section 4, and analyze the dust properties in Section 5. We summarize our findings in Section 6.

1.1. TWA 7 TWA 7, an M3.2 star in the TW Hya association (Herczeg & Hillenbrand 2014;Gagné et al. 2017), is a 6.4 +1.0 -1.2 Myr old star (Binks et al. 2020) that has an estimated stellar mass of 0.46 +0.07 -0.10 M (Stassun et al. 2018) at a distance of 34.10 ± 0.03 pc (Gaia Collaboration et al. 2021). It hosts a circumstellar disk with an infrared excess of L IR /L star = 1.7 × 10 -3 (Kral et al. 2017). Spectral energy distribution (SED) analysis of TWA 7 by [START_REF] Riviere-Marichalar | [END_REF] reveals a bimodal distribution of the disk, which suggests two blackbody rings at 38 au and 75 au. Nontheless, by identifying possible background contamination sources, Bayo et al. (2019) is able to fit the SED with one ring centered at 25 au.

With resolved disk images, we can investigate the discrepancies in SED modeling.

Existing studies of TWA 7 in scattered light reveal the system primarily within ∼2 or 70 au. By assembling a large number of stellar point spread function (PSF) images, [START_REF] Choquet | [END_REF] report the first resolved scattered light image of the TWA 7 system using a 1998 HST/NICMOS observation; using Very Large Telescope (VLT)/SPHERE, Olofsson et al. (2018) identify two rings and a spiral in H-band polarized light. In comparison, ALMA observations of the dust emission at 870 µm show that the disk likely extends beyond ∼2 (Bayo et al. 2019), making it possible that existing studies in scattered light did not have enough instrumental sensitivity in probing the exterior faint regions due to the inverse-square law of stellar illumination.

Being a member of the young TW Hya association, TWA 7 could still possibly host a protoplanetary disk as an M3.2 star. By measuring a 10%Hα width of 111.6 km s -1 for TWA 7 using VLT/X-shooter, which is smaller than the classical threshold of 270 km s -1 (White & Basri 2003), [START_REF] Manara | [END_REF] conclude that this system is not accreting. When we situate the 10%Hα width of TWA 7 in the trend that relates stellar accretion in Figure 3 of Natta et al. (2004), we find that TWA 7 is near the lower boundary of that trend. Assuming that trend can be applied to TWA 7, its accretion rate would be Ṁacc ≈ 10 -12 M yr -1 , which is the modeldetermined upper limit for very low mass objects that have no accretion evidence. In addition, the (0.8-80) × 10 -6 M Earth mass of CO gas observed by ALMA is produced through gas release from exocomets (Matrà et al. 2019). Combining these aspects, TWA 7 is thus more likely a debris disk than a protoplanetary disk.

For an M star, although stellar radiation alone cannot efficiently blow dust out (e.g., Arnold et al. 2019), stellar winds can remove the small sub-micron-sized dust particles in debris disks (e.g., Strubbe & Chiang 2006;Augereau & Beust 2006). To estimate the stellar mass loss rate ( Ṁstar ) for TWA 7, we can convert its Swift X-ray luminosity L X = 4.60 × 10 29 erg s -1 (in 0.3-10 keV; Yang et al. 2012) 1 to a stellar surface flux of F X = 7.56 × 10 6 (R star /R ) -2 erg s -1 cm -2 . Even for a conservative estimate of R star /R = 1 (which can actually reach 0.35), its surface flux exceeds the Wood et al. (2005) threshold of 8 × 10 5 erg s -1 cm -2 by nearly one order of magnitude, a threshold beyond which stellar winds suddenly weaken and thus might not contribute to small dust removal. Nevertheless, similar as AU Mic which is an M1 star that hosts an edge-on debris disk with a comparable F X value, it is unclear whether the Wood et al. (2005) relationship can be applied to small stars like TWA 7 to constrain its stellar wind activity Ren et al. (see AU Mic: Strubbe & Chiang 2006). However, even with no constrained stellar mass loss rate, it is possible to infer that rate by comparing disk radial distribution with models (i.e., Section 4.3 of Strubbe & Chiang 2006): the surface density power law index for the tail of a debris disk varies between -2.5 and -1.5 for different levels of stellar wind activity.

More generally, studying the nearly face-on TWA 7 debris disk can be more informative when contrasting it with circumstellar disks around other M stars. On the one hand, being members of the TW Hya association, the M3.2 star TWA 7 and its M0.5 sibling TW Hya (i.e., TWA 1; Herczeg & Hillenbrand 2014) both have nearly face-on structures, yet TW Hya has a fractional infrared excess L IR /L star = 0.25 and it hosts a protoplanetary disk (e.g., [START_REF] Krist | [END_REF]Weinberger et al. 2002;Debes et al. 2017), while TWA 7 is believed to host a debris disk (e.g., Matthews et al. 2007;Matrà et al. 2019) with L IR /L star = 1.7 × 10 -3 (Kral et al. 2017). On the other hand, the M1 star AU Mic hosts an edgeon debris disk (e.g., Kalas et al. 2004;Boccaletti et al. 2015) with L IR /L star = 3.9 × 10 -4 (Kral et al. 2017), which makes it necessary to study the nearly face-on TWA 7 disk for a more complete understanding of the architectures for debris disks orbiting M stars.

Observations and Data Reduction

We list in Table 1 the exposure information of the three datasets in scattered light for this study: HST/STIS, HST/NICMOS, and VLT/SPHERE. In Figure 1, We present the reduced observations, and calculate the corresponding radial profiles2 assuming the disk is face-on.

HST/STIS (2019)

We observed TWA 7 using the STIS coronagraph on UT 2019 February 3 under GO 15218 (PI: É. Choquet) using 3 HST visits, and its PSF reference star CD-35 6480 3 using 1 visit. The pivot wavelength of STIS is 0.58 µm 4 (wavelength range: 0.2-1.1 µm), the pixel size is 0. 05072 [START_REF] Riley | STIS Instrument Handbook[END_REF]. In each visit, to obtain the best angular coverage, we use the mutually nearly perpendicular WEDGEA1.0 and the BAR5 occulters: the former offers an inner working angle (IWA) of 0. 5, the latter 0. 2 (Debes et al. 2019). In each TWA 7 visit, there are 6 of 146.6 s readouts using BAR5, and 2 of 640.0 s readouts or 2 of 700.0 s readouts using WEDGEA1.0. The three visits have a mutual telescope roll of 20 • , and the total exposure time is 6718.8 s. In the CD-35 6480 visit, there are 3 short exposures on BAR5 that form a 3-point dithering pattern to sample the PSF (step: 0.25 STIS pixel, e.g., Debes et al. 2019), with each exposure having 2 of 150.0 s readouts; there is 1 long exposure on WEDGEA1.0 that has 2 of 670.0 s readouts. The total exposure time for the PSF star is 2240.0 s, with no telescope roll.

We remove the stellar PSF from the target observations using the exposures of the PSF star through classical reference differential imaging (cRDI): we minimize the standard deviation of the residuals in the region of the STIS diffraction spikes. We calibrate the images to units of µJy arcsec -2 using the PHOTFLAM headers as in Ren et al. (2019), see Figure 1c for the final image. To estimate the noise map, we first calculate the standard deviation map of the on-detector images after PSF subtraction, then rotate the standard deviation map according to the on-sky telescope roll of each image, and compute the noise map from the square root of the sum of squared rotated standard deviation maps (e.g., Ren et al. 2019).

HST/NICMOS (1998)

We retrieve the archival NICMOS coronagraphic observations of TWA 7 with the F160W filter using the NIC2-CORON aperture on UT 1998 March 26 under GTO/NIC 7226 (PI: E. Becklin) from the Archival Legacy Investigations of Circumstellar Environments (ALICE; PI: R. Soummer; [START_REF] Choquet | Proc. SPIE[END_REF]Hagan et al. 2018) program5 . The central wavelength of F160W is 1.60 µm (wavelength range: 1.4-1.8 µm), the pixel size6 is 0. 07565, the IWA is 0. 3 [START_REF] Viana | NICMOS Instrument Handbook[END_REF]. There are two telescope orientations for this target, and their telescope roll difference is 29. • 9. Each orientation has three readouts: two 191.96 s ones, and a 223.96 s one. The total NICMOS exposure time is 607.88 s.

We reduce the NICMOS data with multi-reference differential imaging (MRDI): we use multiple images that may come from different stars from the ALICE PSF archive to model a target image. Specifically, we use the non-negative matrix factorization (NMF; Ren et al. 2018a) method: we first select 30% of the most correlated ALICE reference images, then calculate 50 ranked NMF components to model the stellar PSF for TWA 7. We choose the NMF method since it has been shown to better recover faint and extended signals in NICMOS observations (e.g., Ren et al. 2019). We calibrate the NICMOS images using the PHOTFNU parameter for F160W7 as in Ren et al. (2019), and present the final image in Figure 1a. In further analysis, we do not adopt the scaling factor in Ren et al. (2018a) to recover the disk, since we have used 50 NMF components that may not satisfy the requirement that the leading NMF component captures the majority of the disk signal (see Section 2.2.3 of Ren et al. 2018a). We instead adopt a forward modeling strategy for the ring-shaped structures in this system (e.g., [START_REF] Choquet | [END_REF]. Specifically, we subtract disk models from the NICMOS observations, then reduce the observations using the NMF components previously constructed for NICMOS data reduction. The best-fit disk model is the one that minimizes the residuals after such a process. In subsequent forward modeling of the disk, we estimate the noise map by first subtracting disk models from the observations, then rotate the reduced images to on-sky orientation, and calculate the standard deviation of the rotated images. Although this noise estimation procedure can over-estimate the noise, we adopt the noise map to efficiently sample the disk parameters in our modeling procedure.

VLT/SPHERE (2017)

We retrieve the archival H-band SPHERE/IRDIS observations of TWA 7 on UT 2017 March 20 in polarized light under European Southern Observatory (ESO) program 198.C-0209(F) (PI: J.-L. Beuzit) from the ESO Science Archive Facility. The central wavelength is 1.625 µm (wavelength range: 1.48-1.77 µm), and the pixel size is 0. 01225 [START_REF] Maire | Proc. SPIE[END_REF]. The observations used the apodized Lyot coronagraph with a mask radius of 0. 0925 (IWA = 0. 1: Carbillet et al. 2011;Guerri et al. 2011) to suppress the starlight. There are 47 individual exposures, each with 2 of 64.00 s integration: the data have been studied in Olofsson et al. (2018), see their Section 2.1 for the observation details. The total integration time is 6016 s, and the total parallactic angle change is 18. • 0.

We reduce the observations using the IRDAP data reduction pipeline (Version 1.2.4, van Holstein et al. 2020) that performs polarimetric differential imaging (PDI) for IRDIS observations. We follow the calibration procedure in the IRDAP log file, and convert the SPHERE image in units of µJy arcsec -2 by multiplying the image by the ratio between the 2MASS star flux in H-band (1.45 Jy, Cutri et al. 2003) 8and its unocculted detector response9 , see Figure 1b. To reduce the impacts from shot noise in subsequent analysis, we follow Olofsson et al. (2018) and convolve the reduced images with a two dimensional Gaussian that has a standard deviation of 2 pixel.

We use the star-polarization-subtracted Q φ image, which represents the polarized light whose scattering direction is parallel or perpendicular to the radial direction and traces the dust distribution, from the IRDAP output files for analysis. For each radial position, we estimate the noise by calculating the standard deviation within a 3 pixel annulus using the starpolarization-subtracted U φ image, which represents the light that is 45 • from the Q φ light directions and does not trace the dust distribution for optically thin debris disks.

Analysis

Detection

We present the detected features from the disks images in this Section, see Section 4.1 for the corresponding measurements through disk modeling. For the components and their possible motion identified with the aid of disk modeling, see Section 4.2. (The data used to create this figure are available in the "anc" folder on arXiv.)

Ring structure

We identify three ring structures for the TWA 7 system. On the one hand, we confirm the findings in previous NIC-MOS [START_REF] Choquet | [END_REF] and SPHERE (Olofsson et al. 2018) studies within 2 in Figure 1: there is one extended ring peaking at ∼0. 8 in both observations, and a secondary ring at ∼1. 5 in SPHERE data. Our STIS observations confirm the existence of both rings, see the dotted transparent ellipses in Figure 1c. On the other hand, we detect a tertiary ring in our STIS data from 2 to 4 that is outside the views analyzed in [START_REF] Choquet | [END_REF] and Olofsson et al. (2018). By adopting larger views in data reduction and reducing the data with state-of-the-art methods (i.e., IRDAP for SPHERE, and NMF for NICMOS), we are able to recover this tertiary ring in both SPHERE and NICMOS observations.

Ren et al.

The surface brightness radial profile for SPHERE is brighter than that for NICMOS within ∼2 in Figure 1d. This stems from the fact that NMF, albeit recovering the ring beyond 2 for the first time with NICMOS, overfits the disk when scaling factors are not adopted (Ren et al. 2018a, their Figures 2 and3). Overfitting is known to occur in PSF removal for the Karhulen-Loève image projection algorithm (KLIP: Soummer et al. 2012;Pueyo 2016), see, e.g., [START_REF] Choquet | [END_REF] for correcting disk overfitting with forward modeling. Similarly, to study the architecture of the TWA 7 system in subsequent analysis, we perform disk modeling to correct for the overfitting without introducing multiple scaling factors for NMF.

Background star

The southeast point source in the STIS image at (∆RA, ∆Dec) ≈ (3. 9, -1. 5) is a background star, which is located at ∼18 kpc with Gaia EDR3 ID 5444751795151523072 (parallax: 0.057 ± 0.186 mas; Gaia Collaboration et al. 2021). The point source in the NIC-MOS image, see also Figure 2 of Lowrance et al. (2005), at (∆RA, ∆Dec) ≈ (1. 5, -1. 9) is the identical background star in STIS.

The apparent motion of the background source between NICMOS and STIS observations originates from the proper motion of both stars. The proper motion of TWA 7 is 24 times that of the background star (Gaia Collaboration et al. 2021), making TWA 7 motion dominate the relative motion. We do not detect the point source in the SPHERE Q φ image, since star light is not polarized.

Modeling setup

To take into account differences such as pixel size and reduction method for the three datasets, we perform disk modeling for the ring structures in the three observations. In order to prevent the domination of noise from individual instruments, we model the three instruments separately. In addition, aiming at establishing the spatial structure of the system for different dust populations seen in the three instruments, we do not attempt to use one single model to reproduce the observations, which have either different observation wavelengths or different observation techniques.

Geometric structure

Assuming the rings are not dynamically connected with each other, we adopt static geometric disk models to explore the distribution of dust in the system. In cylindrical coordinates, the spatial distribution of the scatterers in a ring fol-lows a double-powerlaw along the mid-plane and a Gaussian dispersion along the vertical axis (Augereau et al. 1999),

ρ(r, z) ∝ r r c -2αin + r r c -2αout -1 2 exp - z hr β 2 ,
(1) where β = 1 for non-flared disks, r c is the critical radius, h is the scale height, α in > 0 and α out < 0 are the asymptotic power law indices when r r c and r r c , respectively. The critical radius can be converted to peak density radius, r max , using Equation (1) of Augereau et al. (1999),

r max = - α in + β α out + β (2αin-2αout) -1 r c . (2)
The TWA 7 disk is nearly face-on, which makes constraining the vertical structure in Equation (1) less meaningful. We thus adopt h = 0.04 from the vertical structure study for debris disks in Thébault (2009). The disk has an inclination of θ inc , which is defined as the dihedral angle between the disk mid-plane and the sky. The position angle of the disk is θ PA , which is defined as the position angle of the disk's semimajor axis-the one which is 90 • counterclockwise from the semiminor axis that is closer to Earth-measured from North to East.

Scattering phase function

Scattering phase function describes the intensity of light as a function of scattering angle, i.e., the angle measured from the incident light ray to the outgoing ray. With a ∼10 • inclination for the TWA 7 disk (Olofsson et al. 2018), we are only able to study the scattering phase function of the system in a limited angle range (i.e., from ∼80 • to ∼100 • ).

Motivated by the fact that scattering phase functions can be described with a single parametric function in such a limited angle range (e.g., Hedman & Stark 2015), we adopt the parametric phase function in Henyey & Greenstein (1941) for the scatterers seen in total intensity in STIS and NICMOS,

I tot (θ) = 1 -g 2 4π(1 + g 2 -2g cos θ) 3/2 , ( 3
)
where θ is the scattering angle that is defined as the angle between the incident and emergent rays, g ∈ [-1, 1] is a parameter that ranges from backward scattering when g ≤ 0 to forward scattering when g ≥ 0.

For the SPHERE Q φ image in polarized light, we have a Rayleigh-like polarization fraction term to modify the Henyey-Greenstein total intensity phase function (e.g., Engler et al. 2017;Olofsson et al. 2018),

I pol (θ) = I tot (θ) × 1 -cos 2 θ 1 + cos 2 θ . (4)
Ren et al.

Approach

The Millar-Blanchaer et al. (2015) code can satisfy the above requirements in geometric structure10 and scattering phase function, we thus use it to produce synthetic images for optically thin disks. To scale the overall brightness of an output disk image to match the observed data, we multiply the output by a parameter f flux . To take advantage of parallel computation, we distribute the calculations on a computer cluster using the DebrisDiskFM package (Ren et al. 2019), and explore the parameter space with a Monte Carlo Marko chain (MCMC) approach using the emcee package (Foreman-Mackey et al. 2013).

To obtain the best-fit disk parameters for the observations, we maximize the following log-likelihood function under independent Gaussian distribution,

ln L (Θ | X obs ) = - 1 2 N i=1 X obs,i -X model,i σ obs,i 2 - N i=1 ln σ obs,i - N 2 ln(2π), (5)
where Θ denotes the set of disk parameters (i.e., θ inc , θ PA , α in , α out , r c , g), X is a disk image with N pixels that is either an actual observation (X obs ) or a synthetic model (X model ), σ is the uncertainty map that has the same dimension of X. Aiming at only establishing the structure for the TWA 7 system in this paper, we do not adjust Equation (5) to address correlated noise using a covariance-based log-likelihood function as in Wolff et al. (2017).

The STIS and SPHERE disk images are not expected to suffer from data reduction biases in cRDI or PDI, while the NICMOS image do in MRDI with NMF. We therefore adopt different modeling approaches to retrieve the disk parameters. For the SPHERE and STIS images, we maximize Equation (5) by directly modeling the disk images. For NICMOS, we use forward modeling: we subtract disk models from the original observations, and perform NMF reduction using the previously calculated NMF components in Section 2.2 to minimize the residuals. Such a negative injection process for NICMOS effectively changes Equation ( 5) by substituting X obs with the residual map after data reduction, and X model with zero.

Disk Modeling

We perform disk modeling to establish the architecture of the TWA 7 system and compare the data across different instruments. To minimize the impact from residual physical structures on disk parameters (e.g., the spiral and secondary ring in Olofsson et al. 2018) for a ring, we exclude certain regions in calculating the likelihood function in Equation ( 5).

After iteratively performing disk modeling and inspecting the residuals with different ignored regions, we detect possibly physical structures with pixelwise average signal to noise ratio (S/N) > 1: the background star, two possible spiral arms, and a clump. We thus converge on masking out the S1 (the arm reported in Olofsson et al. 2018) To model the rings in the TWA 7 system, we have attempted modeling the system by simultaneously exploring all the disk parameters in Equation ( 5) for the three ring components in Figure 1, however, combining multi-dimensional MCMC exploration with disk modeling is computationally prohibitive even with a computer cluster. We have also treated the rings as spatially separated components in different annuli from the star, however the rings could not be depicted by the double-powerlaw description in Equation (1). We detail our adopted modeling procedure below, and present the 50th±34th percentiles for the retrieved disk parameters in Table 2.

SPHERE

We start with modeling the SPHERE image that best resolves the first two rings in Olofsson et al. (2018). Given that the SPHERE observation reveals two radial components that have similar intensity in the star-illumination-corrected radial profile of Olofsson et al. (2018), we model the rings sequentially: we perform disk modeling for one ring, then subtract its best-fit model from the data and fit for another ring. To obtain the best-fit parameters for the system, we convolve the models with a two dimensional Gaussian that has a standard deviation of 2 pixels, and compare the convolved models with the convolved Q φ image.

First, we model the inner ring at ∼0. 8 (i.e., "Ring 1") in Olofsson et al. (2018) within a 184 × 184 pixel region centered at the star. To minimize the impact from the residual light that is close to the coronagraph, we mask out a circular region with a radius of 24 pixel centered at the star (3 times the size of the IWA). In addition, we ignore the region that hosts non-disk signals in calculating the log-likelihood (e.g., Wang et al. 2020), thus reducing the potential impact from the S1 spiral arm in Olofsson et al. (2018) on the retrieved disk parameters.

Second, we remove the best-fit model for the inner ring from the SPHERE image, then model the secondary ring at ∼1. 5 (i.e., "Ring 2") within an annulus that is between 110 and 135 pixel (1. 348 to 1. 658). We assume Ring 2 has identical inclination and position angle as Ring 1, and thus only explore the critical radius and radial power law indices in Equation ( 1), and the Henyey-Greenstein g parameter.

Ren et al. (The data used to create this figure are available in the "anc" folder on arXiv.) Third, we model the tertiary ring (i.e., "Ring 3") after removing the models for the first two rings, assuming identical -0.14

• • • • • • θPA a -94. • 3 +0. • 5 -0. • 5 • • • • • • αin 6.34 +0.07
• • • 9.61 +0.15 -0.16 αout • • • d • • • -3.76 +0.02 -0.02 rc (au) • • • d • • • 92.15
NICMOS αin d • • • • • • • • • αout d • • • • • • • • • rc d (au) • • • • • • • • • g 0.87 +0.02 -0.02 • • • 0.96 +0.02 -0.11
Note.

a Retrieved for Ring 1, and fixed for other instruments and rings.

c Modified by Rayleigh scattering in Equation ( 4) for SPHERE data.

b Peak surface density radius calculated using Equation (2).

d Fixed to SPHERE best-fit values.

e Ring 2 is not recovered in modeling STIS and NICMOS due to resolution and data quality.

f Ring 3 parameters in SPHERE are retrieved from binned image (i.e., 4 × 4 pixel to 1 bin).

inclination and position angle as Ring 1 and exploring identical parameters as for Ring 2. We bin the observation data in a 4×4 pixel region to 1 bin to increase calculation efficiency and S/N. In our first attempt, we could recover a southern clump that is evident in the STIS data, we thus mask out that region in the modeling to reduce its impact on the retrieved disk parameters in disk modeling. We present in Figure 2 (The data used to create this figure are available in the "anc" folder on arXiv.) the major axis reported here is shifted 5 • clockwise 11 . This can be explained by the fact that S1 is located at the southwest region, and the S1 tip have contributed to the fitting of Ring 1 in Olofsson et al. (2018). On the other hand, by additionally ignoring Ring 2, the best-fit power law index for the tail of Ring 1 is -1.7 here, which is steeper than the Olofsson et al. (2018) best-fit of -1.5; this has been expected then. We do not discuss other parameters since they are less impacted by S1 and Ring 2. What is more, we do not compare the credible intervals since we have not adopted the kernel density approach, where the kernel widths depend on user inputs, as in Olofsson et al. (2018).

11 Note the difference in the definition, θ PA here is 180 • -φ in Olofsson et al. (2018).

STIS

We follow the above SPHERE procedure of sequentially modeling the ring components to model the system architecture in STIS. To simulate the responses of the STIS detector, we convolve the models with a TinyTim PSF [START_REF] Krist | Proc. SPIE[END_REF] 12 that has an effective temperature of 4018 K for TWA 7 (Gaia Collaboration et al. 2018). See Figure 3 for the detailed sequential modeling results.

After attempting to vary α out for Ring 1, we do not obtain statistically different result from the SPHERE value, we thus fix it to be the best-fit value from SPHERE for STIS. Despite Ring 2 being visible by eye in Figure 1c, we are limited by residual noises in disk modeling, and thus we cannot recover Ring 2 in STIS modeling beyond 1σ in f flux . We therefore do not report the disk parameters for Ring 2.

In comparison with the SPHERE residuals in Figure 2d, we tentatively identify in STIS an additional spiral, S2, with an average pixelwise S/N of 2.3. If S2 does not originate from residual speckles in Figure 3c, then STIS is able to better resolve it than SPHERE. Such a scenario is likely, since that the exposure time with STIS is 11 times that with SPHERE, and that STIS observes the system in total intensity while SPHERE in polarized light.

NICMOS

NICMOS and SPHERE cover nearly identical wavelengths, we thus use the SPHERE spatial distribution to model the NICMOS image. To address the difference between total intensity and polarized phase functions in Equations (3) and ( 4), we only vary the g parameters and brightness levels for the three NICMOS rings. We rotate a disk model according to the telescope roll angles of the NICMOS observations, then convolve the rotated disk models with a TinyTim PSF for the F160W filter. In forward modeling, we subtract them from the NICMOS observations, and reduce the disk-subtracted observations with previously generated NMF components.

We simultaneously model the three rings in NICMOS to minimize potential overfitting biases in forward modeling, see the best-fit results in Figure 3. Limited by residual noises, we cannot recover Ring 2 in NICMOS beyond 1σ in f flux . For the purpose of illustration only, we adopt the best-fit parameters for Ring 1, perform NMF forward modeling, and present the corresponding results in Figure 3b'. In comparison with the retrieved disk parameters in [START_REF] Choquet | [END_REF], where the disk parameters are loosely constrained using only NIC-MOS data, we are able to better constrain the g parameter for Ring 1 using the spatial distribution inferred from the SPHERE data. Nevertheless, given that the ∼13 • inclination of this system only allows the study of the scattering phase function from 77 • to 103 • , we note that a discussion of the Henyey-Greenstein g parameter, given the large NICMOS pixel size and the contaminations on the rings, would be less informative, see, e.g., Section 5.3.

Spatial Distribution

Using static models, our disk modeling reveals a layered architecture for the TWA 7 debris disk system. The extended ring, Ring 1, whose density peaks at 28 au dominates the overall surface density distribution from ∼10 au to ∼200 au. The secondary and tertiary rings, Ring 2 and Ring 3, peak at 52 au and 93 au, are superimposed onto the outskirts of Ring 1, see Section 4.1.1. In comparison with the locations of the TWA 7 rings from SED analysis (e.g. [START_REF] Riviere-Marichalar | [END_REF]Bayo et al. 2019), we find that blackbody or Mie assumption cannot predict the rings in scattered light images with high accuracy (see also Esposito et al. 2020), and correction factors are needed to account for such discrepancies (e.g., sub-mm to blackbody: Pawellek et al. 2021).

There likely exist two spirals, S1 and S2 that are superimposed on Ring 1 and Ring 2, see Section 4.2.1 for the smoothed residuals. S1 has been tentatively detected in Olofsson et al. (2018), and we now well-resolve it in both SPHERE and NICMOS data after disk modeling. S2 is tentatively resolved in STIS residuals, and it has an average pixelwise S/N of 2.3; yet it is less resolved in SPHERE and NICMOS.

Superimposed onto Ring 3, there exists a clump in all three observations, see Section 4.2.3. We resolve this elliptical clump in all three instruments, and our STIS data best resolves the clump with with an area of ≈150 au 2 . By establishing a 21 yr span in this paper, this dusty clump is a component on the TWA 7 debris disk, since a background galaxy would have otherwise moved with respect to TWA 7. In comparison, this clump has been marginally resolved at ∼2σ by ALMA in Bayo et al. (2019). Our modeling shows that there are three rings in the TWA 7 system, with an extended one centered at ∼0. 8 dominating the overall surface brightness distribution as well as the surface density distribution. For the SPHERE data, we multiply the surface brightness for each location with its corresponding squared stellocentric distance to compare the surface density radial profile, see Figure 4.

Although Ring 1 dominates the surface brightness distribution of the system, we observe in Figure 4 that Ring 2 and Ring 3 have similar surface density in the radial profile for "Residual 1", which corresponds to the residual map after removing Ring 1 best-fit model from the SPHERE Q φ image. When Ring 2 and Ring 3 are superimposed onto the tail of Ring 1, the three rings have similar levels of total surface density that peak sequentially at 0. 8, 1. 5, and 2. 8.

We confirm the existence of Ring 1 and Ring 3 in NICMOS and STIS observations with disk modeling, see Figure 3. In comparison, although we are able to resolve Ring 2 by eye in STIS in Figure 1c, the radially extended spike structures between Ring 1 and Ring 2 may have led to overfitting for Ring 1 in STIS. The overfitting then overwhelms the Ring 2 signal in disk modeling. In addition, we are not able to resolve Ring 2 either by eye or from forward modeling in NIC-MOS. However, we observe a plateau at 1. 5 in Figure 1d that is indicative for the existence of Ring 2, yet this plateau may originate from joint effects from both the spirals and Ring 2. Figure 4. Surface density radial profiles, calculated by multiplying surface brightness radial profiles with squared stellocentric distances, for SPHERE data and models. Ring 1 dominates the surface density distribution of the system. Ring 2 and Ring 3 are of similar surface density: when they are superimposed onto the tail of Ring 1, the three ring components have similar levels of total surface density. See Figure 2 for the corresponding surface brightness distribution, and Table 2 for the parameter values in Equation (1). Note: the surface density radial profiles are normalized by dividing the peak value from the observation.

The power law index for the tail of Ring 1 is α out = -1.739 +0.013 -0.013 from SPHERE disk modeling. The index does not have a significant difference between SPHERE and STIS observations. This suggests that the probed dust particles are likely well-mixed for Ring 1 tail in STIS and SPHERE wavelengths. In comparison, Equation (5) of Thébault & Wu (2008) predicts the dominant particle size as a function of radial separation from the star. Nevertheless, the similar power law indices between SPHERE and STIS here suggest that the larger particles, which are probed by SPHERE and are the dominant sizes for closer-in regions, will continue dominate further out regions.

Using the relationship between α out and surface density power law index Γ out in Augereau et al. (1999), i.e., Γ out = α out + β where β = 1 as in Equation (1), the surface density power law index of the Ring 1 tail is Γ out = -0.739 +0.013 -0.013 , see the annotation on Figure 4. This surface density is steeper than the Olofsson et al. (2018) index of -0.52, since we have our excluded Ring 2 in our fitting. Meanwhile, Schneider et al. (2018) and Ren et al. (2019) have observed similar indices of approximately -0.7 for the tails of the debris disks orbiting A0 star HR 4796 A and F5 star HD 191089.

We discuss the meaning of the power law index of Ring 1 tail below assuming Ring 1 is the birth ring of small particles in the system. However, we caution that with a single birth ring, stellar winds around M stars can create multiple rings in shorter wavelengths, and the location of the brightest ring deviates from that of the birth ring (e.g., Figure 5 of Pawellek et al. 2019).

The classical expectation for the surface density power law index of collision-dominated debris disk tails is -1.5 (e.g., Strubbe & Chiang 2006;Thébault & Wu 2008), while the index is -2.5 when a disk is dominated by corpuscular and Poynting-Robertson drag from strong stellar winds (e.g., Strubbe & Chiang 2006). In Section 4.3 of Strubbe & Chiang ( 2006), the former is expected for 1 Ṁstar / Ṁ 10, while the latter for 10 2 Ṁstar / Ṁ 10 3 . The measured surface density power law index of -0.7 for M star TWA 7 here, however, deviates from the expectations, which makes the involvement of additional mechanisms necessary (e.g., companion radiation from HR 4796 B in HR 4796 A tail, or possible interstellar medium slowdown in HD 191089 tail).

In comparison with M star AU Mic, whose surface density power law index of -1.5 is explained in Strubbe & Chiang (2006) with Ṁstar / Ṁ 10, the existence of Ring 2 and Ring 3 makes it a more complicated scenario for the surroundings of M star TWA 7. Nevertheless, the detected CO gas that extends to ∼100 au, or ∼3 , around TWA 7 in Matrà et al. ( 2019) is likely a non-negligible source: gas can slow down the radially outward motion of the smallest dust particles; under this scenario, the slow-down explanation in Ren et al. (2019) might help interpret the power law index for the tail of Ring 1.

Ring 3 peak

The peak density radius, r max , for Ring 3 is likely located outwards in the STIS data (r max = 97 au) than in the SPHERE data (r max = 93 au). In addition, the α out power law index for Ring 3 is shallower in STIS (α out = -3.76) than in SPHERE (α out = -8.0). Both suggest that the smaller dust particles observed by STIS at ∼0.6 µm are likely located more extended than the larger dust particles observed by SPHERE at ∼1.6 µm.

The extension of small particles in STIS can be explained by the force balance between stellar radiation pressure and gravitational pull. For a particle with size a, the radiation force F radiation ∝ a 2 , while the gravitational pull F gravity ∝ a 3 . The ratio between the two forces is F radiation /F gravity ∝ a -1 . Therefore, in comparison with the relatively large particles seen in SPHERE, the smaller particles observed by STIS likely have relatively (rather than absolutely) higher ratio between the radiation force and the gravitational pull. This relatively higher ratio helps small particles reach higheccentricity orbits and thus arrive at further distances from the central star. Indeed, this is supported by the theoretical expectation that the optical depth is dominated by smaller particles at farther stellocentric separations, see, e.g., Equation (5) of Thébault & Wu (2008).

Possible moving components

Spirals

We confidently resolve S1 in SPHERE after removing the disk models, and in NICMOS after smoothing the residuals after disk model removal, see Figure 5. We tentatively detect S2 in STIS with an average pixelwise S/N of 2.313 after removing the disk model, and marginally observe it in SPHERE after disk removal, see Figure 5b. For S2, our STIS model might be influenced by quasi-static noise near the coronagraph. In addition, S2 is located in the northern region of the TWA 7 system, which is at the far side from Earth in our disk modeling: the marginal detection of S2 in SPHERE and NICMOS is thus likely caused by scattering effects and less observation time. Specifically, in comparison with forward scattering, less light is scattered in backward scattering. Therefore, we only quantify the S1 motion between the 1998 NICMOS and 2017 SPHERE observations.

To measure the motion of the S1 spiral, we follow the procedure in Ren et al. (2020). Specifically, in the surface density distribution map on the disk plane, we first fit Gaussian profiles to the radial profile at each angular position to obtain the peak position for the arm, then use dummy variables as proxies to obtain morphological parameters under polynomial description for the arm, and simultaneously quantify the arm motion rates under two hypotheses (i.e., either planetdriven motion or local Keplerian motion; Ren et al. 2020). In the planet-driven scenario, the entire spiral moves as a solid body, and the spiral pattern motion traces the orbital motion of the driver; in the local Keplerian motion scenario, the motion of the spiral is faster when its location is closer to the star. To minimize systematics in pixel size, we have convolved the residual images with a two dimensional Gaussian that has a standard deviation of 75.65 mas (i.e., 1 NICMOS pixel) for the images in Figure 5. We analyze the deprojected surface density maps on the disk plane to minimize stellar illumination and projection effects. We fit S1 in polar coordinates with p-degree polynomials as in Ren et al. (2020), and obtain minimized Akaike and Schwarz information criteria (AIC and SIC) at p = 1, which corresponds to linear description of an arm in polar coordinates in the solid body motion scenario.

We present the motion rates of S1 under the two mechanisms in Table 3. To account for possible instrumental and data reduction uncertainties (e.g., NICMOS centering, focal Note.

a The driver has a circular orbit along the midplane of the disk around a 0.46 +0.07 -0.10 M central star. b Enclosed mass within 20 au (i.e., star and disk mass combined), as inferred from local Keplerian motion. c r is the stellocentric separation in units of arcsec. plane mask location, NMF overfitting, surface brightness distribution difference in total intensity and polarized light: Ren et al. 2018b), we have adopted a relatively large 12 • uncertainty for NICMOS (i.e., ∼3 times the size of a NICMOS pixel at 1 ) when propogating the errors.

The measured motion direction of the S1 spiral is clockwise, which is consistent with the rotation of the CO gas in Matrà et al. (2019) where the east side of the disk moves away from Earth. To compare the two motion mechanisms for S1, we calculate the AIC and SIC for both mechanisms. Both criteria are modifications to the classical χ 2 statistic by adding penalty terms to dissuade excessive use of free parameters and avoid overfitting, and the model with the smallest AIC or SIC value is adopted as the best-fit model. Although the planet-driven model is slightly preferred, the difference in both information criteria is less than 1. 14 In comparison with the classical threshold of 10 for model selection (e.g., Kass & Raftery 1995), we cannot distinguish the two motion mechanisms here.

To investigate the nature of S1, we recommend follow-up SPHERE observations that probe the system with the same instrument and identical setup. When re-observed after 2021, a >4 year timeline will be established to help determine the motion mechanism for S1. Although it would be a shorter timeline than that presented in this paper, the identical setup and instrument can better constrain spiral motion: Xie et al. (2021) has recently demonstrated that PDI data with 1 year 14 Under independent Gaussian noise assumption. The difference is smaller when the noises are correlated.

Ring 2

We investigate the Ring 2 observations in the three instruments to explore its formation and motion mechanism. On the one hand, given that SPHERE best resolves Ring 2, we can inspect its surface brightness and thus optical depth distribution. On the other hand, given that we have visually detected Ring 2 in both SPHERE and STIS images, and that the radial profile for NICMOS reveals the possible existence of Ring 2, we can study its motion using the observations in all instruments that span from 1998 to 2019. Here we present two possible simple scenarios that can be tested using motion measurements with another SPHERE observation.

Static image: a resonance structure? -We remove the best-fit models for Ring 1 and Ring 3 from SPHERE, correct for not only the inverse-square law of stellar illumination but also the best-fit scattering phase function for Ring 2 to obtain the optical depth for the residuals (e.g., the procedure in Stark et al. 2014, except we here obtain the Henyey-Greenstein g parameter from disk modeling). We then deproject the disk to a face-on view, rotate the deprojected optical depth map to align the major axis with x-axis, and present the optical depth map that is normalized at the peak in Figure 6.

The relative under-density of materials in the northwest and southeast regions resembles resonance structures in-duced by exoplanets on debris disks (e.g., Ozernoy et al. 2000;Stark & Kuchner 2008;Goździewski & Migaszewski 2020). The overall geometry qualitatively resembles Figure 2b of Ozernoy et al. (2000), where the 2:1 and 3:2 resonances are in equal proportions for a system whose planetto-star mass ratio is 0.25 M Jupiter /M . Were TWA 7 system with a similar planet-to-star mass ratio, the planet should be located in the northwest underdensity region at r ≈ 50 au, and its mass is ∼2 M Neptune . We calculate in Section 4.3 the contrast limits for the NICMOS observations, and such a planet is beyound detection using the NICMOS data.

For a qualitative illustration of the possible resonance structures, we inspect the Exozodi Simulation Catalog15 and obtain the distribution of collisionless particles that are shepherded by a 5 M Earth planet located at 10 au from the Sun in Stark & Kuchner (2008). Noticing that we have removed in Figure 6a the Ring 1 model which controls the overall surface density for the TWA 7 system, we first smooth the simulation with a Gaussian kernel that has a standard deviation of 0.7 au, then remove the extended structures in the simulation by subtracting from it a Gaussian-smoothed version of the simulation (standard deviation: 1.4 au, i.e., high-pass filtering the data), see the resulting Poisson-noise-added image in Figure 6b. Under this scenario, although we have spatially separated Ring 2 and Ring 3 in our modeling, the two rings can be physically produced simultaneously by the same perturbing planet (see, e.g., Figure 6 The motion of the dip across the three observations can be explained with a shadow that moves at 19. • 8 ± 1. • 0 yr -1 clockwise. Under this scenario, a shadowing-casting disk under Keplerian motion at 5.3 +0.4 -0.5 au can explain the apparent motion for the underdensity region of Ring 2. Note: the 1998 data are prone to reduction artifacts, and the motion here could completely change with follow-up observations. Despite difficulties in directly imaging a planetary perturber for TWA 7, follow-up SPHERE observations will help investigate the existence of the underdensity regions, as well Ren et al. as their motion. Specifically, the motion of the underdensity regions in Figure 6a should trace the orbital motion of the planet. Nevertheless, we caution that our resonance inference here is solely based on a static image of Ring 2 and we have ignored Ring 1 in our qualitative comparison: a more dynamically motivated approach is needed to study all the rings and explore the resonance structure for the system.

Multi-epoch images: a shadowing structure? -We calculate the azimuthal profiles for Ring 2 to explore the apparent azimuthal motion of the ring (rather than onsite physical motion, see, e.g., Debes et al. 2017). Although we do not recover Ring 2 in HST disk modeling, we do observe plateaus in the corresponding radial profiles. For a stellocentric separation between 46 au and 59 au, in which Ring 2 resides, we calculate the median and standard deviation for the deprojected surface density maps with an angular step of 15 • . To measure the motion for the underdensity region, we first fit cosine profiles to the azimuthal profiles (e.g., Debes et al. 2017), then calculate linear angular motion for the dips, see Figure 7.

The 2017 SPHERE and 2019 STIS data resolve Ring 2 in Figure 1 visually. The reason that we cannot extract Ring 2 from STIS disk modeling is caused by the noise between Ring 1 and Ring 2 in Figure 1c. We calculate an apparent rotation rate of 15. • 5 ± 9. • 9 yr -1 clockwise between the two observations. Although we do not visually detect Ring 2 in the NICMOS data, the corresponding azimuthal profile can be fit with a cosine profile, and the dip in NICMOS is shifted ∼20 • from that in SPHERE, see Figure 7a.

The dip might have rotated for an additional 360 • between NICMOS and SPHERE that establish a 19 yr baseline. Such an additional rotation is likely. When we include this rotation, we can explain the motion of the dip in the three instruments with a 19. • 8 ± 1. • 0 yr -1 clockwise rotation in Fig- ure 7b. Under this scenario, for a 0.46 +0.07 -0.10 M central star, the rotation of the dip traces a shadowing source that is under Keplerian rotation at 5.3 +0.4 -0.5 au, or 0. 155 +0. • 012 -0. 015 that is just beyond SPHERE's coronagraph. Nevertheless, we reiterate that our NICMOS data reduction is prone to overfitting, and thus our measurement is not conclusive unless follow-up SPHERE observations are available: they would allow us to better investigate the motion mechanism of the Ring 2 dip (i.e., the northwest underdensity region of Figure 6a).

Clump

We recover a southern clump that is located at ∼3 from the star in all three instruments. Given that the proper motion of TWA 7 is (pm RA , pm Dec ) = (-118.75 ± 0.02, -19.65 ± 0.03) mas yr -1 in Gaia Collaboration et al. ( 2021), a background galaxy would have a relative motion of (∆RA, ∆Dec) = (2. 4772 ± 0. 0004, 0. 4099 ± 0. 0006) with respect to TWA 7 between NICMOS and STIS observations, the clump is therefore not a background galaxy. In comparison with the ALMA millimeter observations in 2016 in Bayo et al. (2019), the southern elongated 2σ structure in their Figure 4 (from ≈3. 5 to 5. 5, with a position angle of 195 • ) is consistent with the STIS clump in Figure 3c.

We have a clear detection of the ≈150 au 2 clump in STIS with a sensitivity of ∼1 µJy arcsec -2 . Although STIS best resolves the clump, it is expected to map dust particles that are smaller than those in NICMOS and SPHERE. We thus do not compare the STIS data with the other instruments for motion analysis; instead, we compare the NICMOS and SPHERE observations to estimate its motion.

The azimuthal motion of the clump between NICMOS and SPHERE is 3 • ± 1 • using cross-correlation, which corresponds to the Keplerian motion at 130 +50 -20 au from a 0.46 M central star. The derived location is consistent with the location of the clump that starts at ∼3 = 102 au within 2σ. In addition, there is a marginal evidence in Figure 5 that the motion rate of the clump decreases as a function of stellocentric separation. Nevertheless, since the clump extends beyond the edge of the NICMOS detector, and the NMF data reduction method may have altered its surface brightness distribution in NICMOS, we do not perform a clump motion measurement as for the S1 spiral. In addition, we cannot properly constrain the radial motion of the clump if it is moving outwards.

In comparison, Boccaletti et al. (2015) report apparently outward motion of ripple-like features for the edge-on debris disk orbiting M1 star AU Mic, and Chiang & Fung (2017) suggest that these features are sub-micron dust particles that are repelled by stellar wind. Limited by the NICMOS field of view, data quality, and data reduction artifacts, we cannot properly constrain the apparent motion for the clump in TWA 7. Nevertheless, the clump in the TWA 7 system might be under the same mechanism as the ripples in AU Mic. We fit Gaussian radial profiles to the surface density distribution of the clump in HST, and obtain peak location of 116 ± 4 au and 130 ± 5 au for NICMOS and STIS, respectively. Assuming STIS and NICMOS map the same dust source, this offset corresponds to a radial speed of 3.2 ± 1.4 km s -1 outwards. In comparison, the escape speed is 1.8 km s -1 at 130 au for a 0.5 M star. The clump is thus possibly unbound to the star; yet we reiterate that our calculation is based on the multiple assumptions that can pose challenges to our motion measurement.

Improving point source detection

To investigate the impact of disk modeling on point source detection, we calculate for the NICMOS observation the 5σ contrast curves under two scenarios: before and after removing the best-fit model from the observations. We inject the TinyTim PSF of TWA 7 to a physical on-sky position on Improvement Factor Disk Modeling Improvement Figure 8. 5σ detection limit of point sources before and after disk removal for NICMOS, and improvement factor defined by dividing the two curves. The the full width at half maximum areas for the Ring 1 and Ring 3 models are shaded. For a 6.4 Myr system, the mass of hot-start and cold start planets in Spiegel & Burrows (2012) are marked in red and blue, respectively. On the one hand, the disk removed contrast curve is deeper than the original one; on the other hand, the contrast improvement is the best near the location of Ring 1: both prove that disk modeling can help improve contrast.

the observations, perform reference differential imaging using KLIP (using the same number of references and components as NMF), then calculate the S/N for that position following the recipe in Section 3.4 of Debes et al. (2019). For each location, we vary the flux of the injected point source through a binary search in log scale, until when its S/N is within 0.01 from 5, the corresponding planet-to-star flux ratio is then the 5σ detection limit. We define contrast as the flux of the detection limit divided by 1.296 Jy (i.e., the NIC-MOS flux of TWA 7. See Section 5.1 for the calculation.).

For each on-sky radial location, we perform this calculation with an azimuthal step of 30 • , and calculate the mean contrast for all azimuthal positions. We repeat the above process for each radial separation between 5 NICMOS pixel and 53 NICMOS pixel (i.e., 0. 38 to 4. 01), and present the contrast curves in Figure 8. By dividing the two contrast curves, we notice that the disk-removed contrast curve is deeper than the original contrast curve by m = 16% ± 12%. What is more, the contrast improvement can reach m = 45% near the location of Ring 1, or m ≈ 25% near Ring 3. 16 Our comparison of the two curves demonstrates the necessity of disk modeling in detecting fainter planets.

We have assumed that there are no detectable planets in our contrast curve calculation for NICMOS observations. Were there hypothesized planets that should be detectable using the NICMOS data, such planets might have been overfit by the best-fit disk model. An ideal approach is to mask out different positions and perform disk modeling then contrast calculation, since such a procedure could extract these planets (e.g., PDS 70c: Wang et al. 2020) and consequently help detect fainter planets while enhancing the improvement in detection limit. We do not perform such a calculation, since we have not detected point sources with high confidence, and the presented results Figure 8 have demonstrated the necessity of disk modeling in detecting fainter planets.

Dust Properties

We study dust properties for the system using the best-fit models in the three instruments, see Figure 9a for the surface brightness radial profiles used for analysis in this Section. To compare the radial profiles, we have interpolated the best-fit images at a resolution of twice the NICMOS pixel size.

In comparison with Figure 1d, the recovered NICMOS surface brightness is at a level that is brighter than SPHERE, which confirms the expectation that NICMOS data reduction with NMF is under the influence of overfitting. We thus use the best-fit models to minimize such effects in comparing the system at different instruments. Nevertheless, we only recover the plateau at the Ring 2 location in the SPHERE model, we thus do not include the Ring 2 component in our analysis here.

Color

Using the best-fit NICMOS and STIS models, we calculate the color for the scatterers. To account for the differences such as stellar brightness at the two wavelengths and instrument response, we first calculate the unobstructed instrumental responses to the [START_REF] Kurucz | ATLAS9 Stellar Atmosphere Programs and 2 km/s grid[END_REF] star model of TWA 7 using pysynphot (STScI Development Team 2013). The pysynphot inputs are the effective temperature for TWA 7 (4018 K: Gaia Collaboration et al. 2018), its V -band magnitude (10.91: Torres et al. 2006), and its surface gravity log g in units of cm s -2 (4.18 : Stassun et al. 2018). For NICMOS F160W, the ObsBandpass input parameter is 'nicmos,2,f160w', and the corresponding instrument response is 0.953 Jy (yet we adopt a value of 0.953× 2.034×10 -6

1.496×10 -6 = 1.296 Jy to take into account the difference in the PHOTFNU parameters in the two NICMOS eras); for STIS, 'stis,ccd,a2d4', and 0.177 Jy. Then for each instrument, we divide the surface brightness distribution by the corresponding instrumental response, and obtain the relative brightness between the disk and the star. The integrated relative reflectance, F disk /F star , is 3.53 × 10 -3 in STIS, and 1.01 × 10 -3 in NICMOS.

We calculate the relative brightness radial profiles, then convert them to magnitude per square arcsec, and obtain the color radial profile by subtracting the converted NICMOS profile from the STIS one, see the results in Figure 9b. We observe that the relative reflectance in STIS is higher than that in NICMOS, making the scatterers blue. In addition, the system is bluer when stellocentric distance increases: Ring 1 is relatively brighter in STIS with ∆mag = -1.25, while Ring 3 has ∆mag = -1.75. Both phenomena are consistent with the expectations that M stars can retain small dust particles in their debris disks (e.g., Arnold et al. 2019), and that small particles which can scatter more light at shorter wavelengths are pushed further out than larger ones (e.g., Strubbe & Chiang 2006;Thébault & Wu 2008).

Polarization fraction

We investigate the polarization properties for the scatterers since NICMOS and SPHERE cover similar wavelengths: the former probes the system in total intensity at 1.4-1.8 µm, the latter in polarized light at 1.48-1.77 µm. We calculate the linear polarization fraction by dividing the surface brightness profile of the best-fit SPHERE model by that of the best-fit NICMOS model, see Figure 9b. Given that we cannot re-cover Ring 2 with NICMOS modeling, we only focus on Ring 1 and Ring 3 here.

The Ring 1 region has a peak polarization fraction of 85%, while the Ring 3 region peaks at 75%. High polarization is consistent with the explanation that the scatterers are small and/or not compact: on the one hand, high polarization is expected theoretically for small dust particles (0.1 µm -0.2 µm: Dabrowska et al. 2013;Perrin et al. 2015); on the other hand, similar levels of polarization fraction have been experimentally measured for fluffy aggregates [START_REF] Volten | [END_REF]. Given that there could exist non-negligible stellar wind activity around M stars, and that we cannot constrain the activity for TWA 7 using existing indicators such as Xray brightness (Section 1.1) and Ring 1 tail distribution (Section 4.1.2), we do caution that small dust particles created by collisional cascade can be efficiently removed by stellar winds.

Ring 1: Rayleigh scattering in total intensity?

In the Rayleigh scattering regime (i.e., the size of the scatterers is smaller than the observation wavelength by more than one order of magnitude), the scattering phase function in total intensity follows

I tot (θ) = 3 4 1 + cos 2 θ , (6)
while the polarization fraction follows Equation (4). We sample 5000 data points from the NICMOS phase functions from Ring 1 modeling for comparison, see Figure 10. The NIC-MOS phase functions do not have a large overlap with the Rayleigh phase function in total intensity. Nevertheless, we cannot conclude on the scattering regime for the TWA 7 disk seen in NICMOS.

Mathematically, the parametric Henyey-Greenstein phase function in Equation ( 3) is a monotonic function for 0 • ≤ θ ≤ 180 • , while the Rayleigh total intensity phase function is symmetric about θ = 90 • that cannot be depicted by the Henyey-Greenstein function. Therefore, we cannot test whether a system is under Rayleigh scattering in NIC-MOS total intensity using Henyey-Greenstein phase functions. Similarly, we cannot test whether the SPHERE data in polarized light is under the Rayleigh scattering regime with the modeling presented in this paper.

Physically, collisional cascade produces significantly more small dust particles than larger ones, i.e., n(a) ∝ a -3.5 where a is the particle size (e.g., Dohnanyi 1969). The size of the particles in debris disks can reach the observation wavelength in scattered light. The existence of particles with sizes comparable with observation wavelength makes it necessary to model them using Mie theory and its variations, or more complicated radiative transfer modeling methods.

To determine the scattering regime for Ring 1, especially to test the λ -4 dependence of scattered light intensity as a function of wavelength λ, future high resolution multiwavelength observations of TWA 7 in total intensity, as well as the actual Rayleigh scattering phase function, are needed.

Summary

From a favorable nearly face-on vantage point of the M star TWA 7 system, we have performed a multi-instrument and multi-decade characterization of its circumstellar disk in scattered light. By analyzing the TWA 7 observations with HST/STIS, HST/NICMOS, and VLT/SPHERE using disk modeling, we report a layered debris structure around M star TWA 7. We identify three rings (Ring 1, Ring 2, and Ring 3), one spiral (S1), and an ≈150 au 2 dusty clump. We also tentatively detect a secondary spiral (S2) that has low average pixelwise S/N. Our disk modeling using static geometric models shows that the extended Ring 1, which peaks at 0. 8, dominates the overall distribution of dust. The rest of disk components are superimposed onto the outskirts of the extended ring.

M stars have stellar winds bring corpuscular and Poynting-Robertson drag, and consequently their debris disk tail have a surface density power law index from -1.5 to -2.5 depending on the stellar mass loss rate Ṁstar (Strubbe & Chiang 2006). However, if Ring 1 is the birth ring of small particles (see Pawellek et al. 2019 for exceptions), the power law index of -0.7 for its tail is outside the expected range: it is even shallower than that for quiescent stars with Ṁstar 10 Ṁ . The existence of secondary CO gas in this system might help explain the power law index for Ring 1. With additional complicating factors such as the existence of Ring 2 and Ring 3, we cannot infer the stellar wind activity for TWA 7 with our power law measurements. Nevertheless, such constraints might be obtained by fitting the multi-instrument observations and retrieve dust properties such as size and porosity, or by a dynamically motivated study of the spatial distribution of the system (e.g., Schüppler et al. 2015).

We have compared the reference differential imaging detection limit of point sources for NICMOS before and after disk modeling. By modeling the NICMOS disk, we can detect planets that are fainter by ∼15% overall. What is more, at the Ring 1 location where the disk is the brightest, we can improve the detection limit by 45%. Our study demonstrates the necessity of disk modeling in detecting fainter planets for upcoming high-contrast imaging missions, e.g., Roman/CGI, LUVOIR, and HabEx.

We have compared the surface brightness distribution of the system in three instruments. On the one hand, using NICMOS and STIS total intensity images, we observe that the system has a blue color. The system is bluer when the stellocentric distance increases. On the other hand, using SPHERE and NICMOS images, we compare the system in polarized light and total intensity: ∼80% of the light is polarized. The former is consistent with the expectation that small particles scatter light more efficiently at shorter wavelengths, the latter is consistent with the expectation for small or fluffy dust. Nevertheless, we caution that the existence of stellar Ren et al. winds can change our inference. In addition, by comparing the total intensity phase functions between NICMOS and Rayleigh scattering, we find that the Henyey-Greenstein approach cannot be used to test the Rayleigh scattering regime.

We have calculated the rotation of the S1 spiral arm, and the radial motion of the clump, after removing the disk models from the observations. The S1 rotation can be described by both solid body movement and local Keplerian movement. The radial motion of the clump is possibly unbound to the system. Although the observations establish a twodecade timeline for motion analysis, multiple factors including the data reduction method, pixel size difference, stellar alignment, and different field of view of the detectors, keep us from determining the motion mechanisms for S1 and the clump.

We have inspected the morphology and apparent rotation for Ring 2, and investigated two possible scenarios for its formation that can be tested with new SPHERE observations. On the one hand, the optical depth distribution of Ring 2 in SPHERE might result from mean motion resonances with a planetary perturber that is located at the northwest underdensity region. Were the system under this scenario, the Ring 2 and Ring 3 components, which are spatially separated in our disk modeling, would be physically produced simultaneously by the hidden perturber. On the other hand, the azimuthal motion of an underdensity region (or a brightness dip) on Ring 2 can be explained by shadowing effects from a Keplerian rotating disk at 5.3 +0.4 -0.5 au. The layered debris disk around TWA 7 in scattered light reveals that M stars can host complex debris structures, and the motion of these structures can help understand their formation mechanisms. SPHERE re-observations after 2021 will establish a >4 year timeline to help constrain the motion for the S1 spiral and the clump, investigate the existence of the S2 spiral, and examine the morphology and motion for Ring 2. Being debris disks orbiting M stars, the almost face-on TWA 7 system offers a complementary view than the edge-on AU Mic system. Being members of the TW Hya association, the existence of the face-on debris disk around TWA 7, and the face-on protoplanetary disk around M star TW Hya (i.e., TWA 1), shows that circumstellar disks can evolve at different rates even for similar spectral types.

  arXiv:2105.09949v2 [astro-ph.EP] 26 May 2021 TWA 7 Layered Debris Disk Ren et al.

  Ren et al.

Figure 1 .

 1 Figure 1. Surface brightness distribution of the TWA 7 debris disk in log scale, the dotted transparent ellipses are the maximum density radii for the three SPHERE rings in Section 3.3. (a) 1998 HST/NICMOS F160W total intensity, the data are prone to over-fitting. (b) 2017 VLT/SPHERE H-band Q φ linearly polarized light. (c) 2019 HST/STIS total intensity. (d) Radial profiles for the three images, excluding the background star at ∼18 kpc idenified by Gaia in Section 3.1.2. The radial profile for NICMOS is subject to overfitting in data reduction, see Section 5 for that for the best-fit models.

  , S2 (a low S/N structure in STIS), clump, and background star regions in calculating the likelihood in disk modeling, see the annotations in Figures 2 and 3.

Figure 2 .

 2 Figure 2. Modeling of the SPHERE observation in Section 3.3.1. (a) is the SPHERE Q φ image and (e) is the best-fit model. The residual maps in (b), (c) and (d) have the individual ring models in (f), (g) and (h) cumulatively removed from (a). Notes: each ring has only one double-powerlaw component described by Equation (1); see Figure 4 for the surface density radial profiles.

Figure 3 .

 3 Figure 3. Modeling of the STIS and NICMOS observations in Sections 3.3.2 and 3.3.3, see (a) -(f) and (a') -(f'), respectively. For STIS, (a) is the observation and (d) is the best-fit model; the residual maps in (b) and (c) have the ring models in (e) and (f) cumulatively removed from (a). For NICMOS, we have adopted a forward modeling strategy to fit the three rings simultaneously, see see (d') for the best-fit model and (c') for the corresponding residuals. For the purpose of illustration only, (b') and (e') are the best-fit Ring 1 residual and model. Note: there is only one double-powerlaw component in each ring model; Ring 2 fit was attempted for both instruments, but was not confidently recovered.

Figure 5 .

 5 Figure 5. Spiral motion in Section 4.2.1: smoothed residuals after disk modeling in (a) 1998 NICMOS, and (b) 2017 SPHERE. For S1, the dashed lines are the best-fits from solid body motion; for the clump, approximate centers: the white ones are for the corresponding observation, and the transparent ones are for the other. Panel (c) shows the location measurement and motion fitting for S1, and the clockwise motion can be explained by two spiral motion mechanisms. Notes: the spiral locations are measured from surface density maps, the error bars here are 3σ, and the angular separation is calculated counterclockwise from the southeast semi-minor axis of the disk on the disk plane. separation using SPHERE is sufficient in constraining the motion mechanism for the spiral arms in the protoplanetary disk surrounding SAO 206462. Nevertheless, we note that other mechanisms can excite spirals in debris disks (e.g., Sefilian et al. 2021).

  of Stark & Kuchner 2008).

Figure 6 .Figure 7 .

 67 Figure 6. Mean-motion resonance scenario for Ring 2 in Section 4.2.2. (a) Deprojected optical depth of SPHERE residuals after removing the best-fit models for Ring 1 and Ring 3. The residuals resemble resonance structures induced by a planetary perturber (e.g., Figure 2b of Ozernoy et al. 2000). (b) High-pass-filtered and Poisson-noise-added simulation of mean motion resonance for collisionless dust particles whose radiation-to-gravity force ratio is 0.0023 in a system that has a 5 MEarth planet, marked by ⊗, located at 10 au from the Sun (Stark & Kuchner 2008).

Figure 9 .

 9 Figure 9. (a) Radial profiles measured on best-fit models. (b) Radial profiles for color and linear polarization fraction, measured from STIS-NICMOS after correction of instrumental response and SPHERE/NICMOS, respectively. Due to the non-detection of Ring 2 in NIC-MOS modeling, the polarization fraction for Ring 2 is expressed in dotted line. Note: the error bars in both panels are 3σ standard errors.

Table 1 .

 1 Observation log Note. Column 1: instrument. Column 2: target name. Column 3: filter. Column 4: central wavelength. λc for STIS is the pivot wavelength. Column 5: pixel scale. Column 6: HST aperture, or SPHERE coronagraph combination name. Column 7: inner working angle. IWA for STIS is the half-width of the wedge-shaped occulter. Column 8: total exposure time. Column 9: number of individual readouts. Column 10: total parallactic angle difference for HST, and field rotation for SPHERE. Column 11: observation UT date.

	Instrument	Target	Filter	λc	Pixel Scale	Aperture	IWA	Texp	N frame	∆θPA	UT Date
				(µm) (mas pixel -1 )		( )	(s)
	(1)	(2)	(3)	(4)	(5)	(6)	(7)	(8)	(9)	(10)	(11)
	STIS	TWA 7 CD-35 6480	Clear Clear	0.58 0.58	50.72 50.72	BAR5 WEDGEA1.0 BAR5 WEDGEA1.0	0.2 2638.80 0.5 4080.00 0.2 900.00 0.5 1340.00	3 × 6 2 × 20 • 3 × 2 2 × 20 • 1 × 6 • • • 1 × 2 • • •	2019 Feb 3 2019 Feb 3
	SPHERE	TWA 7	H	1.62	12.25	N ALC YJH S	0.1 6016.00 47 × 2	18. • 0	2017 Mar 20
	NICMOS	TWA 7	F160W 1.60	75.65	NIC2-CORON	0.3	607.88	2 × 3 1 × 29. • 9 1998 Mar 26

Table 2 .

 2 50th±34th percentiles of TWA 7 ring parameters

	Instrument Parameter	Ring 1	Ring 2 e	Ring 3 f
	θinc a	12.
	SPHERE

• 91 +0. • 18 -0.

• 19

Table 3 .

 3 Pattern motion measurement for S1

	Motion Pattern	Parameter	Value
		Rotation Rate (yr -1 )	1. • 3 ± 0. • 6
		Driver Location a (au)	33 +19 -10
	Solid Body	Orbital Period a (yr)	280 +240 -90
		AIC	5.68
		SIC	8.95
		Rotation Rate (yr -1 ) (1. • 0 ± 0. • 6) × (1 /r) 3/2c
	Local Keplerian	Enclosed Mass b (M ) AIC	0.3 +0.5 -0.3 5.98
		SIC	9.16

  Normalized total intensity phase functions between TWA 7 Ring 1 in NICMOS and Rayleigh scattering show little overlap. However, we cannot use the parametric Henyey-Greenstein phase function to depict Rayleigh scattering, see Section 5.3. Note: the 77 • -81 • region shaded with dots is contaminated by S1.

	Normalized Phase Function (at 90 )	0.6 0.8 1.0 1.2 1.4 1.6		NICMOS Rayleigh
		80	85 Scattering Angle ( ) 90 95	100
	Figure 10.

L X rescaled to match the Gaia Collaboration et al. (2021) distance.

The errors in this paper are 1σ unless otherwise specified.

An M1Vk star that is selected according to it being 3. • 6 away from TWA 7, and having a ∆(B -V ) = -0.14 and ∆V = -0.345 with TWA 7, see https://www.stsci.edu/hst/phase2-public/15218.pdf.

The effective wavelength in the observation is longer since TWA 7 is an M star.

https://archive.stsci.edu/prepds/alice/

We have ignored the ∼0.9% pixel size difference along X/Y directions inSchneider et al. (2003).

PHOTFNU = 2.03470 × 10 -6 Jy s DN -1 for NICMOS Era 1. See, e.g.,Hagan et al. (2018) for the two NICMOS observation eras.

Converted usinghttps://irsa.ipac.caltech.edu/data/SPITZER/docs/ dataanalysistools/tools/pet/magtojy/.

With transmission and integration time taken into account.

We have adjusted the code using Equation (1).

http://tinytim.stsci.edu

In comparison, Section 3.6 ofDebes et al. (2019) uses average S/N of 1 to calculate the detection limit of extended structures (specifically, rings) for STIS.

https://asd.gsfc.nasa.gov/Christopher.Stark/catalog.php

For a small number m (specifically, |m| < 0.6), we have 2.5 log 10 (1 + m) ≈ ln 10 • log 10 (1 + m) = ln(1 + m) ≈ m. The improvement in magnitude thus has a nearly identical value as the flux ratio in this study.

We thank the anonymous referee for their useful suggestions that make this paper more thorough. We thank Johan Olofsson and Rob van Holstein for discussion on SPHERE/IRDIS data reduction and calibration. Based on observations made with the NASA/ESA Hubble Space Telescope, obtained from the data archive at the Space Telescope Science Institute. STScI is operated by the Association of Universities for Research in Astronomy, Inc. under NASA contract NAS 5-26555. This research has made use of data reprocessed as part of the ALICE program, which was supported by NASA through grants HST-AR-12652 (PI: R. Soummer), HST-GO-11136 (PI: D. Golimowski), HST-GO-13855 (PI: E. Choquet), HST-GO-13331 (PI: L. Pueyo), and STScI Director's Discretionary Research funds, and was conducted at STScI which is operated by AURA under NASA contract NAS5-26555. The input images to ALICE processing are from the recalibrated NICMOS data products produced by the Legacy Archive project, "A Legacy Archive PSF Library And Circumstellar Environments (LAPLACE) Investigation," (HST-AR-11279, PI: G. Schneider). Based on observations collected at the European Organisation for Astronomical Research in the Southern Hemisphere under ESO programme 198.C-0209(F). Part of the computations presented here were conducted on the Caltech High Performance Cluster, partially supported by a grant from the Gordon and Betty Moore Foundation. Facilities: HST (NICMOS, STIS), VLT:Melipal (SPHERE). Software: DebrisDiskFM (Ren et al. 2019), emcee (Foreman-Mackey et al. 2013), IRDAP (van Holstein et al. 2020), nmf imaging (Ren 2020), pysynphot (STScI Development Team 2013), TinyTim (Krist et al. 2011).


Anaïs Begemann

Heinrich Sticht

Amber Begtrup

Antonio Vitobello

Laurence Faivre

Siddharth Banka

Bader Alhaddad

Reza Asadollahi

Jessica Becker

Tatjana Bierhals

New insights into the clinical

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