RZ Piscium Hosts a Compact and Highly Perturbed Debris Disk

RZ Piscium (RZ Psc) is well known in the variable star field because of its numerous irregular optical dips in the past 5 decades, but the nature of the system is heavily debated in the literature. We present multiyear infrared monitoring data from Spitzer and WISE to track the activities of the inner debris production, revealing stochastic infrared variability as short as weekly timescales that is consistent with destroying a 90 km sized asteroid every year. ALMA 1.3 mm data combined with spectral energy distribution modeling show that the disk is compact (∼0.1–13 au radially) and lacks cold gas. The disk is found to be highly inclined and has a significant vertical scale height. These observations confirm that RZ Psc hosts a close to edge-on, highly perturbed debris disk possibly due to migration of recently formed giant planets that might be triggered by the low-mass companion RZ Psc B if the planets formed well beyond the snowlines.


INTRODUCTION
Circumstellar disks of gas and dust are expected to form as a natural consequence of the star formation process, and they are sites that allow astronomers to study how planetary systems form and evolve.In the early gas-rich protoplanetary disk stage, these disks set the initial conditions for planet formation (Miotello et al. 2022;Manara et al. 2022); in the later gas-poor debris phase, these disks bear the imprints of planets and highlight the pathways of their formation and migration history (Wyatt 2008;Krivov 2010;Matthews et al. 2014;Hughes et al. 2018).A particularly interesting sub-class of debris disks is called "extreme debris disks" (EDDs) generally found around young stars that show an exceptionally large amount of warm dust in the terrestrial zone (see Figure 1 for an example).Because these EDDs are mostly found around stars younger than ∼200 Myr, corresponding to the end of terrestrial planet accretion (Chambers 2013;Quintana et al. 2016), the large amount of warm dust is commonly interpreted as the product of terrestrial planet formation (Melis et al. 2010;Meng et al. 2014).The terrestrial planet formation link for Gyr-old systems such as HD 69830 and BD+20 307 is less clear because the extreme dustiness could also be due to (1) intense dynamical stirring via the presence of planets for the first case (Lovis et al. 2006) or binarity for the latter case (Zuckerman et al. 2008), and (2) the consequence of unstable moons (Hansen 2023).

Su et al.
More than half of the EDDs are known to show photometry variability either in their total infrared outputs through Spitzer and WISE monitoring (Meng et al. 2014(Meng et al. , 2015;;Su et al. 2019Su et al. , 2020;;Rieke et al. 2021;Moór et al. 2021Moór et al. , 2022) ) or by irregular dips through routine all-sky optical surveys (Gaidos et al. 2019;Powell et al. 2021;Melis et al. 2021).The characteristics of the dips are best described as star-size dust clumps, presumably generated by large collisions within a few au, passing in front of the stars (de Wit et al. 2013).Because infrared emission directly traces the production and dissipation of freshly generated dust, the combination of infrared variability and optical dips provides a powerful means to extract the details of these large collisions such as the location and the minimum size and mass of bodies involved as illustrated in the HD 166191 system (Su et al. 2022).
The interpretation of the EDDs around young stars as a late stage of terrestrial planet formation is not unique.An alternative possibility is that the high dust levels are the result of transient dynamical clearing of planetesimal regions that might be analogous to the Late Heavy Bombardment in our solar system.Growing evidence has shown that the instability that shaped the reconfiguration of the giant planets in our solar system happened earlier than previously thought (Nesvorný et al. 2018;Clement et al. 2019;de Sousa et al. 2020) and was likely triggered by the dispersal of the gas disk (Liu et al. 2022).In addition to perturbation by newly formed giant planets, external excitements from nearly low-mass companions could also be the cause of instability.In such a scenario, the planetesimal disk might be truncated to a small radius with a significant scale height, a stage that corresponds remarkably well with the observed properties of the RZ Piscium (RZ Psc) system, the focus of this paper.
RZ Psc is a popular, solar-like target among the amateur community (Percy et al. 2010) because it is known to show sharp irregular, optical dips over the past five decades (Zaitseva 1978;Grinin et al. 2010;de Wit et al. 2013;Kennedy et al. 2017).The age of RZ Psc is heavily debated in the literature, ranging from a few Myr-old young system (Grinin et al. 2010) to an evolved giant star (Kaminski ȋ et al. 2000).The most recent comprehensive study by Punzi et al. (2018), including X-ray and multi-epoch optical, high-resolution spectroscopic observations, suggests an age of 30-50 Myr and reveals accretion activity from a reservoir of circumstellar gas.Compared to typical accreting protoplanetary disks, the gas accretion rate is low ( Ṁ ≲ 7 × 10 −12 M ⊙ yr −1 ) as estimated from an empirical magnetospheric accretion model (Potravnov et al. 2017).Optical spectro- scopic monitoring shows that the accretion rate is variable without known periodicity, and can reach one order of magnitude higher during outbursts (Dmitriev et al. 2023).The Galactic motion of RZ Psc calculated with new Gaia DR2 astrometric data suggests a possible membership in the Cas-Tau OB association with an age of 20 −5 +3 Myr (Potravnov et al. 2019).Gaia EDR3 gives a distance of 184.1±0.9 pc to RZ Psc (Gaia Collaboration et al. 2016, 2021), which is ∼5% closer than the value used in earlier studies.We adopt this distance for all the derived quantities in this paper.
Although RZ Psc has all the attributions of the socalled UXOR variability commonly found among premain sequence stars since the 1980s, the infrared excess indicative of a circumstellar dust disk was not identified until 2013 using all-sky WISE and AKARI data (de Wit et al. 2013).The broad-band excess Spectral Energy Distribution (SED) is well represented by blackbody emission of ∼500 K with an infrared fractional dust luminosity of ∼7×10 −2 ; follow-up N-band spectroscopy also shows a prominent 10 µm feature (Figure 1, Kennedy et al. 2017).There is no sign of cold dust (debris outside tens of au), although its presence can not be ruled out due to lack of far-infrared/millimeter observations.High contrast, ground-based imaging obtained with VLT/SPHERE recently shows that RZ Psc hosts a 0.12 M ⊙ companion at a projected separation of 22 au, and concludes that the disk (source of the infrared excess) orbits the primary star (Kennedy et al. 2020).The low-mass component is expected to have significant impact on the disk itself such as truncating the disk outer edge and heavily stirring the remaining planetesimals.The latter is consistent with the suggestion of a dynamically active asteroid belt revealed as numerous optical occultation events as the products of collisions (de Wit et al. 2013).
To further shed light into the nature of the RZ Psc disk, we present multiyear monitoring data from Spitzer and WISE to probe the short-(weekly) and long-term (monthly to yearly) infrared variability.Millimeter observations from SMA and ALMA are also presented to assess the cold disk properties.Section 2 describes the data and basic reduction used in this paper.Further analyses extracted from the data are presented in Section 3 including (1) the disk properties derived from the ALMA image and SED model, indicative of a compact disk with large scale height, and (2) the general implications of the intense activity in the inner zone of the planetary system from the Spitzer monitoring.In Section 4, we discuss the origin of the activity observed in the infrared, the potential sources of the gas reservoir, and the mechanisms that create the large disk scale height by linking with the intense dynamical activity revealed by photometric monitoring.A short conclusion is given in Section 5.

Warm Spitzer/IRAC
Warm Spitzer monitoring data were obtained under two GO programs PID 11093 and 13014 (PI Su) with observations covering from 2015 March to 2018 December.RZ Psc has two Spitzer visibility windows per year, each a total of ∼36 days in length.During the first year of the monitoring, cadences of ∼5±1 and ∼10±2 days were used to search for weekly variation due to material as close as 0.1 au from the star.We later used a finer cadence of 3±1 days in the 2016-2018 observations.A total of 95 sets of observations at both 3.6 and 4.5 µm bands was obtained.We used a frame time of 2.0 s with 15 cycling dithers (i.e., 15 frames per Astronomical Observation Request (AOR)) to minimize the intra-pixel sensitivity variations of the detector (Reach et al. 2005) at both bands, achieving a signal-to-noise ratio of ≳180 in single-frame photometry.These data were first processed with IRAC pipeline S19.2.0 by the Spitzer Science Center.We performed aperture photometry on each of the data sets, following the procedure outlined in Su et al. (2019).The final weighted-average photometry is given in Table A1.We also determined the instrumen- tal photometry repeatability by monitoring two isolated sources in the field of view.A <1% rms in the measured photometry was found for nonvarying sources in these data.
Using the stellar parameters listed in Table 1, we estimated the stellar photosphere in the two IRAC bands to be 34 and 22 mJy with a typical uncertainty of 1.5% limited by the accuracy of the optical and near-infrared photometry.After subtracting the stellar contribution by assuming the star is stable at these infrared wavelengths, the infrared excesses (i.e., the disk fluxes) and associated color temperatures are shown in Figure 2 over the four-year span in the units of Barycentric Modified Julian Date (BMJD).Overall, the infrared output of the system at 3.6 and 4.5 µm shows stochastic variations with a 50-60% maximum peak-to-peak flux difference.No significant periodicity was found using several pe- riodogram tools, except for spurious periods associated with the sampling windows.The median values of the excess fluxes are 24.0 mJy and 28.3 mJy at 3.6 µm and 4.5 µm, respectively (shown as horizontal dashed lines in the upper panel of Figure 2).Because the star dominates the noise in the measured photometry, the typical uncertainty in the two-band color temperatures is ∼40-50 K with an average of 1000 K (the horizontal dashed line in the middle panel of Figure 2).We note that the color temperature derived from the two shortest IRAC bands is not necessarily equal to the dust temperature in the system.Under optically thin conditions, the derived color temperature might reflect the true dust temperature, but the absolute numbers also depend on the assumed photospheric fluxes that are subtracted off from the measurements.A ∼10% change in the photospheric fluxes could result in changes to the derived color temperature on the order of ∼100 K in the case of RZ Psc.

WISE/NEOWISE and Optical Photometry
NEOWISE provides all-sky W1-and W2-band photometry every 6 months since 2014 (Mainzer et al. 2014), a time-domain resource for variability study.We extracted the single-exposure photometry from IRSA to fill the gaps between Spitzer visibility windows and to assess the long-term variability including data from the cryogenic WISE mission.For easy viewing, we binned the single-exposure data using a cadence of 3 days and plotted 1 them in Figure 2 when overlapping with the Spitzer data.Because the WISE/NEOWISE mission gives more than 10 years of W1 and W2 photometry, the multiyear WISE W1/W2 photometry is shown in Figure 3. WISE/NEOWISE data reveal that the overall infrared excess peaked in mid-2014 and gradually declined to the end of 2021, back to the 2010 level.
We also made use of the publicly available optical photometry from ASAS-SN Sky Patrol (https://asas-sn.osu.edu/,Shappee et al. 2014;Kochanek et al. 2017) and AAVSO.For ASAS-SN, we binned/combined and normalized the V − and g−band data following the procedure outlined by Su et al. (2022).In the AAVSO database, there were lots of data from different observers with slightly different filters.We selected data from five different observers who used V -band filters and have enough measurements to establish "zero" points just like treating different cameras used in the ASAS-SN survey.All the normalized optical photometry is shown in the bottom panel of Figure 2 (Qi 2003;Sault et al. 1995) were used to calibrate and generate the final continuum and line maps.Using natural weighting yielded a synthesized beam of 1. ′′ 3 × 1. ′′ 0 at a position angle (P.A.)=-81 • , and an root-mean-square (rms) noise level of 0.6 mJy beam −1 .RZ Psc was not detected in the continuum map.We also imaged the CO 2-1 line (rest frequency: 230.538GHz) with 1.2 km s −1 channel width, which achieved a 1. ′′ 2 × 0. ′′ 97 (P.A.=-80 • ) synthesized beam, and an rms noise level of 44 mJy beam −1 .Similarly, no line was detected.Because these non-detection limits are much less stringent than the ALMA observation (next Section), we will not discuss the SMA data further.

ALMA Observations
ALMA observations were taken on 2021 August 22 and 23 under the project 2019.1.01701.S (PI: K. Su).Two observing blocks each with 36 min of on-source integration were obtained using the Band 6 receiver and 49/44 12m antennas with baselines ranging from 47 m to 11.6 km under a mean precipitable water vapor (PWV) column of 0.3 mm.Three spectral windows were devoted to the continuum at 245 GHz with a 2 GHz bandwidth, while one was centered on the CO(2-1) 231 GHz line with 1920 channels (0.31 km s −1 ).J0238+1636 was the main calibrator for pointing, atmosphere, bandpass, phase and flux calibration.Calibrated data products were produced and provided by the North American ALMA Science Center using CASA version 6.1.1.Because RZ Psc is a young low-mass star, we also checked for potential mm variability due to chromospheric and/or accretion-related events.To search for point source emission, we plotted the baseline-averaged visibility amplitudes as a function of time for each date of observation.Since a point source has a constant amplitude as a function of baseline length, baseline averaging should reveal any point source variability as a function of time within the data set, yet all points were consistent with a constant flux density using the 2 and 30 s average data.Since the known mm flares have all taken place on timescales of minutes or less (MacGregor et al. 2020), there is no mm flaring event during our ALMA observation.
We imaged the continuum by combining all four spectral windows using natural weighting, resulting in an rms of 8 µJy beam −1 with a beam size of 0. ′′ 06×0.′′ 04 at a P.A. of 14.8 • .There appears to be one bright source detected near the expected RZ Psc position with a peak flux of 41.6 µJy beam −1 surrounded by a faint, elongated extended structure (Figure 4).The source extension (defined by the 3σ contour) on the synthesized image is 0. ′′ 20×0.′′ 07 at a P.A. of −13 • .Given the source's location (north of ∼28 • ), the synthesized beam is extremely elongated, but along a different P.A. from that of the extended structure.The comparison between the source extension and the synthesized beam, particularly the P.A. (15 • vs. −13 • ), suggests the source is spatially resolved.We will discuss the fidelity of the extended structure in Section 3.1.
Using the Gaia DR2 proper motion values (pmRA=27.544±0.028mas/yr and pmDEC=−12.521±0.019mas/yr) for RZ Psc, there is a small offset (−7 mas in both RA and Dec) between the detected and expected positions.The pointing accuracy is ∼10 mas with the synthesized beam and signal-to-noise above, i.e., this small offset is insignificant.We conclude that the emission detected in the 1.3 mm map is from the RZ Psc system, and most likely the dust emission around the star because the stellar photosphere is expected to be at 0.3 µJy.While there is a suggestive peak that coincides with the position of the companion RZ Psc B, it is only at the 1.5 σ level and is therefore entirely consistent with expected noise properties.At the distance of 184 pc and the H/K S magnitudes for the M-dwarf companion (Kennedy et al. 2020), the companion is at least 10 times fainter than RZ Psc at 1.3 mm.Although we cannot completely rule out the possibility of excess emission around the companion, it is more likely just a noise spike in the data.
For the CO(2-1) observation, we concentrated on the spectral windows that contain the CO 230.538GHz rest frequency in the two days of data by averaging in 1s time intervals to avoid time smearing.We used the cvel command to convert the velocities from topocentric (TOPO) to kinematic local standard of rest (LSRK).We then converted the measured heliocentric velocity (2 km s −1 , Punzi et al. 2018) to LSRK velocity of 0.91 km s −1 .After subtracting the continuum with the uvcontsub task, we used tclean to generate channel maps centered on the rest frequency, and found no significant emission on a channel-by-channel basis.The zeroth moment map of CO was then constructed by integrating channels within ±10 km s −1 of the star's LSRK velocity.Assuming the disk is unresolved, we determine the 3σ upper limit of CO (2-1) flux is 1.03×10 −2 Jy km s −1 .We will further discuss the upper limit for the CO gas in Section 3.1.

Disk Parameters from the ALMA Observation
To further evaluate the significance of an extended structure in the ALMA data and constrain the disk extent, we fit the observations in the visibility domain with simple parametric models.We adopt two models that are centered on the star: (1) the disk emission has a Gaussian-like profile, and (2) the disk emission comes from a narrow (unresolved) ring.The common parame-ters for both models are the P.A. of the disk major axis, the disk inclination, and total flux, and the different ones are in terms of disk size where the former one only constrains the Gaussian core size σ, which is related to the Full-Width-Half-Maximum (FWHM) of the disk emission (FWHM = 2.354 σ), and the latter estimates the peak radius (r 0 ) of the narrow ring.For the Gaussianlike model, we found the Gaussian σ=0.′′ 03 +0.′′ 03 −0.
′′ 02 , i.e., the best-fit FWHM of the disk emission is 0. ′′ 07 (just slightly larger than the beam) and at a P.A. of −13 • , consistent with the 3σ contour.The residual map and associated uv-data are also shown in Figure 4.For the narrow-ring model, we found r 0 =0.′′ 03 +0.′′ 04 −0.′′ 02 , similar to the Gaussiandisk model.The disk is really compact with the peak radius ∼0.′′ 03 (6 au at 184 pc), and only slightly larger than the beam with the maximum radius less than 0. ′′ 07 (13 au).Given the compactness of the disk, the other geometric parameters like the disk inclination and P.A. are less constrained.By combining the two models, we estimate the disk inclination angle is 83 • with a total flux density of 42±8 µJy at 1.3 mm.The large uncertainties associated with the derived parameters corroborate the fact that the disk is only marginally resolved (i.e., mm emission peaks at 6 au with a maximum radius less than 13 au) and the inclination and P.A. are not well constrained.Finally, the disk is unresolved vertically, and we can place an upper limit on the mm disk scale height <28 au (3 σ at 184 pc).
We can convert the 1.3 mm flux to a total dust mass (including grains up to mm-cm sizes) by assuming a dust opacity of 2.3 g cm −2 (Beckwith et al. 1990) and a dust temperature of 20-100 K 2 , which corresponds to 6.4-40×10 −3 M ⊕ .Using the 3-σ upper limit of CO (2-1) flux (Section 2.4) and an excitation temperature of 40 K, the corresponding upper limit of the cold gas mass is 1.0×10 −5 M ⊕ .The ALMA observation indicates a very low gas-to-dust mass ratio < 10 −3 in the system, confirming that RZ Psc has evolved out of the gas-rich protoplanetary stage and into the debris phase.

SED model
To test whether the observed broad-band SED from infrared to millimeter is largely consistent with disk parameters extracted from the ALMA observations mainly the disk outer radius), we employ SED modeling of a single disk component with the publicly available radiation transfer code, MCFOST (Pinte et al. 2022).For simplicity, we assume an azimuthally symmetric, wedgelike (i.e., a fixed opening angle) disk and fix the majority of the parameters in the model, and mainly explore four disk parameters: inner and outer disk radii, disk scale height and total dust mass.The disk surface density is assumed to follow a r −3/4 distribution where r is the radius.Because there is no/very little gas detected by ALMA, we set the disk flaring index to 1, i.e., the disk scale height, H, as H = H o (r/r o ) where r o = 10 au.The star with the parameters listed in Table 1 is the only heating source.We also fix the grain properties by adopting Astrosilicates as the composition in a powerlaw size distribution n(a) ∼ a −3.5 (where a is the grain radius) from a min = 0.1 µm (to match the prominent 10 µm silicate feature) to a max = 1 cm (to account for the millimeter emission).An initial search on a wide parameter space suggests that the inner radius of the disk is close to 0.1 au in order to account for the emission between 2-5 µm.The temperature at that radius is close to the dust sublimation temperature for sub-µm silicate-like grains.We, therefore, also fix r in = 0.1 au.
We employ the Python package emcee (Foreman-Mackey et al. 2013) to derive best-fit values for the disk outer radius, scale height and total dust mass.The results are shown in Table 1 and the best-fit SED is shown in Figure 5. Interestingly, the model millimeter flux is ∼3 times lower than the ALMA measurement with r out = 12 +1 −2 au.This is why the SED derived dust mass is much lower than the direct mass conversion from the millimeter flux (Table 1).A larger outer disk radius would increase the millimeter flux (i.e., the mass), but such models have too much emission in the 2-5 µm region.This suggests that the disk structure might not be as simple as we assume.For example, the disk might have multiple components (i.e., a separate cold component that emits mostly in the far-infrared and millimeter wavelengths) or the disk is not axi-symmetric (i.e., the radial temperature difference between the apses of an eccentric ring might account for the missing disk millimeter emission).A two-component model (inner and outer) would certainly fit the SED as long as the outer cold component does not contribute much of the infrared flux (i.e., either the cold component lacks sub-µm grains and/or the location of the cold component is outside of ∼5 au).We note that a cold component extending to 3 times larger than the ALMA size would also fit the observed SED, really illustrating the degeneracy of SED modeling.
We also note that the favored inclination angle is ∼42 • in the SED search.Although within the possible range measured by ALMA, such models produce too much optical polarization (∼5% level at V-band), >5 times higher than the observed averaged value (Shakhovskoi et al. 2003).In fact, we find that the polarization percentage can be significantly reduced if the disk is more inclined like i ∼75 • while allowing similar fits to other photometry points (i.e., inclination does not have a big impact on the SED models within the search range because the disk is not very optically thick).Because

Su et al.
i ∼75 • is more consistent with the values derived from the ALMA data and the optical transits, we adopt it as the best-fit value.Optical polarimetry is very sensitive to the disk inclination, therefore, polarization measurements provide much stronger constraints on the disk inclination angle.We will explore further constraints on the disk structure (inclination and the degree of asymmetry) using time-series optical polarimetric measurements in a coming paper (Anche et al. in prep.).
Finally, a scale height of 12 +2 −2 au at the radius of 10 au indicates that the disk has a very wide opening angle, ∼50 • .Such a large opening angle is in stark contrast with the nominal angle (a few degrees) measured in typical debris disks.Although the SED derived disk scale height is degenerate with many other parameters (i.e., the exact quantity is uncertain), a larger than normal debris disk scale height is evident in comparison of the SEDs shown in Figure 1.Both the 20 Myr-old systems, RZ Psc (an EDD) and HIP 61049 (a typical young debris disk), show solid-state features in the midinfrared that trace small grains in the optically thin part of the disk (usually on the surface of the disk above the midplane).The prominence of the mid-infrared features (i.e., the equivalent width of the feature) is known to be sensitive to the disk scale height in the SEDs of protoplanetary disks, which has a similar effect in the SED models of a dusty disk around white dwarfs (see Figure 4 in Ballering et al. 2022).We emphasize that the prominence of the mid-infrared features are also affected by other parameters in the model, particularly the amount of small grains.The disk scale height derived from our SED model should be taken with a grain of salt.Nonetheless, a linear size of 12 au at 184 pc would have an angular size of 0. ′′ 065, consistent with not being resolved by the ALMA observation.Based on the marginally resolved ALMA map and the SED behavior, we can conclude that RZ Psc hosts a highly inclined (∼75 • ), compact (∼0.1-13 au) but vertically extended debris disk.We will discuss the possible mechanisms to create this specific structure in Section 4.3.

Implication from Spitzer Observations
While RZ Psc has been well-known for its optical variability both in broad-band photometry and polarimetry over decades, its infrared variability was recently noticed by Kennedy et al. (2017) using initial WISE and AKARI data.It is clear that the infrared variability cannot come from the star itself because the star is not large enough to account for the infrared flux change, i.e., the changes must come from the circumstellar environment either due to (1) new debris generated in the collisions and aftermaths between large bodies such as ≳100 km size asteroids and planetary embryos or (2) heating variation in the existing circumstellar material such as the dust clumps seen partially blocking the star.Although both are likely, the second case would make it difficult to explain the increasing infrared flux, and the source of occultation (naturally produced by the first case) requires further explanation.For simplicity, we only consider that the infrared variability arises from the changes in the dust mass and distribution.
In a typical debris system where its infrared emission is sustained by the collisional cascades of km-size planetesimals, the production and loss of dust debris is balanced and no short-term variation is expected.Whether a source is variable or not is subject to the detectability, largely driven by the instrument sensitivity and the frequency of sampling.Because the change in flux is directly related to the change in the total cross section of the emitting material (∆Σ), one can quantify the detectability of variations for a given a sensitivity.For example, a flux change of 1 mJy at 4.5 µm (∆F 4.5 = 1 mJy) around a solar-type star at a distance of 200 pc requires a minimum change in the surface area of ∼12 stellar sizes (8×10 −4 au 2 ) assuming a dust temperature of 600 K. ∆Σ is scaleable as ∆F4.5  1 mJy [ 200 pc d ] 2 where d is the distance, and is sensitive to the dust temperature so that the hotter the temperature the smaller the surface area.In other words, within the current infrared sensitivity limits from WISE and Spitzer, the detectable infrared variability has to come from large-scale events, such as violent impacts between large bodies, to provide enough change in the surface area.
In addition to the infrared flux variability, the change in the color temperature derived from W1/W2 or IRAC1/IRAC2 bands can also shed light on the evolution of newly-generated debris.Under the optically thin condition, there are two mechanisms to manifest as a temperature change: one is a change in the dust location which goes as T d ∼ r −0.5 , and the other is a change in the dominant grain size which goes as T d ∼ a −1/6 .To extract the physical quantity, one needs to know the baseline level, which is difficult to establish for a constantly varying system like RZ Psc.
It is interesting to note that the Spitzer high-cadence (∼3-5 days) data reveal that the debris emission (i.e., the excess flux) varies on a weekly scale, and the longest period without a large degree change is about two weeks (∼4 sequential Spitzer measurements).Over the fouryear Spitzer monitoring period, the highest flux of the dust emission occurred in 2016 (near the display dates of 428 and 644 in Figure 2), ∼35.5 mJy at 4.5 µm, while the lowest flux occurred near the beginning of 2015 (the display date of 43) and the end of 2018 (the display date of 1380), ∼23.4 mJy at 4.5 µm.The variation behavior is consistent with the long-term behavior observed in WISE (Figure 3).Compared to the WISE/NEOWISE data, the highest flux observed by Spitzer is ∼86% of the peak flux in 2014, and the lowest is at a similar level ever observed (in 2010, 2016/2017).Using the average flux as the baseline, the typical change of flux at 4.5 µm is +7 and −5 mJy for increasing and decreasing cases, respectively.
Putting the observed infrared variability into context, one can estimate the minimum change in the emitting dust surface area (hence the minimum change in the associated dust mass) given an object's distance (184 pc for RZ Psc) and dust temperature (∼1000 K from Section 2.1).Adopting the typical flux changes observed at 4.5 µm (5-7 mJy), the change in the dust emitting surface area is (3.9-5.4)×10−4 au 2 , i.e., the associated dust mass is (0.3-1.1)×10 21 g assuming a dust density of 3.5 g cm −3 and a Dohnanyi-like size distribution (Dohnanyi 1969) with a minimum grain size of 0.1-0.5 µm.We note that these numbers are lower limits because of the assumption that the dust emission is optically thin, and are sensitive to the adopted dust temperature: ∼10% change in the assumed dust temperature would result in ∼50-100% change in the derived dust mass.A total dust mass of 1.1×10 21 g is roughly equivalent to the mass of a 90-km diameter asteroid.Such a flux change was recorded at least, if not more than, once per year.The RZ Psc system has experienced intense collisional activity, i.e., equivalent to total destruction of a 90-km asteroid every year (mass loss rate of ∼3.5×10 13 g s −1 ).It is unlikely such a high mass loss rate has persisted over the entire age of the system (20-50 Myr), i.e., requiring a mass equivalent to ∼4-9 M ⊕ .In other words, if such a mass loss rate persisted over a period of 10 3 (10 6 ) yr, the associated dust mass is roughly equivalent to the mass of the Moon (Earth).

Origin of High Collision Rate in the Inner Disk
ALMA 1.3 mm data combined with SED modeling show that the RZ Psc disk lacks cold gas and is composed of dust in a compact form (∼0.1-13 au radially), in stark contrast to typical UXOR-type young disks.As discussed extensively by Kennedy et al. (2017), the dust occulting events that dim the star in the optical must arise from a small fraction of the planetesimal population and lie at relatively high vertical positions above the main disk, either due to their high orbital inclinations or because they have been ejected into the inner region.The color change in the star during an occultation is similar to that expected from typical interstellar dust (Shakhovskoi et al. 2003), indicating that sub-µm particles play a dominant role in extincting the stellar emission.The presence of sub-µm grains is also consistent with the distinct 10 µm feature (Kennedy et al. 2017).The measurements with IRAC at 3.6 and 4.5 µm show a color temperature of ∼1000K, suggesting material exists inside an au.These observations indicate that much of the dust that dominates both the occultations and the 3-5 µm emission/variability is from grains small enough (≲ 1 µm for a solar-like system, Arnold et al. 2019) that they are ejected by radiation pressure force.This conclusion is consistent with the rapid fading of the events in the infrared discussed above and the failure of the occulting clumps to survive for a full orbital period of the inner disk (Kennedy et al. 2017).RZ Psc represents an extreme level of chaotic variability at 3-5 µm among all known debris disks, requiring episodic dust creation.Taking these observations together, they indicate a very high level of collisional activity among the planetesimals in the inner disk, that can produce large clouds of very finely divided dust and in extreme cases eject it well above the main planetesimal disk plane.
Can the short-term infrared variability result from gravitational stirring by the binary companion directly?The stability of planets and planetesimals in a binary system has been explored in many studies (Cuntz et al. 2007;Nesvold et al. 2016;Quarles et al. 2020), and strongly depends on the properties of the binary (the mass ratio and detailed orbital parameters of the companion), as well as the starting orbital parameters of the planetesimals/planets relative to the binary orbit.The current available properties of RZ Psc B (mass and projected distance) suggest that planetesimals in the outer ∼10-au region might be greatly influenced by the companion through the Kozai-Lidov Effect (Naoz 2016), but the ones within ≲ an au are well within the stable range.This result is intuitively reasonable, since the gravitational force from the companion star is less than 1/10,000 that from the primary star in the inner region.A number of items of circumstantial evidence also support this conclusion.First, binary stirring is likely to have persisted for a significant fraction of the life of the star, but the mass loss rates computed in Section 3.3 seem unreasonably high if the perturbation has been active over half of the lifetime of the system.Second, other binary systems do not have debris systems nearly as extreme as that of RZ Psc (Trilling et al. 2007;Yelverton et al. 2019).
An interesting alternative possibility is that the extreme debris activity is a result of planetary migration or planet-planet scattering.RZ Psc is of an age where planetary migration and scattering are thought to have Su et al. been rather common.Migration proceeds by scattering of planetesimals, which deflects the small-body population into highly eccentric and inclined orbits (Levison et al. 2007).During planet-planet scattering, "the planets roam around the system on eccentric and inclined orbits for an extended period" quoted from Marzari (2014), creating a period of chaotic evolution.Either mechanism shifts large numbers of planetesimals onto colliding orbits and will result in large amounts of collisionally produced dust.
The fate of this dust depends on its size (and composition).Grains ≲ 1 µm in radius will be ejected via radiation pressure force (Arnold et al. 2019), exhibiting rapid drops in the infrared flux (Su et al. 2019).Grains larger than the radiation blowout size would also drift outward radially with the aid of tenuous gas (Kenyon & Bromley 2016;Najita & Kenyon 2023).Larger grains will experience drag forces from the stellar wind (particularly important for young solar-like systems) and Poynting-Robertson drag.The conditions for this process around RZ Psc appear to be very similar to those for V488 Per discussed in Rieke et al. (2021).We adapt the estimate in that paper to find that where τ wind is the time for a grain size expressed as β (a ratio of radiation pressure force to gravitational force) to be dragged into the star from an initial distance of r 0 .For example, a grain with β = 0.01, corresponding to a radius up to ∼ 30 µm (Arnold et al. 2019), will be drawn from 0.50 to 0.49 au in of order a year.The change in orbit will expose it to an increased probability of a collision, likely resulting in a rapid rise in the infrared flux.In other words, the region in the inner au will be in turmoil, with objects of various sizes on different trajectories, i.e. rubble from recent collisions on eccentric and inclined orbits and moderate sized grains spiraling rapidly inward toward the star.The result will be a very high rate of collisions and production of debris dust, as observed.The behavior is consistent with expectations for the planetary scattering or migration hypothesis.

The Origin of Gas Reservoir
As mentioned earlier, the presence of circumstellar gas is inferred by the complex, optical line profiles observed in Hα, Ca II, Na I lines.The Hα profiles suggest that the emission structures have velocity widths of ≳300 km s −1 , originating within a few stellar radii of RZ Psc (Punzi et al. 2018).Using a magneto-accretion model designed to model gas accretion in T Tauri-like stars, Dmitriev et al. (2023) model the line profiles taken during the outburst in 2013 and find that the accretion rate is one order of magnitude higher during the outburst, and the magnetosphere extends from ∼5 to 10 stellar radii (∼0.03-0.06au) much larger than typical value for T Tauri-like stars.It is interesting to note that the derived inclination angle of RZ Psc is 43 • ±3 • with a relatively weak dipole field strength of ∼0.1 kGs.Given the disk inclination angle (∼75 • ), there appears to be a large obliquity in the RZ Psc system.If confirmed, the RZ Psc system would be the sixth solar-like system that exhibits large (≳30 • ) obliquity between the disk and the star as discussed by Hurt & MacGregor (2023).Dmitriev et al. (2023) speculate that the extended and weak magnetic field of RZ Psc could result in a long-lived primordial accretion disk.Assuming the circumstellar gas is primordial (i.e., H 2 dominated and an ISM-like CO-to-H 2 mass ratio of 10 −4 ), the mass of the gas reservoir (≲10 −7 M ⊙ , estimated from the ALMA CO upper limit in Section 3.1) could sustain the low accretion rate over 10 5 years.Given the age of the system and lack of cold gas as traced by CO observation, there is no clear evidence that the hot gas is primordial as suggested by the magneto-accretion model.
Tenuous gas has been found in some young debris disks (Hughes et al. 2018).Punzi et al. (2018) suggested that the hot gas revealed from the optical spectroscopy could originate from the accretion of material from a hot Jupiter-like planet swallowed by the star.The aftermath of such an engulfment might have been observed recently in ZTF SLRN-2020 (De et al. 2023).If true, the detailed evolution of the hot gas in RZ Psc might provide an additional test bed for such a phenomenon.The gas could also originate from out-gassing comets scattered toward the star, collisions of volatile-rich asteroids, or stripped atmospheres from giant impacts between protoplanets.Najita & Kenyon (2023) recently review these processes and conclude that all of them are capable of generating 10 −3 -10 −2 M ⊕ of secondary gas for stars younger than a few 100 Myr.If the gas is secondary in origin, the ALMA CO gas limit (< 10 −5 M ⊕ ) strongly disfavors the giant impact origin as the estimated mass of stripped atmosphere is on the order of 10 −4 M ⊕ for each giant impact (Najita & Kenyon 2023).It is difficult to differentiate the remaining two mechanisms because both asteroids and comets could store volatile material when they formed, and would later release it in the inner region.However, inward-scattering of comets is a low efficiency process (Bonsor et al. 2012) and often requires a chain of low-mass planets and/or a massive reservoir of icy planetesimals.Even accounting for the entire mm flux, there is no massive reservoir of icy planetesimals (material outside ∼10 au) as the potential source of scattered comets.Vigorous collisions between volatile-rich asteroids, consistent with the numerous observed optical dips, is more likely the source of the secondary gas in RZ Psc.4.3.The Mechanisms for the Puffed-up Disk Scale Height For gas-poor debris disks, the disk scale height is one of the observables that can constrain a system's dynamical excitation, and hence provide information about the processes and bodies shaping it under the influence of gravity (Wyatt & Dent 2002;Quillen et al. 2007;Thébault 2009).Similar to gas-rich protoplanetary disks, the disk vertical structure is also strongly stratified in terms of grain sizes -the larger the grains, the closer to the midplane in gas-poor debris disks.The disk scale height inferred from small grains and/or gas species has less diagnostic power in constraining the size of the stirring bodies due to the higher influence of nongravitational forces such as the radiation pressure and gas drag (Thébault 2009;Olofsson et al. 2022).A disk scale height derived from SED modeling would have different behavior from that described above (not an observed quantity).As stated in Section 3.2, it is uncertain to pin point the exact number for the disk scale height due to the marginally resolved ALMA image and SED model degeneracy, but a larger than nominal disk scale height appears to be present in the RZ Psc system.The question is what causes it?
The newly discovered, low-mass companion, RZ Psc B, is expected to play an important role in the behavior of the debris system, as discussed by Kennedy et al. (2020).The Hill radius of the companion is ∼10 au assuming a mass of 0.12 M ⊙ at a separation of 22 au.The disk outer radius (≲13 au derived from the ALMA observation, and 12 +1 −2 au derived from the SED model) is consistent with truncation by the companion.Furthermore, as discussed by Nesvold et al. (2016), a stellarmass perturber orbiting exterior to and inclined to the disk midplane can excite planetesimals into high eccentricity via the Kozai-Lidov mechanism (Naoz 2016), creating debris in significantly extended vertical structures compared to those caused by an interior planet.Based on single epoch data, little is known about the orbital parameters of the companion (except for an estimated mass from its brightness, and projected distance), but it locates at a P.A. of ∼75 • that is very different from the disk's P.A. (−16 •+44 • −20 • ) measured by ALMA.It is probable that the relative inclination between the external companion and the disk falls into the allowable range where the eccentric Kozai-Lidov mechanism is active (Naoz 2016).Because the Kozai timescale is inversely proportional to the orbital period of planetesimals, a hallmark for the Kozai perturbation will be highly eccentric particles found in the outermost region closer to the perturber, while particles closer to the star have relatively lower eccentricities.This is in line with the discussion in Section 4.1 that the collision activity in the inner au region is less affected by RZ Psc B.
Nonetheless, the presence of the external companion might provide the initial trigger for the planet migration hypothesis responsible for the inner activity.In this case, the birth place of the migrating planets is likely to be well beyond the snowline and close to the outer region of the disk where the influence of the companion is most effective.Furthermore, the potential misalignment between the companion and the disk midplane would create a potential warp or spiral structure, manifest as a puffed-up disk scale height when viewing close to edge on.In this case, one would expect asymmetric structure between the two disk ansae.That is, if the intense activity in the inner zone is caused by fully formed giant planets that are in low eccentricity orbits, the debris disk in RZ Psc is expected to be axi-symmetric like many of the Kuiper-belt analogs revealed by ALMA (Matrà et al. 2023) and vertically thin (Nesvold et al. 2016).If the low-mass companion played a significant role in shaping the disk dynamical activity, we should expect an asymmetric disk (Stuber et al. 2023).Future higher resolution disk images will provide further evidence to differentiate the two scenarios.
Can planet migration alone create a large disk scale height?The small bodies in our own Kuiper belt have a bimodal distribution (dynamically cold with low inclinations and dynamically hot with high inclinations, Brown 2001).The bimodal distribution in the small bodies is thought to be a telltale sign of Neptune's migration (Nesvorný 2015).Among the exoplanetary systems, β Pic is the only one known to show a clear bimodal distribution using both data in optical scattered light (tracing small grains) and millimeter emission (tracing large grains) (Golimowski et al. 2006;Matrà et al. 2019).The vertical structure along with the disk asymmetry observed in scattered light and the cold CO gas (Dent et al. 2014) in the β Pic system are likely caused by the migration of a yet-to-be discovered super Neptune (Matrà et al. 2019).
The large vertical scale height in the RZ Psc system is consistent with the suggestion in the preceding section (Section 4.1) that the intense collisional activity in the innermost part of the system arises through the effects of planetary migration.The asteroid belt in our solar system has experienced significant depletion (99%) com-pared to its primordial phase (Weidenschilling 1977) and the early depletion was linked to the giant planet migration as suggested by the Grand Tack model (Walsh et al. 2011).If there were giant planets in the RZ Psc system, they must have fully formed already because of the lack of cold gas in the system.In our solar system, fully formed giant planets played a critical role in shaping the terrestrial planets (Raymond et al. 2014).In this case, RZ Psc might represent the beginning stage in cleaning its primordial asteroid belt under the influence of newly formed giant planets.In this case, the observed disk scale height is directly linked to the required mass of the scattering planet.If the RZ Psc disk does have an opening angle as large as 50 • (unlikely as we discuss earlier due to the limitation of SED modeling), the planet needs to be, at least if not more than, as massive as Jupiter (with an escape velocity of 60 km s −1 ) to form a scattered disk at ∼5-10 au.Such a planet (Jupiter-like at <5 au) should have detectable signal in radial velocity data.Sporadic radial velocity variability was detected by Punzi et al. (2018), but the origin is unknown and a future detailed, long-term radial velocity study would provide better constraints on the presence of giant planets in the RZ Psc system.

CONCLUSION
We present multiyear infrared monitoring data from Spitzer and WISE to track the activities of inner debris production/destruction in RZ Psc.Millimeter observations along with the SED modeling provide a good assessment of the overall disk properties.These observations confirm that RZ Psc hosts a close to edge-on, highly perturbed disk, consistent with numerous optical dips caused by dust clumps in the past five decades.
ALMA 1.3 mm data obtained in 2021 reveal a compact, marginally resolved disk with a synthesized beam of 0. ′′ 06×0.′′ 04.Parametric models in fitting the visibilities suggest the millimeter emission has a radial extent less than 13 au, inclined by 83 • .The total 1.3 mm flux is 42±8 µJy, corresponding to a dust mass of 6.4-40 ×10 −3 M ⊕ .No cold CO(2-1) line was detected, suggesting <1.0×10 −5 M ⊕ for the CO gas mass.The ALMA data indicate a very low gas-to-dust mass ratio <10 −3 , confirming that RZ Psc hosts a gas-poor debris disk and has evolved out of the gas-rich protoplanetary disk stage.Adopting an axe-symmetric disk geometry, SED models suggest the disk inner radius is near ∼0.1 au, close to the dust sublimation radius, in order to produce enough infrared excess in the 3-5 µm region.The disk outer radius is at 12 +1 −2 au (consistent with the derived ALMA value) with a large scale height, suggesting a wide open-ing angle of ∼50 • .Although the disk scale height and inclination are degenerate in SED models, the detected percentage of the optical polarization (∼1%) further suggests that the disk is inclined by ∼75 • , consistent with the ALMA derived value and that a large disk scale height is necessary to fit the prominent 10 µm feature.The SED derived millimeter flux is 3 times lower than the ALMA value, suggesting the disk might have a more complex structure than axi-symmetric.An additional cold component that only emits at far-infrared or millimeter wavelengths, or an eccentric disk where grains extend further out radially both can account for the missing flux in the SED.
The four-year, high-cadence (∼3-5 days) 3-5 µm data reveal that the RZ Psc system exhibits stochastic weekly infrared variability with a typical flux change of ±5-7 mJy at 4.5 µm over a year if not more.Adopting a dust temperature of ∼1000 K under the optically thin condition, the 4.5 µm flux variation corresponds to a minimum change in the dust surface area of ∼5×10 −4 au 2 or 10 21 g.The intense collisions in the inner planetary zone totally destroy a mass equivalent to a 90-km-size asteroid on yearly basis.If such a high rate has persisted over a period of 10 3 (10 6 ) years, a Moon-(Earth-)size object would have been destroyed.
Although there is no cold CO gas detected in the RZ Psc system, the complex, variable optical line profiles do suggest the presence of close-in circumstellar gas with a very low accretion rate (10 −12 M ⊙ yr −1 ) if assuming an H2-dominated, primordial disk.Alternatively, the hot gas could be the remnant of a Jupiter-like planet engulfment as suggested by Punzi et al. (2018).Furthermore, the gas could also be secondarily created by out-gassing comets, collisions of volatile-rich asteroids, or stripped atmospheres from giant impacts between protoplanets (Najita & Kenyon 2023).Collisions between volatilerich asteroids are likely the source of the secondary gas given that (1) the low efficiency in inward scattering of comets and that RZ Psc does not have a massive cold reservoir of Kuiper-belt zone outside ∼10 au, where Jupiter-family comets in our solar system originate, and (2) stripped planetary atmospheres are expected to release more gas than observed.
It is evident that RZ Psc hosts a debris disk highly perturbed either by recently formed giant planets and/or by its low-mass companion RZ Psc B. The intense collision activity in the inner au region is likely a result of planetary migration or planet-planet scattering, both deflecting large numbers of planetesimals onto colliding orbits and resulting in large amounts of collisionally produced dust manifest as stochastic infrared variability and numerous optical dips.RZ Psc might be in the begin-ning stage of giant planet migration that has been initiated by gravitational effects caused by the companion if the migrating planets were formed beyond the snowline.Furthermore, the large disk scale height might also be related to the low-mass companion.In this case, one would expect the disk to show warp or spiral structures induced by the external perturber similar to the ones found around HD 106906 (Fehr et al. 2022).Follow-up observations to better determine the orbital parameters of the companion, and future high-resolution disk images will shed more light into the nature of the highly perturbed RZ Psc disk.

APPENDIX
A. WARM SPITZER PHOTOMETRY FOR RZ PSC Table A1 shows the 3.6 and 4.5 µm photometry of the RZ Psc system obtained during the Spitzer warm mission as described in Section 2.1.The excess emission and its uncertainty were derived by the subtraction of the expected photospheric value (34 and 22 mJy at the 3.6 and 4.5 µm bands, respectively) with a typical uncertainty of 1.5% of that value added in quadrature.This research has made use of the NASA/IPAC Infrared Science Archive, which is funded by the National Aeronautics and Space Administration and operated by the California Institute of Technology.This work has also made use of data from the European Space Agency (ESA) missionGaia (https://www.cosmos.esa.int/gaia),processed by the Gaia Data Processing and Analysis Consortium (DPAC, https://www.cosmos.esa.int/web/gaia/dpac/consortium).Funding for the DPAC has been provided by national institutions, in particular the institutions participating in the Gaia Multilateral Agreement.
The Submillimeter Array is a joint project between the Smithsonian Astrophysical Observatory and the Academia Sinica Institute of Astronomy and Astrophysics, and is funded by the Smithsonian Institution and the Academia Sinica.

Figure 1 .
Figure1.SED comparison for three different types of debris disks around solar-like stars: RZ Psc (red, extreme debris disk), HIP 61049 (green, typical young debris disk) and HD 69830 (blue, highly stirred warm belt) with the latter two normalized to the distance of RZ Psc.Photometry points are shown as various symbols, mid-infrared spectra as the solid-line, and model disk emission as the dashed line.The amount of infrared excess emission in EDDs is exceptionally large compared to typical debris disks.

Figure 2 .
Figure2.Spitzer monitoring data of the RZ Psc system: top for the 3.6 and 4.5 µm excess fluxes, middle for its corresponding color temperature, and bottom for the optical light curves.The arrow-like symbols near the top of the first panel mark the year boundaries of the data.Horizontal dashed lines in the upper two panels are the median value of the derived quantity.Optical light curves include ASAS-SN V -(orange) and g-band (green) data, and AAVSO V -band (various colors) data (for details see Section 2.2).

Figure 4 .
Figure 4. (a) A naturally weighted ALMA image of the 1.3 mm continuum observation centered at RZ Psc (star symbol).Contours are [1, 3, 5]×σ, where σ= the rms noise = 8 µJy beam −1 .The synthesized beam dimensions, represented by the red ellipse in the lower left corner, are 0. ′′ 06×0.′′ 04 with a position angle of 14.8 • .The position of RZ Psc B (Kennedy et al. 2020) is marked by the red plus for reference.The middle and right panels are the corresponding residual maps after the subtraction of a Gaussian-disk model in (b) and the associated uv-data in (c) (Details see Section 3.1).

Figure 5 .
Figure5.SED models of the RZ Psc system including the star and a single-component disk (see Section 3.2 for details).Various symbols are publicly available broad-band photometry points from Simbad, 2MASS(Cutri et al. 2003), ALLWISE(Cutri & et al. 2014), and AKARI(Ishihara et al. 2010), and the yellow dashed line is the N -band grism spectrum fromKennedy et al. (2017).

Facilities:
Spitzer (IRAC), WISE, ALMA, IRSA, AAVSO We thank the referee for providing constructive comments that improved the clarity of this manuscript.This work has been supported by NASA ADAP programs (grant No. NNX17AF03G and 80NSSC20K1002).H.B.L. is supported by the National Science and Technology Council (NSTC) of Taiwan (Grant Nos.111-2112-M-110-022-MY3).The paper is based on observations made with the Spitzer Space Telescope, which was operated by the Jet Propulsion Laboratory, California Institute of Technology.This paper has made use of data from AAVSO, and we acknowledge with thanks the variable star observations from the AAVSO International Database contributed by observers worldwide and used in this research.The MCMC part of the work was supported by the University of Arizona High-Performance Computing (HPC) resources.This paper makes use of the following ALMA data: ADS/JAO.ALMA#2019.1.01701.S. ALMA is a partnership of ESO (representing its member states), NSF (USA) and NINS (Japan), together with NRC (Canada), MOST and ASIAA (Taiwan), and KASI (Republic of Korea), in cooperation with the Republic of Chile.The Joint ALMA Observatory is operated by ESO, AUI/NRAO and NAOJ.

Table 1 .
Parameters and references for RZ Psc