# JWST-TST High Contrast: Asymmetries, dust populations and hints of a
collision in the $\beta$ Pictoris disk with NIRCam and MIRI.

Isabel Rebollido Centro de Astrobiología (CAB CSIC-INTA) ESAC Campus Camino
Bajo del Castillo, s/n, Villanueva de la Cañada, 28692, Madrid, Spain  Space
Telescope Science Institute, 3700 San Martin Drive, Baltimore, MD 21218, USA
Christopher C. Stark  NASA Goddard Space Flight Center, Greenbelt, MD 20771,
USA  Jens Kammerer European Southern Observatory, Karl-Schwarzschild-Straße 2,
85748 Garching, Germany  Space Telescope Science Institute, 3700 San Martin
Drive, Baltimore, MD 21218, USA  Marshall D. Perrin Space Telescope Science
Institute, 3700 San Martin Drive, Baltimore, MD 21218, USA  Kellen Lawson
NASA Goddard Space Flight Center, Greenbelt, MD 20771, USA  Laurent Pueyo
Space Telescope Science Institute, 3700 San Martin Drive, Baltimore, MD 21218,
USA  Christine Chen Space Telescope Science Institute, 3700 San Martin Drive,
Baltimore, MD 21218, USA  Dean Hines  Space Telescope Science Institute, 3700
San Martin Drive, Baltimore, MD 21218, USA  Julien H. Girard Space Telescope
Science Institute, 3700 San Martin Drive, Baltimore, MD 21218, USA  Kadin
Worthen William H. Miller III Department of Physics and Astronomy, John’s
Hopkins University, 3400 N. Charles Street, Baltimore, MD 21218, USA  Carl
Ingerbretsen William H. Miller III Department of Physics and Astronomy, John’s
Hopkins University, 3400 N. Charles Street, Baltimore, MD 21218, USA  Sarah
Betti Space Telescope Science Institute, 3700 San Martin Drive, Baltimore, MD
21218, USA  Mark Clampin Astrophysics Division, Science Mission Directorate,
NASA Headquarters, 300 E Street SW, Washington, DC 20546, USA  David
Golimowski Space Telescope Science Institute, 3700 San Martin Drive,
Baltimore, MD 21218, USA  Kielan Hoch Space Telescope Science Institute, 3700
San Martin Drive, Baltimore, MD 21218, USA  Nikole K. Lewis Department of
Astronomy and Carl Sagan Institute, Cornell University, 122 Sciences Drive,
Ithaca, NY 14853, USA  Cicero X. Lu Gemini Observatory/NSF’s NOIRLab, 670 N.
A’ohoku Place, Hilo, HI 96720, USA  Roeland P. van der Marel Space Telescope
Science Institute, 3700 San Martin Drive, Baltimore, MD 21218, USA  Emily
Rickman European Space Agency (ESA), ESA Office, Space Telescope Science
Institute, 3700 San Martin Drive, Baltimore, MD 21218, USA  Sara Seager
Department of Earth, Atmospheric and Planetary Sciences, Massachusetts
Institute of Technology, Cambridge, MA 02139, USA  Department of Physics and
Kavli Institute for Astrophysics and Space Research, Massachusetts Institute
of Technology, Cambridge, MA 02139, USA  Department of Aeronautics and
Astronautics, MIT, 77 Massachusetts Avenue, Cambridge, MA 02139, USA  Rémi
Soummer Space Telescope Science Institute, 3700 San Martin Drive, Baltimore,
MD 21218, USA  Jeff A. Valenti Space Telescope Science Institute, 3700 San
Martin Drive, Baltimore, MD 21218, USA  Kimberly Ward-Duong Department of
Astronomy, Smith College, Northampton MA 01063, USA  C. Matt Mountain
Association of Universities for Research in Astronomy, 1331 Pennsylvania
Avenue NW Suite 1475, Washington, DC 20004, USA

###### Abstract

We present the first JWST MIRI and NIRCam observations of the prominent debris
disk around Beta Pictoris. Coronagraphic observations in 8 filters spanning
from 1.8 to 23 $\mu$m provide an unprecedentedly clear view of the disk at
these wavelengths. The objectives of the observing program were to investigate
the dust composition and distribution, and to investigate the presence of
planets in the system. In this paper, we focus on the disk components,
providing surface brightness measurements for all images and a detailed
investigation of the asymmetries observed. A companion paper by Kammerer et
al. will focus on the planets in this system using the same data. We report
for the first time the presence of an extended secondary disk in thermal
emission, with a curved extension bent away from the plane of the disk. This
feature, which we refer to as the “cat’s tail”, seems to be connected with the
previously reported CO clump, mid-infrared asymmetry detected in the southwest
side, and the warp observed in scattered light. We present a model of this
secondary disk sporadically producing dust that broadly reproduces the
morphology, flux, and color of the cat’s tail, as well as other features
observed in the disk, and suggests the secondary disk is composed largely of
porous, organic refractory dust grains.

††facilities: JWST††software: python, astropy, spaceKLIP

##  1 Introduction

Beta Pictoris ($\beta$ Pic) is perhaps the most iconic and well studied
planetary system. Since the landmark discovery of $\beta$ Pic’s bright, edge-
on debris disk (Smith & Terrile, 1984)—making it the first exoplanetary system
ever imaged—$\beta$ Pic has been studied across the electromagnetic spectrum
using both ground- and space-based observatories. Two decades ago, high
angular resolution HST STIS coronagraphic observations revealed an inclined
inner disk, the outer edge of which was used not only to argue for the
presence of a putative perturbing companion, but also to constrain its mass
and orbital distance (Heap et al., 2000). Since then, ground-based, near-
infrared, high contrast imaging using VLT/NaCo has revealed the presence of
$\beta$ Pic b, a 10 MJup companion with a semi-major axis of $\sim 10$ au
(Lagrange et al., 2009), suggested to be consistent with the observed inner
warp of the disk (Dawson et al., 2011). More recently, long-term HARPS radial
velocity monitoring revealed the presence of a second planet, $\beta$ Pic c,
interior to the first with a mass 9 MJup and a semi-major axis of $\sim 2.7$
au (Lagrange et al., 2019), and subsequently confirmed by interferometry with
VLTI/GRAVITY (Nowak et al., 2020).

Multi-wavelength imaging is a powerful strategy for studying disks, as their
appearance changes with wavelength, potentially revealing different components
(e.g., Hughes et al., 2018). Millimeter continuum images are primarily
sensitive to thermal emission from large grains whose orbits are largely
unaffected by radiation pressure and drag effects. As a result, ALMA images
are believed to reveal the planetesimal birth rings. Visual observations, on
the other hand, are primarily sensitive to scattered light from the smallest,
sub-micron-sized grains. These grains are sensitive to radiation pressure and
are expected to be blown out more efficiently around higher mass stars. As a
result, HST STIS scattered light observations have revealed very extended
disks, extending to hundreds of astronomical units from their parent
planetesimal belts (Schneider et al., 2018). Mid-infrared observations at
$10-20~{}\mu$m are sensitive to thermal emission from warmer particles, and
can reveal disk components that elude observation at either visible or
millimeter wavelengths, especially with the great leap in sensitivity and
resolution enabled by JWST. For instance, mid-infrared observations of
Fomalhaut with JWST spatially resolved the warm inner disk and revealed an
unexpected new intermediate dust belt, providing additional clues to dynamical
processes within that complex system (Gáspár et al., 2023).

Detailed modeling of the history of our Solar System indicates that, over its
lifetime, the dust production rate was dominated by individual collisions.
However, whether individual collisions can affect the appearance of a debris
disk is not well understood. Early modeling of debris disks assumed that the
systems were in steady-state, constantly producing dust via collisions or
outgassing and removing dust via radiation pressure and corpuscular stellar
wind drag (e.g. Strubbe & Chiang, 2006). However, more recent Spitzer IRAC
time monitoring of extreme debris disks (with optical depths
$\tau=L_{\mathrm{IR}}/L_{*}>0.01$) discovered that large collisional events in
massive disks may be common (Meng et al., 2015). Spitzer IRAC observations
have also revealed fast and dramatic increases in the near to mid-infrared
excess that have been attributed to collisional destruction of asteroid sized
objects in the terrestrial planet zone (Su et al., 2022). Since these
collisions occur very close to the star, the dust they produce is very
challenging to spatially resolve. In parallel, high SNR HST scattered light
imaging has discovered density enhancements in the more distant Kuiper Belt
regions of debris disks (Stark et al., 2014; Jones et al., 2023).

The dynamical outcome of a collision or outgassing event largely depends on
the grain sizes produced. Whether a grain is radiatively driven out of the
system or gravitationally bound is determined by $\beta$, the ratio of the
force due to radiation pressure compared with that due to gravity (Artymowicz,
1988). Small grains tend to be more affected by radiation pressure and have
larger values of $\beta$, while large grains can be relatively unaffected and
have $\beta\sim 0$. However, grains with the same $\beta$ may have a variety
of different grain sizes, porosities, and compositions (Arnold et al., 2019).
The near- and mid-infrared imaging enabled by JWST is expected to be broadly
sensitive to grains near the blow out size limit ($\beta\sim 0.5$) for typical
debris disk host stars.

Several lines of evidence suggest that a significant dust production event may
have recently occurred in the $\beta$ Pic system. A clump of excess mid-IR
emission detected on the southwest side of the disk (Telesco et al., 2005; Li
et al., 2012) has been interpreted as potentially due to production of dust
from a massive collision of very large planetesimals ($\gtrsim 100$ km).
Moreover, high contrast observations from the ground presented in Skaf et al.
(2023) report a complex structure inwards of $\sim$ 60 au, with different
clumps of material. Recent publications argued whether the main clump’s
position changes over time as if in Keplerian motion around the star (Skaf et
al., 2023), consistent with the recent collision hypothesis, or is instead
stationary (Han et al., 2023). Collisions and outgassing events should release
gas as well, and indeed ALMA observations resolved a compact clump of gas
spatially co-located with the mid-IR excess clump, detected first in 12CO
(Dent et al., 2014) and subsequently in C I (Cataldi et al., 2018). Modeling
of gas production and photodissociation rates for the clump, as well as the
lack of azimuthal dispersion of the gas and dust, imply the production event
was recent (within the last few thousand years or less; Cataldi et al., 2018).

Even a well-studied system such as $\beta$ Pic is still capable of surprises.
The Mid-Infrared Instrument (MIRI) coronagraphy presented in this work has
revealed previously unseen additional components of the system: a sharply
curved feature which extends up dramatically away from the disk plane, an
extended asymmetric disk on the opposite side, and complex surrounding
nebulosity. We label the newly-seen curved feature of the $\beta$ Pic disk as
the “cat’s tail” based on its shape, and present below a dynamical hypothesis
for its origin from radiation pressure acting on small particles in the
aftermath of collisions.

The remainder of this paper is organized as follows. Section 2 presents the
observations and data reduction processing, including high contrast Point
Spread Function (PSF) subtractions and PSF deconvolutions. Section 3 presents
results on disk morphology, describing the observed asymmetries and
substructures seen in the disk, including the “cat’s tail”, followed by
Section 4 which presents results on the spatially resolved disk spectral
energy distribution (SED) and colors. In Section 5 we present a dynamical
hypothesis to explain the origin of the “cat’s tail” extended feature. We
discuss the system comparing results presented here to previous observations
in Sect. 6; and we summarize and concude in section 7.

###  1.1 Programmatic context

Our collaboration is conducting several complementary studies to characterize
the $\beta$ Pic system with JWST. A companion paper by Kammerer et al.
(submitted) studies the planet $\beta$ Pic b using the same NIRCam
coronagraphy data presented herein. $\beta$ Pic b is not visible in the MIRI
coronagraphy data, because the thermal emission from the planet is hidden by
the much greater thermal emission from the disk at 15 and 23 µm. Separate MIRI
MRS observations of $\beta$ Pic (program 1294; PI C. Chen) have yielded
thermal emission spectra of the disk, and also the planet at shorter mid-IR
wavelengths where the disk is not yet so bright. These observations will be
presented in Worthen et al. (submitted) and Chen et al. in prep. NIRSpec IFU
observations for additional near-IR spectroscopic characterization of the disk
were just obtained at the end of 2023 (program 1563; PI C. Chen). The
characterization of the planetary system will be continued by NIRISS Aperture
Masking Interferometry (AMI, Sivaramakrishnan et al., 2023) observations
(program 2297; PI T. Stolker), which aim to study the presence of clouds in
the atmosphere of the inner planet “c” in particular. Collectively these
programs are bringing a large fraction of JWST’s full observing capabilities
toward the detailed characterization of this key exoplanetary system.

This paper is part of a series to be presented by the JWST Telescope Scientist
Team (JWST-TST)111https://www.stsci.edu/~marel/jwsttelsciteam.html, which uses
Guaranteed Time Observer (GTO) time awarded by NASA in 2003 (PI Matt Mountain)
for studies in three different subject areas: (a) Transiting Exoplanet
Spectroscopy (lead: N. Lewis); (b) Exoplanet and Debris Disk High-Contrast
Imaging (lead: M. Perrin); and (c) Local Group Proper Motion Science (lead: R.
van der Marel). A common theme of these investigations is the desire to pursue
and demonstrate science for the astronomical community at the limits of what
is made possible by the exquisite optics and stability of JWST. The first
results in each of these areas are presented in Grant et al. (2023), (Ruffio
et al., 2023) and Libralato et al. (2023), respectively.

The program on $\beta$ Pic presented here was designed in particular to
exercise JWST’s coronagraphy modes over their full available wavelength range
from NIRCam F182M to MIRI F2300C, on a key target that has played a central
role in the study of exoplanetary systems since the field’s inception.

##  2 Observations and Data Reduction

We describe in turn the observations themselves, the reduction and
calibrations of images individually and combined via PSF subtractions, and
deconvolutions to reduce blurring of details by the JWST coronagraphs’ off-
axis PSF. $\beta$ Pic’s disk is sufficiently bright to be immediately visible
at very high SNR in all the observations, obvious in raw uncalibrated files
before reduction or PSF subtractions. The disk is even spatially resolved in
the MIRI target acquisition images, taken through a neutral density filter and
with an exposure time of 2 seconds. Nonetheless we perform carefully optimized
data reductions and image processing to maximize the visibility of fine
details in the disk.

###  2.1 Observations

$\beta$ Pic was observed with MIRI (Rieke et al., 2015; Wright et al., 2023)
and the Near Infrared Camera (NIRCam; Rieke et al., 2005, 2023) on board JWST
(Gardner et al., 2023). Observations were performed as part of GTO program
1411 (PI C. Stark) on 2022 Dec 13 (MIRI) and 2023 March 18 (NIRCam)222The
NIRCam observations were initially attempted on December 22, 2022, but failed
due to a guiding issue. The observations were subsequently rescheduled roughly
3 months later, at a position angle rotated by approximately 90 degrees from
the original plan.. The observations used the coronagraphic imaging modes of
both instruments (Boccaletti et al., 2022; Beichman et al., 2010; Girard et
al., 2022). A summary of instrument settings is shown in Table 1.

Our NIRCam observations targeted $\beta$ Pic in six different filters from
$\sim 1.7$–$5~{}\text{\textmu m}$, with filter coverage chosen to be sensitive
to the presence of potential disk constituents such as water and CO ices. We
used NIRCam’s 210R coronagraphic mask for short wavelength (SW) observations;
these were the first-ever coronagraphic science observations taken with SW
channel. We used the 15 and 23 µm coronagraphic filters for the MIRI
observations, chosen for maximum sensitivity to thermal emission from the
disk.

We adopted the common high contrast strategy of observing two roll angles on
the science target plus a PSF reference calibrator. We observed the bright
star $\alpha$ Pic as the PSF reference and included small grid dithers (SGD,
Lajoie et al. (2016)) for the MIRI 4QPM/F1550C and all NIRCam observations.
SGD were deemed not necessary for the Lyot/F2300C combination, as the target
acquisition accuracy was good enough, and the PSF is fairly extended. For
MIRI, we obtained images of a nearby empty patch of sky to subtract
astrophysical backgrounds as well as the “glow-stick” thermal background
artifact discovered post-launch (Rigby et al., 2022). To maximize the spatial
extent of the disk visible, we positioned the disk’s major axis at $\sim
45\deg$ to the 4QPM quadrant boundaries and the MIRI Lyot mask support strut
using a position angle constraint. Similarly, we placed the disk diagonally
across the NIRCam coronagraphic field of view. For all filters we observed the
disk at two different roll angles separated by $\sim 10^{\circ}$. Observations
were grouped in one non-interruptible sequence per instrument, to minimize the
effects of temporal variations.

The total on-source exposure time was roughly one hour per each of the 8
filters used. We chose integration times and readout patterns to avoid
saturation based on ETC calculations, which is particularly an issue at
shorter wavelengths where the star is brighter. Integration times and total
exposure duration on the brighter PSF calibrator were roughly half that on the
science target.

Instrument | Obs. dates | Coron. Mask/Filter | $\lambda$ | $\Delta\lambda$ | Readout Pattern | Exp. Time | Aperture PAs
---|---|---|---|---|---|---|---
|  |  | $(\micron)$ | $(\micron)$ |  | (s) | (degrees)
NIRCam | 2023 Mar. 18 | MASK210R/F182M | 1.84 | 0.24 | RAPID, 4, 90 | 1885.47 $\times 2$ | 84.3$\arcdeg$, 94.3$\arcdeg$
MASK210R/F210M | 2.09 | 0.21 | RAPID, 4, 90 | 1885.47 $\times 2$
MASK335R/F250M | 2.50 | 0.18 | BRIGHT2, 10, 80 | 1797.62 $\times 2$
MASK335R/F300M | 3.00 | 0.32 | BRIGHT2, 10, 80 | 1797.62 $\times 2$
MASK335R/F335M | 3.36 | 0.35 | SHALLOW4, 10, 35 | 1871.53 $\times 2$
MASK335R/F444W | 4.42 | 1.02 | SHALLOW4, 10, 35 | 1871.53 $\times 2$
MIRI | 2023 Dec. 13 | 4QPM/F1550C | 15.50 | 0.90 | FASTR1, 100, 55 | 1331.18 $\times 2$ | 353.8$\arcdeg$, 3.9$\arcdeg$
Lyot/F2300C | 22.75 | 5.50 | FASTR1, 100, 55 | 1799.49 $\times 2$
Table 1: Observing log for JWST observations of the $\beta$ Pic disk in
program 1411. The $\lambda$ column gives the pivot or center wavelength for
each filter, and $\Delta\lambda$ gives the filter bandwidth. The readout
column gives detector readout pattern, number of groups, and number of
integrations. The exposure times listed are per roll, times two rolls as
indicated for the total combined observing time on target. PSF calibrator and
background observations were obtained with similar settings as the above.

###  2.2 Data reduction through stage 2

We downloaded the raw uncalibrated data (“uncal” files) from the MAST archive
and processed them with spaceKLIP333https://github.com/kammerje/spaceKLIP
(Kammerer et al., 2022; Carter et al., 2023). SpaceKLIP invokes the official
jwst stage 1 & 2 pipelines444https://github.com/spacetelescope/jwst (used in
its 1.12.1 version Bushouse et al. (2023), with calibration files CDRS_CTX
$jwst\\_1090.pmap$ and $jwst\\_1088.pmap$ for MIRI and NIRCam respectively)
with parameters optimized for coronagraphy (discussed in Carter et al.
(2023)), and provides further image processing and analysis steps to perform
background subtraction, bad pixel cleaning, removal of 1/f noise, and PSF
subtraction using either principle component analysis (PCA)-based or classical
methods as described below. Achieving high-quality PSF subtractions depends
critically on careful attention to pixel-level calibrations and bad pixel
handling, and on image alignment, so substantial effort has been spent in
optimizing these steps and their parameters.

For the MIRI data, as in Carter et al. (2023), we skipped the pipeline’s dark
subtraction step because the available darks for the coronagraph subarrays are
of poor quality (low signal-to-noise ratio and a large number of bad pixels
from e.g., cosmic rays). We set the pipeline’s cosmic ray jump rejection
threshold to 8, since this value was found to yield the best contrast
performance by Carter et al. (2023), and adjusted the saturation correction
step parameters so that only the left/right/upper/lower neighbors of a
saturated pixel get flagged as bad pixels, as opposed to an entire 3 by 3
pixel box around a saturated pixel. In addition to these changes already
implemented by Carter et al. (2023), we also skipped the pipeline’s flat-
fielding step because the available flats do not yet correctly describe the
throughput loss due to the four-quadrant phase mask edges and the Lyot mask
support, which created bright bar-like artifacts in the reduced images. By
skipping the flat-fielding step, a flat-fielding error of typically less than
1% for the Lyot mask and less than 2% for the four-quadrant phase masks (as
long as further than $\sim 15$ pixels from the subarray edge) is introduced
and propagated through our data reduction process. We discarded the first
frame of each MIRI integration due to an increased level of noise caused by
reset switch charge decay (Carter et al., 2023). Then we identify and clean
bad pixels as described in more detail in the next subsection. After bad pixel
cleaning, we subtracted the corresponding background images from the science
target and reference star images to remove the glow stick artifacts
(Boccaletti et al., 2022; Carter et al., 2023).

For the NIRCam data, we handled dark subtraction, jump rejection, and
saturation correction similarly as for the MIRI data, but with a jump
rejection threshold of 4 instead of 8. Moreover, we defined a custom border of
four reference pixels around the edges of the NIRCam subarrays to remove
offsets caused by detector bias drifts as discussed in Carter et al. (2023).
SpaceKLIP now includes a custom pipeline step for correlated 1/f stripe noise
removal on coronagraphic data. It works on ramp data after jump detection by
first making a preliminary estimate of the slope, subtracting that from each
group to yield a residual image, fitting a Savitsky-Golay filter to the
residual image to estimate the 1/f noise contribution in that group, and
subtracting the noise model group-by-group prior to the ramp fit. Given the
high SNR of this dataset, subtracting the 1/f noise yields a mostly cosmetic
improvement in noise levels far from the star, but has negligible effect on
the bright regions. Next, we median-subtracted each frame to remove any
remaining background and cleaned bad pixels in a similar way as for the MIRI
data.

###  2.3 Bad pixel cleaning

We paid particular attention to bad pixel cleaning. We combined the DO_NOT_USE
pixels from the pipeline’s data quality extension with a custom bad pixel map
of hot pixels. For MIRI, these additional hot pixels were identified by visual
inspection of the MIRI background images. Through this procedure, we flag 28
additional hot pixels for the F1550C subarray and 48 ones for the F2300C
subarray. For NIRCam, we first run the pipeline until and including the KLIP
PSF subtraction step with the default DO_NOT_USE bad pixel map. Then, we
inspect the PSF-subtracted frames by eye to identify additional hot pixels and
re-run the entire pipeline, this time with an additional custom bad pixel map.
For the short wavelength subarray, we flag no additional bad pixels, and for
the long wavelength subarray, we flag 28 ones on top of the default DO_NOT_USE
bad pixels. In both cases, all bad pixels are cleaned using a succession of
three different algorithms. The first algorithm replaces pixels that are only
bad in some (but not all) integrations with their median value in the good
integrations. The second algorithm replaces remaining bad pixels with the
median of their four direct neighbors (only considering the non-bad
neighbors). The third algorithm replaces all remaining bad pixels with an
image plane median filter of size five pixels using the
SciPy.ndimage.median_filter function.

###  2.4 Image Alignments

PSF subtraction depends sensitively on proper image registration of science
and PSF reference observations. We recentered and aligned the NIRCam data as
follows. For the first science frame of each dataset, we determined the true
position of the host star $\beta$ Pic behind the coronagraphic mask using
phase cross-correlation with a simulated model PSF from
WebbPSF_ext555https://github.com/JarronL/webbpsf_ext (Leisenring, 2021; Girard
et al., 2022). For this procedure, the bright core of the PSF was masked out
as it is typically not represented well in simulations. Then the same Fourier
shift and subtract least squares minimization was used to align all subsequent
science and reference frames to the first science frame of each dataset. This
alignment procedure has sufficient precision to recover the injected few-
milliarcsecond small-grid dithers for the reference star and the $\sim 1$ mas
jitter ball from the pointing stability of _JWST_ (Rigby et al., 2023). As
discussed further in the companion paper by Kammerer et al. (submitted), for
this program the alignment analyses show nominal performance for the short
wavelength observations (F182M & F210M), but the long wavelength observations
(F250M, F300M, F335M & F444W) have an unusually large $\sim 50$ mas pointing
offset between the science and reference stars relative to the coronagraph
mask. This slightly impacts the PSF subtractions at small separations,
relevant to the planet $\beta$ Pic b, but has negligible impact on the disk
seen at wider separations.

We used the measured offsets to shift interpolate the PSF reference data to
the science data prior to subtraction. To avoid ringing artifacts (i.e., Gibbs
phenomena) from interpolating undersampled images, the F250M and F300M data
were first convolved with a Gaussian filter. The filters had a full-width half
maximum of 1.96 pixels and 1.71 pixels, respectively, chosen to equal 125% of
the Nyquist sampling criterion of $\lambda/2.3D$, where $D=5.2$ m is the
effective pupil diameter in this mode behind the NIRCam round mask Lyot
stop666https://jwst-docs.stsci.edu/jwst-near-infrared-camera/nircam-
instrumentation/nircam-coronagraphic-occulting-masks-and-lyot-stops.

For MIRI the situation was much simpler. Following Carter et al. (2023), we
simply assumed that the CRPIX position recorded in the FITS WCS header is the
coronagraph mask center and also the star location for both science and
reference targets, and aligned the images based on that. This simple
assumption neglects the uncertainties in MIRI coronagraph target acquisition
and in the precise knowledge of the exact mask center. These uncertainties are
of order $\sim 10$ mas, at least an order of magnitude less than the PSF FWHMs
at these wavelengths. In any case the high brightness of the disk relative to
the star at these wavelengths means that this simple approach suffices.

###  2.5 PSF Subtraction strategies

We tested and compared several methods for PSF subtractions, including
classical and KLIP reference differential imaging (RDI), angular differential
imaging (ADI), and model-constrained reference differential imaging (MCRDI).
We ultimately found that for MIRI, classical RDI yields excellent results,
while for NIRCam, the best disk images were achieved via MCRDI.

Figure 1: The $\beta$ Pic debris disk as seen with NIRCam and MIRI. The top 6
panels for NIRCam result from MCRDI PSF subtraction as described in section
2.5.2, while the bottom two panels for MIRI are from classical PSF
subtractions as described in section 2.5.1. In the NIRCam data, looking
closely the planet $\beta$ Pic b can be seen just to the left (NE) of the
occulted star (indicated with an arrow in the F300M data), along with the
edge-on debris disk and its subtle warp already known from prior observations.
In the MIRI data, a highly asymmetric curved feature we dub the “cat’s tail”
is seen to the west (SW), apparently bending sharply away from the disk plane
(see Fig. 4 and Sect. 3.2). A mask has been applied to the F2300C image only
to block the pixels behind the opaque classical Lyot occulter and its
supports. All images are shown rotated 60.9∘ to align the main disk major axis
horizontally, on identical asinh stretches and identical spatial scales,
though with a larger vertical field of view shown for the more extended
nebulosity seen with MIRI.

####  2.5.1 Classical PSF Subtraction

For MIRI, the high brightness of the disk, relative faintness of the star, and
extremely stable PSFs at these long wavelengths allow excellent results from
straightforward classical PSF subtraction (Krist et al., 1997). After
alignment, we scaled the reference PSF images to match the observed stellar
flux and then subtracted them. We scaled the reference images by minimizing
the residuals in the PSF-subtracted images in regions where only speckle noise
from the stellar PSF, but no significant disk flux was present. We then
rotated the subtracted images from each roll to place north up. Finally, we
averaged the images while weighting individual pixels by the coronagraphic
mask throughput of the 4QPM for F1550C and the opaque Lyot occulter and its
supports for F2300C, yielding the final images that are shown in the lower
panels of Figure 1. Uncertainties after PSF subtraction were simply computed
from those generated through the JWST pipeline, and first order propagations
were applied accordingly.

We performed similar classical subtractions for the NIRCam data, but
ultimately the MCRDI approach described in the next section was found to yield
superior images.

####  2.5.2 Model-Constrained Reference Differential Imaging

Compared to classical RDI, techniques like KLIP (Soummer et al., 2012) and
LOCI (Lafrenière et al., 2007) may better reproduce the stellar PSF’s
morphology by combining a set of reference images to minimize residuals with
the science data. When circumstellar flux is present, such techniques
necessarily overestimate the brightness of the starlight and therefore remove
astrophysical signal as a side effect of subtracting the stellar PSF estimate
from the data (e.g., Pueyo, 2016; Ren et al., 2018). This effect, known as
“oversubtraction”, varies both spatially and spectrally (e.g., Pueyo et al.,
2012; Betti et al., 2022), and thus introduces significant barriers for disk
studies. Model Constrained RDI (MCRDI; Lawson et al., 2022) remedies this by
suppressing the impact of circumstellar flux during construction of the
stellar PSF model using a synthetic model of the circumstellar scene.

To optimize a model of circumstellar signal (CSS) for use in MCRDI, a standard
forward modeling procedure is used (as in Lawson et al., 2023). In this
approach, an initial (unconstrained) RDI reduction of the data is carried out,
with each trial estimate of the CSS then being forward modeled and compared
with these residuals to determine goodness-of-fit. For this purpose,
optimization (of both the CSS and starlight models) considers only the region
with stellocentric separations $0\leq r<2\farcs 34$. PSF-subtraction is
performed by determining the linear combination of the available reference
frames that minimizes the (squared) residuals within this region. For all
filters, the reference integrations at each dither position are median
combined to yield a PSF library of five reference images.

We modeled the CSS as the superposition of a planet and a simple ring-like
disk shown in Fig. 2. To model the disk, we adopted a simplified version of
the GRaTer code (Augereau et al., 1999) as implemented in Vortex Image
Processing (VIP; Gomez Gonzalez et al., 2017) and assumed a scattering phase
function (SPF) that is a linear combination of two Henyey-Greenstein (H-G)
SPFs (Henyey & Greenstein, 1941). We varied disk parameters including
inclination, position angle (PA), fiducial radius ($r_{0}$), scale height, the
radial density power law indices interior and exterior to $r_{0}$, the disk
flare, the vertical density exponent, the two H-G asymmetry parameters for the
SPF, and the relative weight of the first H-G component. Following the
procedure of Lawson et al. (2023), for each roll angle of the data, we rotated
a given raw disk model to the appropriate position angle, multiplied it by a
coronagraphic transmission map, and then convolved it with a grid of
(normalized) spatially-varying PSFs drawn from WebbPSF-ext (Leisenring, 2021;
Girard et al., 2022). Throughout this process, we accounted for the position
of the coronagraph center in each exposure after alignment (see Section 2.4).
The grid of PSFs samples the origin as well as five logarithmically spaced
radial positions at each of four linearly spaced azimuthal positions (for 21
PSF spatial samples in total).

Nominally, the planet model component in the CSS would be governed by three
parameters: two for position and a third for brightness. To simplify the
optimization, we determined the planet’s position a priori using the F444W
data. To do this, we first optimized the disk model while masking the region
within $5\lambda/D$ of the planet’s estimated position. Subsequently, we tuned
the planet’s position by fitting the disk-model-subtracted residuals with PSFs
drawn from WebbPSF (Perrin et al., 2014). The planet is then fixed to this
position for all filters, reducing the planet model component to a function
only of planet brightness. The overall CSS model is scaled post facto to
minimize the residuals with the data. As such, the planet brightness parameter
effectively governs the brightness of the planet relative to the disk.

We optimized the composite CSS model of 12 parameters separately for each
NIRCam filter using the modified Powell algorithm as implemented in the LMFit
package (Newville et al., 2022). For the short wavelength channel data only
(F182M and F210M), we median combined the integrations for each exposure for
the optimization process. This significantly reduces run time, as these
filters’ exposures have both more detector pixels and more integrations than
the long wavelength filters. Once we identified an optimal CSS model, we
performed MCRDI using that model to mitigate oversubtraction but otherwise
using identical PSF subtraction settings as for the initial unconstrained
reduction. Specifically, we construct the model of the stellar PSF by
subtracting the optimized CSS model from the data when computing the optimal
linear combination of the reference images. The resulting stellar PSF model is
then subtracted from the original data. An example of the best-fit MCRDI model
components is shown for the F335M data in Figure 2. In order to create an
uncertainty map, we computed uncertainties for the stellar PSF model using
standard error propagation for a linear combination of the reference images,
and propagated uncertainties for subtraction of the PSF model from the science
data. Finally, we used the standard uncertainty propagation for a median
combination to get uncertainties for the final PSF-subtracted science images.

Figure 2: An example illustrating the components of the MCRDI model for the
F335M data. Left: the reduced MCRDI image. Middle left: the best-fit PSF-
convolved disk model. Middle right: the composite disk and planet model.
Right: the final residual map for this model. The field of view is chosen to
show the model within approximately the $r<2\farcs 34$ optimization region
used for MCRDI. As in Fig. 1, the images were rotated by 60.9∘ to align the
main disk with the horizontal axis.

We note that the MCRDI process involves fitting models to the disk for use in
avoiding oversubtraction, but we do not here attempt to interpret these as a
rigorous or physically-correct model of the disk, as this model framework is
intentionally simple.

For both the MCRDI and classical subtraction methods, the subtracted images
from each roll angle are rotated and averaged together in the usual fashion.
Figure 1 shows the $\beta$ Pic system as seen in all filters after these PSF
subtractions rotated 60.9∘ (given a position angle of 29.1 ∘ from Apai et al.,
2015) to align the main disk horizontally.

###  2.6 Deconvolution

We applied a deconvolution step to remove the blurring effect of the complex
NIRCam and MIRI coronagraphic off-axis point spread functions, which arise due
to the diffractive effects of the coronagraph Lyot stop pupil masks. To
deconvolve the data, we adopted the Richardson–Lucy (R–L) algorithm
(Richardson, 1972; Lucy, 1974) with a small alteration to accommodate
coronagraphy. Specifically, for each iteration of the R–L algorithm, before
convolving the prior iteration with the PSF, we multiplied it with a
coronagraphic transmission map. Effectively, this includes the coronagraph’s
effects in the current estimate of the blurring at each step, providing better
performance near the coronagraphic IWA. Though there are more recent methods
designed to accommodate spatially varying PSFs (e.g., Sakai et al., 2023), we
find that using a single PSF is sufficient for these NIRCam and MIRI Lyot data
so long as the coronagraphic transmission is also considered as described.
While the shape of the PSF for these data does change with detector position,
the changes are sufficiently small that such an approach is not confounded;
Appendix A demonstrates this through application to a synthetic disk model
convolved with consideration for these spatial variations, where surface
brightness measurements for the deconvolved result match those of the ground-
truth (unconvolved) disk model to within $\lesssim 1$% beyond the IWA and
$\sim 10$% within the IWA.

We note that for the MIRI F1550C data, the 4QPM introduces significant spatial
variations in the PSF. Moreover, suitably precise coronagraphic transmission
maps for the 4QPM are not yet available. As such, our deconvolution approach
is least well-suited for these data and induces significant artifacts,
particularly at the edges of the field of view and the 4QPM boundaries. The
deconvolved F1550C image of Figure 3 is therefore provided as a work-in-
progress result and is not the sole basis for any morphological analysis
herein.

Since bright circumstellar flux persists to the edges of the field of view for
both MIRI filters, the deconvolution approach can introduce edge artifacts
that increasingly infringe on the region of interest over successive
iterations. To address this, prior to deconvolution we padded each initial
MIRI image with a 35 pixel wide border of NaNs and then iteratively replaced
the NaN values with the median of values within a $5\times 5$ pixel window
until none remained. Following deconvolution, this padded border is cropped
from the final image. This approach is effective in preventing the
aforementioned edge artifacts and induces no appreciable difference in the
regions free of artifacts, though some small effect should be expected within
a few $\lambda/D$ of the image edges.

Uncertainties for the deconvolved products are the same as those for the
convolved images.

Figure 3 shows the deconvolved disk images, on the same spatial scale and
extent as Figure 1.

Figure 3: Same data as Figure 1 after deconvolution of the off-axis
coronagraphic PSFs as described in section 2.6. The presence of a vertical
warp from the inclined secondary disk is seen more clearly in the deconvolved
NIRCam data. Many complex asymmetric substructures are visible in the
deconvolved MIRI data, consistently detected in both filters. Note that the
deconvolution of the F1550C data is more challenging and thus has higher
resulting uncertainties as discussed in the text. Also the masked region
around the F2300C Lyot occulter and its supports has been expanded to block
edge regions that are not reliably deconvolved. Similarly to Fig. 1, the
images were rotated by 60.9∘ to align the main disk horizontally.

##  3 Disk Morphology

In this section we quantify the morphology of the observed asymmetries and
surface brightness in order to facilitate comparisons with literature data and
the model described in Sect. 5. Figure 4 shows the main features discussed in
this section.

###  3.1 Two misaligned disk components

The JWST observations detect the presence of the secondary disk component,
clearly observed in the MIRI images, both before and after deconvolution, and
also present in the NIRCam images, though significantly fainter (see Figs. 1
and 3). In the NIRCam data, the appearance of the secondary disk is broadly
similar to prior Hubble (e.g. Heap et al., 2000; Golimowski et al., 2006) and
ground based (Milli et al., 2014) observations: slightly inclined with respect
to the main disk, superimposed on it within the innermost $\sim 5$ ″, and
roughly symmetric between the northeast (NE) and southwest (SW) sides. In the
MIRI data, particularly the F1550C, within the same few arcseconds the disk
appears more vertically extended or thicker, plausibly due to thermal emission
from the inclined secondary disk.

Figure 4: Complex nebulosity and substructures in the disk, as seen in F2300C.
The stretch is adjusted to bring up fainter features observed in the outer
regions. Relevant features are annotated including the main and secondary disk
planes and cat’s tail. The nebulosity to the west (top) includes both
irregular clumps or knots, and also some apparent linear striations or
tendrils. To the east (bottom) the clumps are fainter and fewer but still
present, including one particular bright clump at bottom center. These clumps
are consistently detected in identical positions on the sky seen in the two
distinct rolls, confirming they are astrophysical in origin rather than any
instrumental artifact.

At 15 and 23 µm, the inner secondary disk appears to connect outwards to two
distinctly different asymmetric features extending on either side. On the NE
side, a long linear extension of the secondary disk extends out to at least
18″($\sim$ 340 au, the edge of the MIRI Lyot field of view). This results in a
bifurcated or forked appearance of the disk to the NE side, such as has been
previously observed for a few other disks777For instance HD 111520; see figure
4 of Crotts et al. (2022).. On the SW side, instead we see a curved filament
of material seemingly streaming off of the disk at a high angle ($\sim$45∘).
On account of its raised, slightly S-curved appearance we refer to this
feature as the “cat’s tail”, described in more detail in the following
subsection and Sect. 5.

####  3.1.1 Position angle of the main disk

In order to better analyze the location and angle of the secondary disk, we
performed a profile-fitting analysis similar to the one described in
Golimowski et al. (2006). We investigated the projected vertical profile in
the deconvolved images at different distances from the star, to extract a
profile that traces both the main and secondary disks. First we measured the
position angle (PA) of the main disk’s midplane such that we can rotate the
image and align the midplane of the main disk with the horizontal. We put the
images through a high pass Gaussian filter with a 3 pixel sigma. We measured
the average value of all the pixels in the high pass filtered image, and
selected only those above with fluxes 10 times or larger than the average in
order to include only high SNR points. We then fit a line to the main plane of
the disk, obtaining slightly different PAs for each filter, listed in Table 2.
Uncertainties were obtained from error propagation of the reduced images and
the line fitting process.

Given the significant offset between NIRCam and MIRI we investigated an
alternative method to measure disk PA. To do so, we sequentially rotated the
images over a range of possible PAs and fit a single Lorentzian component to
vertical cuts in the disk. We then plotted the location of the peaks of the
fitted function vs the distance to the star. That gives a line that we then
use to validate our initial rotation angle vs the PA. If the rotation angle
does _not_ correspond to the average disk PA, then the location of the peak of
that best fit single component appears diagonal (say above the horizontal line
on the NE and below on the SE). If the rotation angle does correspond to the
average disk PA then the location of the single component peak is centered on
a horizontal line —albeit with some oscillation due to the presence of a
second component, see below— and appears symmetrical from a NE vs SW. More
formally, the measured disk average PA is defined as the rotation angle that
minimizes the absolute difference between the NE and SW locations of the peak
of a single component model.

Using this method we find results broadly consistent with Table 2, and again
indicating variation in PA with wavelength: $29.1$ $\pm$ $1.1\arcdeg$ for the
NIRcam short wavelength channel, $29.6$ $\pm$ $1.4\arcdeg$ for NIRcam long
wavelength channel, and $30.4$ $\pm$ $1.5\arcdeg$ for MIRI filters. The PA we
measure in the NIRCam data is consistent with the value found by prior authors
($29.1\pm 0.1\arcdeg$ in Apai et al., 2015, $29.0\pm 0.2\arcdeg$ in Milli et
al., 2014), although we note the best fit PA found in $L^{\prime}$ SPHERE
observations in Milli et al. (2014) is 30.8∘. Literature reporting on mid-IR
observations (e.g. Li et al., 2012; Telesco et al., 2005; Han et al., 2023)
suggest a PA $>$30∘. The MIRI data are consistent with that, supporting the
inference that the apparent major axis PA does vary with wavelength. At least
a partial explanation to this can be found in Skaf et al. (2023), where they
show in their Fig. 6 how the best PA varies for the two different sides of the
disk. Besides the IR, the gaseous CO disk in Dent et al. (2014) was also
inclined $\sim$ 4 ∘ with respect to the dusty disk, observed in the continuum.

Due to the variation in apparent major axis position angle between filters,
for consistency of image display we decided to use the $29.1\arcdeg$ angle
obtained in the NIRCam short-wavelength channel analysis for the de-rotation
of the images in all figures where we are depicting the major axis horizontal.
Data analyses however used the individual best-fit position angles per each
filter.

####  3.1.2 The secondary disk

After rotating the images to align the main disk with the horizontal, we
binned the data in the horizontal direction to increase the SNR of vertical
cuts. For the NIRCam data we adopted a horizontal bin width of four pixels
while for the MIRI data we adopted a bin width of two pixels. We then
performed a fit to each vertical cut using two Lorentzian curves
simultaneously: the first one to match the position of the main disk, forcing
the peak position to be within a few pixels of the main plane of the disk, and
the second one free to trace the location of the secondary disk/fork. Results
for the peak positions as a function of circumstellar distance for both
components are shown in Appendix B, Figs.20 & 21.

We performed a fit for the secondary disk in both NIRCam and MIRI, including
ranges from -100 to 100 au, where the disk is detected. The angle of the
secondary disk with respect to the main disk for the different filters are
summarized in Table 2.

The inner secondary disk seen using MIRI connects with the asymmetric
extension creating a “fork” on the NE side, and the “cat’s tail” on the SW
side, as shown by Fig.4 and 21. In order to further characterize these
observed features we extracted their inclination angles with respect to the
main disk. We distinguished two different ranges in the SW side for the MIRI
images: (1) base of the cat’s tail, from 120 au until 190 au; and (2) tip of
the cat’s tail, from 190 au until 230 au. The angles with respect to the main
disk of the two regions can be found in Table 2.

| Main disk∗ | Components of secondary disk†
---|---|---
Filter | PA | Sec. | NE | Cat’s tail | Cat’s tail
|  | disk | fork | base | tip
| (deg) | (deg) | (deg) | (deg) | (deg)
F182M | 29.6$\pm$ 1.5 | 5.8 | – | – | –
F210M | 29.6$\pm$ 1.1 | 5.9 | – | – | –
F250M | 28.0$\pm$ 1.2 | 6.2 | – | – | –
F300M | 28.0$\pm$ 1.9 | 5.9 | – | – | –
F335M | 28.0$\pm$ 1.9 | 5.0 | – | – | –
F444W | 28.9$\pm$ 1.7 | 5.1 | – | – | –
F1550C | 31.6 $\pm$ 1.0 | 5.9 | 5.5 | 45.0 | 21.5
F2300C | 33.2 $\pm$ 1.9 | 5.6 | 4.8 | 44.7 | 21.1
Table 2: PA of the main disk, and angles measured for the different components
of the secondary disk. ∗PA angles of the main disk refer to the
counterclockwise angle measured from the north vector. † Angles of the
components of the secondary disk refer to the relative location with respect
to the main disk.

Overall, the secondary disk is prevalent in all images on the NE side, but
harder to detect in the SW side. We also see a dependence with wavelength in
the NIRCam images, where the secondary disk becomes fainter at longer
wavelengths, indicating the grains are efficient at scattering light at visual
and short NIR wavelengths. As the thermal emission starts to dominate, the
secondary disk then picks up again in brightness at MIRI wavelengths, as well
as becoming more spatially extended to the NE creating the appearance of a
fork.

###  3.2 The cat’s tail in $\beta$ Pic

The most conspicuous feature observed for the first time in this JWST dataset
is the prominent curved structure on the SW side of the disk, resembling a
“cat’s tail”, seemingly attached to the disk. We observe it in both MIRI
filters, at 15.5 $\mu$m and 23 $\mu$m, and it is present in both rolls.

The structure extends from $\sim$ 6.5 ″(from the location of the star) to
$\sim$ 13.5 ″, corresponding to $\sim$125 to 265 au. Details of how the
structure were characterized are in the previous section, Table 2 and Fig. 21.
The base of the cat’s tail separates from the disk near a projected separation
of $\sim 6$ ″(115 au) and extends for $\sim 3$ ″(60 au), curving away from the
main disk and reaching a position angle of $\sim$45 ∘ with respect to the main
disk midplane (see Table 2). The tip of the cat’s tail is considerably
shorter, $\sim$1.5 ″(30 au), and shows a very similar inclination in both
filters of $\sim$21∘.

The cat’s tail has not been reported before, despite deep HST observations
with several instruments (see e.g. Burrows et al., 1995; Golimowski et al.,
2006; Apai et al., 2015; Ballering et al., 2016). There is also no evidence
for such structure in ALMA data (Dent et al., 2014; Matrà et al., 2019), nor
any of the ground-based mid-IR observations (e.g. Telesco et al., 2005; Li et
al., 2012; Skaf et al., 2023). There may be faint hints of material near this
location in the outermost part of the $L^{\prime}$ band imaging presented in
Milli et al. (2014), but only at the edge of the FOV where PSF subtraction
residuals make any potential detection highly uncertain. Nor do the JWST
NIRCam observations show any sign of the cat’s tail, setting $3\sigma$ upper
limits of 0.23 and 0.01 mJy/arcsec2 in F210M and F444W respectively for any
excess emission at the location of the center of the tail.

In the upcoming sections we analyze the color (Sect. 4) and possible origin
(Sect. 5) of the cat’s tail. Additionally, further MIRI observations have
already been performed at the time of submission of this paper for the program
GTO 1241 (PI Choquet), including the filters F1140C and F1065C. These
observations will help unveil more details about the nature of the cat’s tail.

###  3.3 The nebulosity in the west side of the disk

The MIRI F1550C and F2300C images show in the west region of the FOV (above
the midplane of the main disk as shown in Figs. 1 and 3) faint small
structures that we will refer to as “nebulosity” from here on. The eastern
region of the image (below the midplane of the main disk), also shows a few
fainter such features. Figure 4 shows the F2300C deconvolved image with
relevant features annotated and a stretch that favors the observation of the
nebulosity. These structures are not observed to have any symmetry axis with
respect to the main disk or the star, and apparently, they do not have a
preferential orientation. Some of this nebulosity has a striated or tendrilous
appearance, in particular a few parallel linear clumps above the “cat’s tail”
region, and a few rough lines to the north and not quite parallel with the
disk plane.

We confirm these features are not instrumental artifacts or residual stellar
speckles. These features are observed consistently in both rolls for both MIRI
filter observations (i.e., four independent exposures, using two completely
different sets of coronagraph optics within MIRI). The fact that these
features are fixed to the plane of the sky, independent of observatory
position angle, unambiguously confirms they are astrophysical in origin. Also
there is no sign of such nebulosity present in the reference star images. In
Sect. 5 we provide a hypothesis of its possible origin as part of the $\beta$
Pic system, and its connection with the secondary disk and the “cat’s tail”.

We cannot yet strictly rule out that some portion of these clumps may be mid-
IR background sources in chance alignment; however the absence of any shorter-
wavelength counterparts in the HST (e.g. Golimowski et al., 2006; Apai et al.,
2015) and NIRCam data rules out most classes of potential astrophysical
background object. The proximity to $\beta$ Pic of so many background sources,
all of which are so extremely red as to be only visible at $\geq 15$ microns,
would be very low probability. Further the apparent alignment of the few
linear streaks near the cat’s tail and others roughly aligned to the disk
plane suggests physical association of some sort. We anticipate that proper
motion analyses using multiple epochs of MIRI observations will conclusively
rule out background status for all or nearly all of this nebulosity.

###  3.4 Surface brightness

Measurements of the surface brightness of the disk can reveal asymmetries in
the dust distribution, and provide an overall view of the shape and flux of
the disk with respect to the position of the star.

In order to investigate the surface brightness, we used the deconvolved images
shown in Fig. 3 with the main disk major axis rotated horizontal according to
the measured PA. We perform an integration in the vertical direction (i.e.,
perpendicular to the main disk) without binning. We obtained values for the
disk surface brightness along the whole FOV with the exception of the location
of the star. The IWA varies for the different instruments and modes, according
to the technical documentation 888See technical details in: https://jwst-
docs.stsci.edu/jwst-mid-infrared-instrument/miri-observing-modes/miri-
coronagraphic-imaging and https://jwst-docs.stsci.edu/jwst-near-infrared-
camera/nircam-observing-modes/nircam-coronagraphic-imaging. This is reflected
in the measured SNR when extracting the profiles. We chose to only report
measurements for which SNR$>$15 for all datasets. As a result, for NIRCam we
can go down to 0.2″ and still retrieve some structure, while for MIRI the IWA
goes up to 2.7″ for the Lyot stop at 23 $\mu$m (4.6$\lambda$/D), in agreement
with the expected performance (Boccaletti et al., 2015, 2022). For the 4QPM at
15.5 µm we excluded a larger portion of the inner disk, due to the PSF
artifacts present even after deconvolution, leaving an effective IWA of 1″,
instead of 0.5″. The cutoff limit for the detection of the disk’s outer edge
was set to SNR$>$5.

We present the surface brightness for both the NE and SW sides, the contrast,
and contrast adjusted for the projected stellocentric distance in Fig. 5.

The observed surface brightness decreases with wavelength for the NIRCam
filters, as the grains become less effective to scatter the stellar light. But
as the thermal emission picks up in the MIRI filters, the surface brightness
of the disk starts to increase with wavelength (top left and right panels).
Contrast levels also increase significantly with wavelength, as the star
becomes less bright middle panels; and when corrected for the effect in flux
of projected stellocentric distance (r2), the contrast remains constant for
NIRCam data from 30 up to $\sim$150 au, where the main disk becomes fainter
bottom panels. The distance corrected contrast for MIRI, however, shows an
abrupt decrease after $\sim$100 au in the SW side (bottom left panel),
specially for the 15.5µm filter. This can also be observed in the images (e.g.
Fig.3), as the cat’s tail picks up. It’s worth noting as well, that the
surface brightness profiles for NIRCam filters show a small excess in the NE
side at $\sim$10 au likely due to the presence of $\beta$ Pic b.

Figure 5: Surface brightness vertically integrated profiles of the disk in the
different JWST filters, as labeled. Only regions with SNR$>$15 were
investigated. Top panels show the surface brightness versus stellocentric
distance to the star for the SW and NE side left and right,respectively.
Middle and bottom panels show the contrast to the stellar flux and the
contrast corrected for the projected distance dependence for the SW (left) and
NE(right) sides. Shaded grey areas correspond to 3$\sigma$ range, color coded
for each filter. The vertical dashed line in second panel (NE) indicates the
location of $\beta$ Pic b.

When comparing the surface brightness of both sides of the disk, we find that
the surface brightness levels along the disk are similar, yet the flux in the
SW side decays faster than in the NE side after $\sim$100 au in all filters.
This shorter extension of the SW side is well established in the literature
(see e.g. Fig. 8 in Janson et al., 2021). There is no clear evidence of large
asymmetries or brightness peaks/clumps in the NIRCam data.

However in both MIRI filters there is significant spatial variation between NE
and SW sides, and presence of localized peaks or clumps in brightness between
$\sim$40 and $\sim$100 au, shown in Figs. 6 and 7. This asymmetry, with larger
flux in the SW side of the disk, seems to be co-located with the CO and dust
clumps previously reported in the literature (e.g. Telesco et al., 2005; Dent
et al., 2014; Cataldi et al., 2018). We observe some structure in the
asymmetry, with one brighter peak, close to 50 au on the SW side, and at least
a second peak at $\sim$80 au with a significantly lower brightness. In order
to further test the connection of the observed asymmetry to previous
observations, we estimate the location of the brightest “clump” through a
simple Gaussian fit (see Fig. 6, dashed lines). We obtain peak locations at
55.5$\pm$3.9 au for F1550C, and 63.4$\pm$2.9 au for F2300C, where
uncertainties are from the fit results using the per-pixel uncertainty maps
from the data reduction process as weights. In the F2300C data, the observed
peak location is not far outside of the area blocked by the large Lyot
coronagraph occulting spot; the apparent difference in peak location between
the two filters may arise at least in part from instrumental systematics. In
order to compare with previous measurements in the literature, we include in
Fig. 6 the central locations measured in Skaf et al. (2023) and Han et al.
(2023). Fig. 7 top panel compares the MIRI F1550C data to the locations of the
four clumps C1-C4 reported by Skaf et al. (2023) based on 12 µm observations
from the ground (their C1 corresponds to the main brightness clump as detected
by several prior works). The C1 and C2 clumps are clearly detected in the MIRI
15 µm data. The other two are at sufficiently smaller separations that PSF
subtraction and deconvolution residuals make identification difficult. The
presence and nature of the mid-infrared brightness asymmetries is further
discussed in Sect. 6.

Figure 6: Difference between the surface brightness of the SW and NE sides of
the disk. Asymmetric peaks or clumps are detected with MIRI between 40 and 100
au, with both the peak positions and brightnesses varying significantly
between 15 and 23 µm. Dashed lines show the Gaussian fit used to estimate the
central location of the brightest peak in the asymmetry at each wavelength.
Shaded grey areas correspond to 3$\sigma$ ranges. Color coded shaded areas
correspond to results from prior observations at 10–12 µm. No significant
variation between SW and NE sides is seen in the NIRCam data for this range of
separations. Figure 7: The MIRI coronagraphy data in F1550C and F2300C (15.5
and 23 µm) presented with a linear stretch and scaled by the projected
distance from the star squared. This stretch better visually emphasizes (a)
the brightness asymmetry between the two sides, and (b) the presence of
localized peaks or clumps of brightness particularly to the SW. The white
arrows indicate the estimated clump locations from the Gaussian fitting shown
in Fig. 6, at 55.5$\pm$3.9 au for F1550C, and 63.4$\pm$2.9 au for F2300C. The
cyan annotations C1-C4 indicate the angular separations of the clumps reported
in ground-based VLT/NEAR 12 µm data (Skaf et al., 2023). Both C1 and C2 are
clearly detected in the MIRI F1550C data.

##  4 Disk Spectral Energy Distribution and Colors

Spectral coverage across the near- and mid-IR enables color measurements that
can help place constraints on the physical properties of disk constituents. In
this section we explore the differences in the colors along the disk.

The main disk of $\beta$ Pic does not present any conspicuous color
differences between the NE and SW sides with the notable exception of the
difference in the location of the dust clump (Sect. 3 and Fig. 6). However, if
we compare the secondary disk, particularly in the MIRI filters, we see that
it is significantly bluer than the main disk. Figure 8 combines F1550C and
F2300C data in a false color image, in which the contribution of the secondary
disk stands out as bluer. Even though some instrumental effects remain in the
F1550C image (annotated in the figure), the nebulosity and the cat’s tail also
exhibit a much bluer color than the main disk, possibly indicating a higher
temperature of the grains. Figure 9 shows the peak flux at different
stellocentric distances for the cat’s tail. In order to avoid contamination
from the main disk, we started the measurement of the cat’s tail at distances
larger than 5.2″($\sim$ 100 au), where the vertical separation between peak
flux of the cat’s tail and the main disk is $\sim$1 ″at 15.5µm. We note that
some contamination from the vertical component of the main disk might still be
present, but removing it completely would require precise and extensive
modeling, which is out of the scope of this paper. The left panel of Fig. 9
shows the ratio of the peak flux of cat’s tail to the main disk; and the right
panel shows the absolute peak flux. The beginning of the cat’s tail is marked
in both panels with a red shaded area. The cat’s tail is much brighter with
respect to the main disk at 15 µm, and seems to increase with distance until
peaking at $\sim$210 au, then decreasing until the detection is no longer
significant. In contrast at 23 µm the cat’s tail is fainter than the main
disk, and the the ratio is approximately constant with separation. In the
right panel, an increase in flux is observed in both filters around $\sim$130
au. Beyond that, both filters show a different behaviour of the flux with
distance, with an exponential decrease in F1550C and a quadratic decrease in
F2300C. Interestingly, there is an increase both in total flux (right panel)
and contrast (left panel) at $\sim$200 au, just after the change in the
orientation angle, which may indicate there is some internal structure in the
cat’s tail itself.

Figure 8: False color image generated from the F1550C (15.5 µm, cyan) and
F2300C (23 µm, red) images. At these wavelengths, the main disk is redder,
while the secondary disk, cat’s tail, and extended nebulosity are bluer. The
several linear striations or tendrils just above the cat’s tail are clearly
seen, and consistently blue. As noted in the text, the deconvolution of the
F1550C is challenging and leaves residual artifacts along the 4QPM boundaries,
as marked.  Figure 9: Cat’s tail surface brightess in both MIRI filters. Left
panel shows the contrast with the main disk, where the cat’s tail at 15.5 µm
(black line) is consistently brighter, sometimes closer to 2 orders of
magnitude. Right panel shows the surface brightness of the cat’s tail, higher
at 23 µm (grey line), and decreasing constantly with the distance to the star.
The red shaded area marks the beginning of the cat’s tail feature.

###  4.1 Spectral energy distribution at different stellocentric distances

We constructed SEDs, normalized in contrast relative to the star, combining
the measurements from scattered light and thermal emission for different
stellocentric distances, shown in Fig. 10. Starting from the disk contrast
profiles obtained in Sect. 3, the measurements for each individual distance
were made averaging three pixels centered in each of the locations for each of
the contrast profiles. For the full disk measurement, we integrated the full
contrast profile. In scattered light, the disk brightness relative to the star
decreases with wavelength from 2 to 4 µm (i.e. the disk-scattered light is
blue relative to the star), reaching a minimum in contrast around 4–5 µm, and
then rapidly brightening at longer wavelengths as thermal emission rises.

The SED shows that the colors are consistent in both sides of the disk, and
the spatial extent of the NE side is slightly larger, as we retrieve an
additional data point at 250 au.

Figure 10: Contrast spectral energy distribution for thermal and scattered
light in both sides of the disk. Data points have been obtained as explained
in Sect. 4.1. Different colors and symbols refer to different positions on the
disk, as labeled. Data below 2 µm, corresponding with NIRCam data, shows a
decrease in the scattered light contrast; while 15.5 µm and 23 µm,
corresponding with MIRI data, show how the thermal emission picks up as the
brightness from the star decreases.

###  4.2 No evidence for water ice

Water ice, a fundamental ingredient for planet formation, is expected to
sublimate in the inner regions of protoplanetary disks, and therefore be
present only at long stellocentric distances in debris disks and mature
planetary systems, past the so called snowline. Around a $\sim$1.8M☉ star,
such as $\beta$ Pic, the snowline is expected to be at $\sim$7 au (based on a
simple black body temperature estimation for crystalline water ice). However,
due to UV photo-dissociation, the evaporation of water ice can extend to
100-150 au (Grigorieva et al., 2007).

The set of NIRCam filters used was specifically chosen to provide coverage of
the 3 $\mu$m water ice feature. This allows us to do an estimation of whether
there is a significant amount of water at different distances from the central
star. Following Sect. 4.4 of Kim et al. (2019), we compared the surface
brightness values obtained in and out of the expected feature range. Given we
do not have observations at the wavelengths they selected, 2.8 and 3.2 µm, we
chose filters at 2.5 and 3.35 µm, which are still expected outside and inside
of a hypothetical water ice feature, respectively. Additionally, we repeated
the calculation for 2.5 and 3.0µm, in order to verify our results. Simply
calculating a ratio of surface brightness from Fig. 10 at said wavelengths, we
obtain values between 1.52$\pm$0.10 and 1.89$\pm$0.52 for distances between 50
and 150 au. Uncertainties were calculated from the uncertainty maps
corresponding to the images obtained from the data reduction, and the
corresponding error propagation.Comparing to Fig. 16 in Kim et al. (2019), we
see no evidence of water ice inside of 150 au. Given our method is purely
based in photometric measurements, future IFU observations will provide a more
precise determination (e.g. GO 1563, PI Chen).

##  5 Sporadic dust production as a source for the disk’s substructures

The MIRI coronagraphic data reveal several features in the debris disk,
described in detail in Sects. 3 & 4: the “fork,” which appears from a radial
extension of material on the NE side of the secondary disk; the “cat’s tail,”
which appears as a spatially extended, curved stream of material on the SW
side of the disk; and nebulosity, especially to the west. These features
exhibit unique traits, some not seen before in other disks:

  * •

The fork appears only to the NE side of the disk, suggesting an eccentric
secondary disk.

  * •

The cat’s tail appears to curve away from the secondary disk, reaching an
apparent position angle of $\sim 45^{\circ}$ with respect to the main disk,
then has a “crook” near the tip that bends back toward the main disk.

  * •

Over the extent of the cat’s tail the surface brightness changes by a factor
of $\sim 5$ from base to tip, similar to the main disk (see Fig. 9).

  * •

The nebulosity appears tendrilous, with apparent arcs of material curving off
of the disk and tail.

  * •

As shown in Figure 8, all of these features appear blue compared to the main
disk, suggesting a common origin that differs in composition or grain size
distribution and is hotter than the main disk.

These unique features motivate the need for a dynamical explanation and can
help constrain models of different dynamical scenarios. Here we examine
several possible causes. Motivated by the extended secondary disk that is seen
only on the NE side of the disk, we start by discussing a steady-state
collision point model that naturally produces eccentric disks, but ultimately
does not adequately explain the color, flux, or crook of the cat’s tail. We
then perturb this model by invoking sporadic dust production and show that it
can reproduce all of the features described above.

###  5.1 Steady state collision point model

Figure 11: Left: orbital trajectories of grains from a collision point as a
function of $\beta$ (inspired from Figure 1 of Krivov (2010)). Middle: steady-
state collision point model with $\beta<2$ producing an eccentric disk of dust
and a fan of blow-out grains. Right: collision point model with periodic dust
production events at 100 year intervals. Individual dust production events
leads to curved streams of material that travel outward while forming a
spiral, with bound grains eventually populating the main disk.

The appearance of an asymmetric extension on only one side of the observed
disk suggests the presence of an eccentric secondary disk. Indeed, Figure 8
shows a population of bluer (hotter) grains that appear to reside in this
eccentric disk. Generating an eccentric disk of dust is not dynamically
challenging—in fact, at some level, it’s unavoidable. The collision point
models of Jackson et al. (2014) predict that a single massive collision will,
over time, lead to a belt of material that continually returns to the same
“collision point,” where small dust grains are continuously produced. The
dynamics of these small dust grains are controlled by $\beta$, the ratio of
radiation pressure force to the gravitational force. The left panel in Figure
11, inspired by Figure 1 from Krivov (2010), illustrates the trajectory of
grains as a function of $\beta$ when launched from a progenitor on a circular
orbit. Bound dust with $\beta<0.5$ is immediately blown onto orbits with
larger semi-major axes and eccentricities, resulting in an eccentric,
asymmetric halo of small dust with periastron located at the collision point
(see lower panel of Figure 15 in Jackson et al., 2014). These collision point
models have been recently invoked to explain the morphology of a multitude of
debris disks seen in scattered light (Jones et al., 2023). We note that over
long timescales, collisions will eventually resorb the eccentric disk (Kral et
al., 2015). Its presence should therefore place an upper limit on the age of
the original destructive event, which we do not model in this paper.

As shown in Figure 11, unbound grains with $\beta>0.5$ leave the system
following hyperbolic or anomalous hyperbolic orbits on relatively short
timescales. When continually produced, these unbound grains produce a fan of
material leaving the disk (see upper panel of Figure 15 in Jackson et al.,
2014). The middle panel of Figure 11 shows an example of steady state dust
production from a collision point model. For this illustrative example, we
assumed a maximum value of $\beta=2$ (ultimately this is determined by the
properties of the star and dust grains) and assumed all dust production
occurred at the same point, originating from parent bodies on circular orbits.

For this steady-state scenario, we posit that the extended secondary disk to
the NE, the HST-observed inner warp/secondary disk (Heap et al., 2000;
Golimowski et al., 2006), and the CO clump (Dent et al., 2014) all belong to a
single secondary disk as suggested by Apai et al. (2015). The CO clump may
trace the location of the collision point, as suggested by Jackson et al.
(2014), which would place periastron to the SW side of the disk and the
eccentric halo to the NE, producing the observed extension of the secondary
disk on that side.

Figure 12: Orientation of the secondary disk for the steady-state model and
our favored sporadic dust production model, which each require a slightly
different inclination of the secondary disk. Note that the difference in
assumed orientation is for the secondary disk only; in both cases the main
disk orientation remains nearly edge-on with the NW side nearer to us, as has
been inferred from brightness asymmetry attributed to forward scattering (e.g.
Golimowski et al., 2006; Milli et al., 2014; Millar-Blanchaer et al., 2015).
The rotation of both disks (main and secondary) is such that the SW side’s
orbital motion is toward us, unambiguously known from radial velocity
observations of the CO gas (Dent et al., 2014) and the planet $\beta$ Pic b
(Snellen et al., 2014).

Bound grains can only explain the appearance of the extended secondary disk.
To explain the cat’s tail and nebulosity, we must invoke unbound grains. In
this scenario, blowout grains emanate from the collision point on hyperbolic
orbits that curve toward the observer. If the secondary disk were seen edge-
on, these orbits would simply appear radial, as all forces invoked so far are
radial. However, if we invoke a small inclination ($\sim 6^{\circ}$) for the
secondary disk relative to the main disk and ensure a component of the line of
nodes is within the sky plane, such that the secondary disk is not edge-on
from our perspective, then the curvature of the hyperbolic orbits can create
the observed apparent position angle of $\sim 45^{\circ}$. The nebulosity and
cat’s tail could therefore be explained as belonging to a fan of blowout
grains. Importantly, this model suggests that while the base of the cat’s tail
is at $\sim 120$ au from the star, the tip of the cat’s tail is at a physical
distance of $\sim 1000$ au. The left panel of Figure 12 illustrates the
orientation and rotation of the disk and cat’s tail in this steady-state
collision point scenario.

Figure 13: Synthetic 23 $\mu$m thermal emission image of a two-component “toy”
steady-state collision point model with $\beta=0.7$ blowout grains and without
external forcing (e.g., ISM drag), as discussed in section 5.1. The left panel
shows the projected view oriented the same as Figure 4, and the right panel
shows a face-on view looking down onto the disk from the NW. Invoking
preferential production of $\beta\sim 0.7$ via porous silicates can produce a
cat’s tail-like stream of material, but seems to be at odds with the surface
brightness profile, color, and crook of the cat’s tail.

Figure 13 shows a synthetic thermal emission image of a “toy” model of this
scenario, as applied to $\beta$ Pic. We assumed parent bodies on circular
orbits occupying a ring with 10% fractional width. At a single azimuthal
location in the parent body ring, we launched bound grains with $\beta$
uniformly and randomly distributed from $0.01$ to $0.49$. We integrated their
orbits using a leap-frog integrator for a full orbital period, periodically
recording their locations to simulate constant production. We assumed
$\beta\propto s^{-1}$, where $s$ is grain radius, and weighted them according
to a Dohnanyi size distribution (Dohnanyi, 1969). We then used dustmap (Stark,
2011) to synthesize the 23 micron thermal emission image shown in Figure 13,
adopting optical constants predicted by Mie theory for 95% porosity
astronomical silicates (Li & Greenberg, 1997). For the cat’s tail feature, we
adopted a stream of grains with $\beta\sim 0.7$. The disk is inclined by 4∘
from edge-on with the line of nodes rotated out of the sky plane by 30∘ and
has a position angle of 5∘. The right panel of Figure 13 shows the same model
face-on. We note that we arbitrarily scaled the number of $\beta\sim 0.7$
grains relative to the bound dust. In reality, this is not allowed, as the
relative surface brightness is governed by the size distribution of
collisional fragments and collision rate of the disk. Prior modeling efforts
suggest a steep drop in the size distribution at the blow-out size (Thebault &
Kral, 2019).

We identify four challenges for this steady state model. First, to explain the
cat’s tail we must invoke a near-singular value of $\beta$ for all blowout
grains. Adopting a range of $\beta>0.5$ would result in a broad fan-like
feature that does not resemble the observed stream. One potential explanation
for this is that a broad range of grain sizes may _all have the same value_ of
$\beta\sim 0.7$. Arnold et al. (2019) determined the $\beta$ value for
different grain sizes and compositions of dust around a $\beta$ Pic like star
using both Mie theory and the discrete dipole approximation method. Figure 4
in Arnold et al. (2019) shows that Mie theory predicts highly-porous silicates
have a near-constant $\beta\sim 0.7$ over several orders of magnitude in grain
size. While in theory this should create a stream of material resembling the
cat’s tail, in practice we find this unlikely, as $\beta$ is inversely
proportional to density, and thus small variations in density among grains
would diffuse this stream substantially.

Second, the stretch of our synthetic image shown in Figure 13 is substantially
larger than the data set shown in Figure 4; the base of the cat’s tail is much
brighter than the tip in our model. Put simply, porous silicate grains at 1000
au are substantially cooler and fainter in thermal emission than at 100 au.
Future studies may find a mixture of materials, which we do not investigate
here, that can mitigate this issue.

The third challenge of our model is the crook near the tip of the cat’s tail.
While our model explains the apparent curvature _away_ from the main disk, it
does not explain the observed curvature back _toward_ the main disk. To
address this, we must invoke either a difference in initial conditions, such
as a mass dependent initial velocity distribution for collisional ejecta, or
some sort of external forcing, such as ISM drag (Debes et al., 2009). Indeed,
the proper motion of $\beta$ Pic is primarily toward the north999From
Hipparcos: pm${}_{\mathrm{RA}}=4.65$ mas/yr; pm${}_{\mathrm{Dec}}=83.10$
mas/yr, potentially providing an ISM drag force that points roughly parallel
to the main disk midplane. Further, the highly-porous silicates that could
explain the selective $\beta$ production may have low enough density to be
altered by the ISM. However, the porous silicates in the cat’s tail would
occupy a broad range of grain sizes, such that ISM drag would splay the tip of
the cat’s tail, effectively size-sorting it into a fan-like structure, which
is not observed.

Finally, the cat’s tail appears to have a finite extent, with its flux
dropping rapidly at the tip. In a steady-state scenario, the cat’s tail should
extend to much greater distances and the flux of the cat’s tail should drop
smoothly until it falls below the sensitivity limits of our observations.

###  5.2 Sporadic dust production model

Here we present our preferred model that addresses all of the above
shortcomings of the previous model by invoking a simple change: sporadic dust
production, in which dust is produced at irregular intervals from irregular
locations within the disk. We start by discussing the concept of discretized
dust production events.

When produced continuously, bound dust forms a broad eccentric disk and
blowout grains form a fan, as shown in the middle panel of Figure 11. However,
we should expect that, at some level, a disk is generated by many individual
collisional or outgassing events. The right panel of Figure 11 shows the
expected geometry of a collision point model in which dust is produced
periodically in regular, discretized intervals. Here the geometry of an
individual production event appears as an arc of debris spreading outward.
Importantly, this arc is $\beta$-sorted, such that particles near the base
have $\beta\sim 0$ while particles near the tip are unbound.

To explain the cat’s tail and tendrilous nebulosity observed in the MIRI data
set, we posit a quasi-steady-state eccentric secondary disk with _sporadic_
dust production events. In our model a single recent collision does not
produce the entire secondary disk. Rather, the secondary disk (including its
extension to the “fork”) may be relatively old, the product of many evolved
dust production events over its history as predicted by the collision point
model, while the cat’s tail, nebulosity and other asymmetric features like the
CO clump are due to recent individual dust production events. While dust
production must occur to the SW on average to explain the eccentricity of the
secondary disk, we do not require a single point of production when modeling
the cat’s tail or nebulosity. Rather, we allow dust to be produced
sporadically in time and location. By doing so, we are able to reproduce the
geometry and finite extent of the cat’s tail, as well as arcs that roughly
resemble some of the observed nebulous features.

From the perspective of an observer, the geometry of a sporadic dust arc could
take on many forms, from linear features to eccentric spirals, depending on
the distribution of $\beta$ values, orbit/orientation of the progenitor, and
time since the collision. To investigate these structures, we integrated the
equations of motion for 10k particles covering a broad range of $\beta$
values, releasing them from different orbits and orientations and tracking
their evolution with time. We examined inclinations for the secondary disk
ranging from $0-10^{\circ}$ and position angles ranging from $0-6^{\circ}$
with respect to the main disk. Given that the longitude of ascending node
cannot point toward the observer, as this would produce an edge-on secondary
disk with strictly radial features, we limited the line of nodes to lie in the
sky plane—future modeling efforts could explore this additional free
parameter. In our efforts to reproduce a feature resembling the cat’s tail, we
found the following:

  * •

Assuming the secondary disk orbits in the same direction as the CO gas (Dent
et al., 2014), we must swap the near and far sides of the secondary disk, as
shown in the right panel of Figure 12. This change affects the secondary disk
only, not the main disk, and is a relatively minor change in inclination given
the nearly edge-on alignment of the disks.

  * •

Assuming the line of nodes of the secondary disk lies in the sky plane, the
original point of dust production must occur on the far side of the disk near
the secondary disk’s projected minor axis; the orbital anomaly of the
collision point is likely within $\sim 5$ degrees of the projected minor axis.

  * •

We require a maximum $\beta$ of $\sim 10$ to reproduce the curvature and
extent of the cat’s tail. These high-$\beta$ grains move outward faster, which
allows grains to reach $\sim 1000$ au in a fraction of the progenitor’s
orbital period, producing the proper curvature and extent of the cat’s tail.

The left panel of Figure 14 shows $\beta$ calculated as a function of grain
size for several commonly-invoked compositions and porosities using Mie theory
and the Bruggeman mixing rule. Dotted lines show non-porous dust while solid
lines show dust with 95% porosity. The dotted red curve shows the best fit
model to the main disk’s halo component from Ballering et al. (2016), while
blue, orange, and black lines show organic refractory (Li & Greenberg, 1997),
astronomical silicate (Draine & Lee, 1984), and amorphous carbon grains (Zubko
et al., 1996), respectively. To first order, high porosity shifts these curves
to larger grain sizes and broadens the peak, resulting in larger and more
abundant blowout grains; the fact that these grains are visible at 23 $\mu$m
suggests that they may have high porosity. To achieve $\beta\sim 10$ dust for
$\beta$ Pic, we likely must invoke refractory organics or amorphous carbon. We
emphasize that independent of color and flux, _the shape of the cat’s tail_
appears to inform the dust composition.

Figure 14: Left: $\beta$ as a function of grain size for Ballering et al.
(2016) best fit composition (dotted red) and other compositions/porosities.
Our dynamic model requires high $\beta$ and larger grains, consistent with
fluffy refractory organic material. Right: Normalized spectral energy
distributions for the best-fit Ballering et al. (2016) halo component for the
main disk (red), as well as three highly porous compositions. Only fluffy
organic refractory material can produce a bluer mid-IR color (higher ratio of
F1550C/F2300C) than the Ballering et al. (2016) composition. Solid, dashed,
and dotted lines indicate different particle size ranges and/or locations, as
labeled in gray.

In addition to requiring a high $\beta$ value, the dust must appear bluer than
the main disk. The right panel of Figure 14 shows the spectral energy
distributions for highly porous (i.e., “fluffy”) grains around $\beta$ Pic
calculated via Mie theory. Solid lines correspond to large grains (50–500
$\mu$m) near the base of the cat’s tail (100 au), while dashed lines
correspond to small grains (0.1–5 $\mu$m) near the tip of the cat’s tail (1000
au). The dotted red line shows the best fit composition to the main disk halo
from Ballering et al. (2016), valid only for small grains at 100 au. Vertical
gray bars mark the approximate bandpasses of the F1550C and F2300C filters.
Only one of the compositions (blue curves) allows the dust to appear bluer
than the main disk (red curve): organic refractory material. We therefore
focus solely on fluffy organic refractory material for our sporadic dust
production model.

Prior to modeling multiple sporadic events, we first focus on reproducing the
single most prominent feature: the cat’s tail. Figure 15 shows a synthentic 23
$\mu$m thermal emission image of a two-component sporadic dust production
model without any external forcing. To produce the eccentric disk component,
we simulated bound grains originating from a single progenitor on a circular
orbit at 85 au. We adopted $0.01<\beta<0.49$ and integrated the equations of
motion for an orbital period, calculated their grain size from the curve shown
in Figure 14, and weighted them according to a Dohnanyi size distribution at
the collision point. We repeated this process for the cat’s tail component,
with $0.5<\beta\lesssim 10$, and adopted a steeper-than-Dohnanyi size
distribution $dN/ds\propto s^{-3.8}$ for the ejecta. This size distribution
better reproduces the flux over the extent of the cat’s tail and is consistent
with the aftermath of recent dust production event that has not yet reached
collisional equilibrium (Leinhardt & Stewart, 2012). For the cat’s tail, we
launched grains from an orbital anomaly of $5^{\circ}$ to the NE of the
secondary disk’s apparent minor axis and stopped the orbital integration after
150 years. We adopted an inclination of $6^{\circ}$ from edge-on and a
position angle of $2.5^{\circ}$ for both components.

Figure 15: Synthetic 23 $\mu$m thermal emission image of a two-component model
without external forcing, as proposed in section 5.2, consistent of a quasi-
steady-state secondary disk and a single 150 yr old collisional event. This
reproduces the shape of the cat’s tail, with exception to the crook near the
tip. The left panel shows the projected view oriented the same as Figure 4,
and the right panel shows a face-on view looking down onto the disk from the
NW.

As shown in Figure 4, the cat’s tail features a crook near the tip, in which
it bends back toward the main disk plane. To explain this crook, we must
include additional physics in our model. We see two possible options. First,
we could include external forcing primarily directed toward the SW, such as an
ISM drag component (Debes et al., 2009). Alternatively, we could posit that
the collisional or outgassing event produced a size-dependent velocity
distribution in the SW direction, such that small grains near the tip of the
cat’s tail received a much larger velocity kick than their larger
counterparts. Size-dependent ejection velocity is a known outcome of cometary
outgassing, as small grains are quickly entrained by the gas tail and
accelerated at a higher rate until the gas tail dissipates (e.g., Moorhead,
2021). Here we focus on the former, ISM drag, as an example, though we do not
rule out alternative explanations.

Figure 16 shows a synthetic 23 $\mu$m thermal emission image for our preferred
sporadic dust production model including ISM drag. We adopted the same
parameters as our previous model, but added ISM drag, which extends the tip of
the cat’s tail to the SW. Here we have included multiple sporadic dust
production events that qualitatively reproduce some of the geometry of the
observed nebulosity and arc-like features seen in the MIRI data, as marked in
Figure 4. These fainter “tendrils” cover a broad range of times since
collision and are all generally created near the projected minor axis, which
notably corresponds to roughly where the primary and secondary disk planes
cross. Events consistent with the SW nebulosity occur on the far side of the
disk and are $\sim 150-200$ years old, while events consistent with the arc-
like features to the NE occur on the near side of the disk and are $\gtrsim
1000$ years old. We note that all ages are approximate, as they are somewhat
degenerate with the orbital period of the parent body and $\beta$ values for
the dust.

The ISM drag we modeled can successfully replicate the crook in the cat’s
tail. Unlike the steady-state model shown in Figure 13, ISM drag only affects
the tip of the cat’s tail because the stream is size-sorted with smaller
grains near the tip. For the model shown in Figure 16, we adopted a density of
$0.09$ g cm-3 for our fluffy refractory organic grains (5% of 1.8 g cm-3; see
Zubko et al. (2004)). Such a low grain density leads to potentially plausible
ISM parameters. For Figure 16, we adopted an ISM density of 30 H atoms cm-3
and a relative ISM velocity of 4 au/yr in the $+x$ direction and 1 au/yr in
the $+z$ direction, where $+x$ points to the SW along the secondary disk’s
midplane, and $+z$ points to the NW along the normal to the secondary disk’s
midplane. This ISM velocity is $\sim-3$ au/yr both in the directions of right
ascension and declination, broadly consistent with the proper motion of
$\beta$ Pic. We did not attempt to exhaustively explore the ISM drag parameter
space and simply note that ISM density and relative speed are somewhat
degenerate; ISM density can be decreased for a relatively small increase in
relative speed. We also note that we did not include circumstellar gas drag,
which may complicate our model.

Figure 16: Synthetic 23 $\mu$m thermal emission image of a multi-component
sporadic dust production model with an ISM drag force, which broadly
reproduces the shape, color, and flux of the cat’s tail ($\sim 150$ years
old), as well as the presence of additional arcs to the SW ($\sim 150-200$
years old) and NE ($\gtrsim$1000 years old), as marked in Figure 4. The left
panel shows the projected view oriented the same as Figure 4, and the right
panel shows a face-on view looking down onto the disk from the NW.

The same porous refractory organic grains invoked to explain the shape of the
cat’s tail, as well as the crook in the cat’s tail, also help to explain the
flux and color of the cat’s tail. To measure the flux along the cat’s tail, we
had to isolate and subtract the “background” disk flux (i.e., the main disk)
from the cat’s tail. To do so, we created a mask that followed the peak of the
cat’s tail with a half width of 3 pixels. We estimated backgrounds on a per-
pixel basis within this region by interpolation using a high-order two
dimensional polynomial fit to the regions located 6-12 pixels on either side
of the cat’s tail. After subtracting backgrounds, we measured the F1550C and
F2300C average surface brightness along the cat’s tail in steps of
$0.44\arcsec$. We find that the F1550C and F2300C surface brightness decreases
by a factor of $4-5$ from a projected separation of $8.25\arcsec$ to
$11.4\arcsec$.

We then used dustmap to synthesize F1550C and F2300C images of our cat’s tail
model shown in Figure 16. We used Mie theory and adopted optical constants
appropriate for 95% porous organic refractory material (Li & Greenberg, 1997).
We approximated the bandpass-integrated flux of each filter by integrating
over 100 individual wavelengths. We approximated each bandpass as a top-hat
function ranging from $15.1-15.9$ $\mu$m and $20.5-25.0$ $\mu$m for the F1550C
and F2300C filters, respectively. No background subtraction was needed for the
model, as the cat’s tail component was modeled in isolation.

The left panel of Figure 17 shows the F2300C flux as a function of separation
for the main disk midplane as well as the background-subtracted cat’s tail.
Empty blue circles show the cat’s tail without deconvolution, while filled
blue circles show the cat’s tail with deconvolution. The black line indicates
the model flux. While our model does not reproduce the shape of the cat’s tail
flux as a function of separation, it is on par with the total decrease in flux
from base to tip.

The right panel of Figure 17 shows the ratio of F1550C to F2300C surface
brightness along the cat’s tail for the convolved and deconvolved reduced data
sets, with the main disk’s color shown for reference. The black line indicates
the color of our sporadic dust production cat’s tail model. Porous organic
refractory grains produce a color on par with, or in excess of the
observations over the entirety of the cat’s tail. This is due to the excess
heating of small grains, as evidenced by the solid and dashed blue lines in
Figure 14. The differences between the model and data in Figure 17 suggest
that further modeling efforts are needed to better constrain composition,
porosity, and grain size distribution, and that the composition may be a
layered mixture of organic refractory material and silicates as suggested by
Li & Greenberg (1998).

Figure 17: Left: F2300C flux as a function of apparent separation measured
along the main disk’s midplane, as well as the measured and modeled cat’s
tail. Right: Color (F1550C/F2300 flux ratio) of the components. Our model with
porous organic refractory grains largely reproduces the flux and color over
the length of the cat’s tail.

To estimate the mass of the cat’s tail, we scaled our optically-thin model
such that the total flux of the cat’s tail matched observations. We adopted a
density of $0.09$ g cm-3, appropriate for highly porous refractory organic
material, and calculated the size-integrated mass. We find that our model
requires $2\times 10^{19}$ kg of dust smaller than 1 mm, or roughly the
equivalent of a large main belt asteroid in the solar system. A total mass
estimate for the original parent body is highly uncertain, as we have no
constraints on the size distribution or density of larger objects. However,
extrapolating that size distribution while maintaining the density for porous
grains, and assuming only 10% of the original object is fragmented, suggests a
total mass potentially $\sim 100$ times larger, corresponding to the mass of a
small moon in the solar system. The fainter tendrils to the SW and NE have
$\sim 5-100$ times less flux than the cat’s tail and would require
correspondingly less mass.

###  5.3 Collision with a body on a highly inclined orbit

We note that one seemingly obvious explanation for the apparent high
inclination of the cat’s tail could be simply that the progenitor was on a
highly inclined orbit. We cannot rule out this scenario. However, in addition
to invoking an unseen object on a highly inclined orbit, this scenario has
several shortcomings. We would expect a collision of this nature to travel
outward as a detached cloud. To remain attached to the main disk over $\sim
100$ au, the collision would have to produce a cloud of material with a very
broad velocity dispersion. Further, the observed CO emission appears to reside
within the secondary disk. If created by a recent collision with a highly
inclined progenitor, we would expect the CO clump to show a more significant
out-of-plane extension. Finally, it is not clear how this scenario could
explain the complex tendrilous nature of the surrounding nebulosity.

##  6 Discussion

###  6.1 Structure and morphology

Forty years of study have yielded an extensive literature about the morphology
and composition of the $\beta$ Pic system. Observations in scattered light,
both from ground (e.g. Golimowski et al., 1993; Larwood & Kalas, 2001; Milli
et al., 2014; Esposito et al., 2020) and space (e.g. Burrows et al., 1995;
Heap et al., 2000; Golimowski et al., 2006; Apai et al., 2015) have
characterized the disk. Almost all of them recover the presence of a warp,
linked to an inclined secondary disk component, and an asymmetry in the
surface brightness profiles of the NE and SW sides of the disk. Within the
inner few arcseconds the NE side of the disk is fainter, but at larger
separations the NE side is brighter and more extended (a result which can be
traced all the way back to the original pioneering observations of Smith &
Terrile, 1984, and is seen especially clearly in the deep VLT/FORS1 images of
Janson et al., 2021). The scattered light observations with NIRCam are broadly
consistent with these prior results.

There are some structures, like the SW dust and CO/C gaseous clumps at around
50-60 au, that are only visible in thermal emission, and are therefore only
recovered at longer wavelengths, i.e. mid-IR (Telesco et al., 2005) and sub-
millimeter (Dent et al., 2014; Cataldi et al., 2018). In recent work, Skaf et
al. (2023) identified several “clumps” in mid-IR images at locations closer
than 5″ to the star. They suggested the main clump, which they locate at
$\sim$52-56 au, has moved when comparing observations spanning 16 years.
However, another paper by Han et al. (2023) suggest that the change in
position is not significant considering the uncertainties. Our results (see
Fig. 6) indicate the asymmetry between the NE and SW sides of the disk has a
clumpy structure, but however the location of the clumps we detect is not
coincident with the ones reported in Skaf et al. (2023) with the exception of
their C1 clump. This corresponds to the brightest clump we observe in both
MIRI filters. The measurements of peak position described in Sect. 3.4
indicate that the location is wavelength dependant, with the clump at 23 µm
being located further from the star. This would be consistent with each
wavelength tracing material at different temperatures, althoug we should not
rule out the possibility of instrumental artifacts, given the location of the
F2300C mask. Our measurement at 15.5 µm is fairly consistent with the
locations reported in Skaf et al. (2023) and Han et al. (2023). Fig. 6
includes shaded regions with some of the estimations from Skaf et al. (2023)
and Han et al. (2023).

In the same region, millimeter observations with ALMA by Dent et al. (2014)
and later Cataldi et al. (2018) revealed the presence of CO and C I gas peaks,
respectively, as well as a significant asymmetry in the continuum, potentially
corresponding to the hypothesized dust clump. Analysis of the gaseous
component in Dent et al. (2014) revealed that the gaseous disk is inclined
with respect to the continuum counterpart by $\sim$5∘, similar to the findings
for the secondary disk (e.g. Heap et al., 2000; Golimowski et al., 2006).
Their velocity analysis also showed the rotation of the disk, with the SW side
approaching the observer, and the NE side moving away; as shown in Fig. 12.
More recently, the re-analysis of ALMA band 6 data (Matrà et al., 2019),
showed no evidence of a NE/SW asymmetry. While they claim that their result is
consistent with no evidence of large dust grains causing the asymmetry due to
the lack of detection at 24 µm in Telesco et al. (2005) and in Vandenbussche
et al. (2010), there is a possibility that the PSF FWHM in those works is too
large to characterize the dust clump. Further JWST MIRI observations could
investigate the structure of the disk at 25 µm, but the lack of a high
sensitive space based mission in the far-IR makes it impossible to contrast
results from _Herschel_.

The most sensitive observations in scattered light close to the star are
possibly those presented in Apai et al. (2015). They performed a comprehensive
study of $\beta$ Pic as seen with HST/STIS, including two epochs over 15 years
based on which they report no time variations in the brightness or morphology
of the disk. They recovered the asymmetry between the NE and SW sides in the
main plane of the disk, but they didn’t find significant evidence of the dust
clump, reinforcing the conclusion that the clump cannot be detected in
scattered light, only in thermal emission. Apai et al. (2015) suggest the
asymmetries observed might be caused by resonances with an unseen planet, or,
as we suggest here, a major collision or outgassing event in the disk.
Interestingly, when they performed a study of the vertical profile, they found
that the peak values of the secondary component left the plane of the main
disk in the SW side, towards much higher vertical distances (see Figs. 6, 7
and 9 in Apai et al., 2015). This resembles our result for the secondary disk
characterization in the MIRI data, where the cat’s tail dominates the
secondary component at distances larger than $\sim$100 au (5″). Their Figure
8, showing asymmetries after subtracting an estimate of the main disk,
highlights that on the SW side the emission is brighter and wider, but on the
NE side faint nebulosity extends outwards from the inclined secondary disk,
seen about 5°counterclockwise from the main disk plane and visible out to
$\sim 9$ ″separation. This is spatially coaligned with the inner portion of
the more extended “fork” seen in the MIRI data, and is consistent with being
the inner portion of the eccentric dust disk originating from steady-state
collisions at the collision point.

The main disk is almost but not precisely edge on, but the inclination appears
to vary between disk components at different radii. For instance, Milli et al.
(2014) reported the presence of a bow in the disk towards the NW, also
reported in the optical by Golimowski et al. (2006) and Apai et al. (2015),
and attributed to forward-scattering in a not completely edge-on disk.
Likewise, in all NIRCam filters the NW side of the inner disk ($\lesssim
3^{\prime\prime}$ or 60 au) is brighter than the SE side, consistent with the
NW side being slightly inclined toward us. On the other hand, at much wider
separations ($\gtrsim 25^{\prime\prime}$ or 500 au) Janson et al. (2021) find
the far outer wings of the disk are bowed in the opposite sense, which may
require invoking an additional wider disk component or outer warp to explain.

The understanding we have about the overall structure of the disk invokes now
several components, including some already known such as the main disk and a
secondary axisymmetric disk responsible for the warp in scattered light; and
some new, reported here for the first time, such as the fork, “cat’s tail” and
nebulosity, all apparently indicative of a collisional scenario. The asymmetry
observed between the NE and SW sides in both MIRI filters, consistent of
possibly various clumps —some of them already reported in the literature—
might also be related to the sporadic dust production model, but a more
detailed multi-wavelength analysis is probably required to further investigate
its nature, and whether its linked to a collision (Han et al., 2023) or rather
to planetary interactions(Skaf et al., 2023).

###  6.2 Composition

Li & Greenberg (1998) fit a model of the $\beta$ Pic disk to photometry and
spectra from a broad range of data sets spanning from 2 to 1300 $\mu$m.
Because these observations were unresolved, degeneracies between separation,
temperature, density, and size distribution were unavoidable in disk-
integrated fits. Nonetheless, Li & Greenberg (1998) concluded that the
particles must be very fluffy (potentially 97.5% porosity) and have an organic
refractory mantle surrounding a silicate core to reproduce silicate features
as well as the mid-IR flux, broadly consistent with our sporadic dust
production model of the cat’s tail and secondary disk presented in Sect. 5. Li
& Greenberg (1998) also found that the size distribution must have enhanced
small particles compared to comet Halley (including those that contribute to
the mid-IR spectrum). While Li & Greenberg (1998) did not adopt a simple power
law size distribution, this roughly translates into $dN/ds\propto s^{-3.1}$
(based on their Figure 2), shallower than the size distribution we invoke for
the cat’s tail. This difference may be due to the fact that we are fitting an
isolated structure, whereas Li & Greenberg (1998) fit a disk-integrated
spectrum.

Contrary to the results of Li & Greenberg (1998), Golimowski et al. (2006)
modeled resolved scattered light observations of $\beta$ Pic and found the
best fit composition to the main disk to be _non-porous_ or moderately porous
silicate and graphite grains. They ruled out high porosity on the basis that
such grains would be significantly bluer in scattered light than observations.
Ballering et al. (2016) also investigated the dust composition in $\beta$ Pic,
simultaneously fitting both spatially resolved scattered light _and_ thermal
emission observations. Similar to Golimowski et al. (2006), their best fit
suggests a composition of _non-porous_ grains comprised of 60% astronomical
silicates and 40% organic refractory material. They found no strong evidence
for water ice, consistent with our measurement in Sect. 4.2. Ballering et al.
(2016) concluded that the $\beta$ Pic halo contained sub-blowout sized grains,
consistent with the conclusions of Li & Greenberg (1998) and with recent dust
production events like we invoke for the cat’s tail.

How can fits to spatially resolved and scattered light data sets conclude no
porosity while fits to the MIR SED suggest very high porosity? As Figure 8
suggests, the $\beta$ Pic disk may simply be composed of two distinct
populations of grains—one of these may dominate scattered light while the
other dominates the mid-IR flux. As modeled in Section 5, the main disk may be
comprised predominantly of low porosity grains, while the secondary disk is
primarily composed of high porosity organic refractory material. Figure 18
shows the SEDs for porous (blue) and non-porous (orange) organic refractory
grains for various grain sizes and circumstellar distances. The 95% porous
grains scatter $\sim 300$ times less starlight than their non-porous
counterparts near 1 $\mu$m, yet can be $\sim 100$ times brighter in the mid-
IR. Such grains would be faint in scattered light observations and bright in
the mid-IR, potentially explaining the lack of a detection of the cat’s tail
in scattered light. Indeed, near-IR surface brightnesses predicted for the
cat’s tail using the SEDs for 95% porous organic grains are well below the
sensitivity limits of the NIRCam observations. Golimowski et al. (2006)
concluded that the two disk components must be composed of different grain
populations, noted that the secondary disk is tentatively bluer than the main
disk, and found that the SW extension of the main disk reddens more quickly
than the NE extension, all of which is consistent with our model of an
eccentric secondary disk composed of porous organic refractory grains with
periastron to the SW. The MIRI bandpasses may provide a unique window into
this population of dust.

Figure 18: Spectral energy distribution for identical disks with porous (blue)
and non-porous (orange) organic refractory material. Porosity reduces the
scattered light by a factor of $\sim 300$ and increases the MIR flux by a
factor of $\sim 100$. A secondary disk and cat’s tail composed of porous
organic refractory material would appear faint in scattered light and bright
in the MIR, potentially explaining the relative absence of these features in
scattered light as well as the varying best fit compositions appearing in the
literature.

###  6.3 Origin and evolution

Possible origins of the secondary disk and its asymmetries have been discussed
at length in the literature and have included a recent giant impact, steady-
state collision point models, secular forcing by a planet, resonant trapping
by a planet, and tidal disruption (e.g., Dent et al., 2014; Jackson et al.,
2014; Nesvold & Kuchner, 2015; Matrà et al., 2017a; Cataldi et al., 2018).
Here, we posit that the secondary eccentric disk is indeed generated in a
manner consistent with a long-term collision point model, such that much of
the CO and C emission may be frequently generated, but the cat’s tail and the
CO/C _clumps_ are the result of an individual recent dust production event
superimposed on the secondary disk. We contend that individual, sporadic
collisions like what created the cat’s tail are the mechanism that populates
and maintains the quasi-steady state secondary disk.

Dent et al. (2014) argued against the idea of a recent collision noting the
short CO photodissociation timescale within the clump, with literature
estimates varying from 50 to 300 years (Visser et al., 2009; Dent et al.,
2014; Matrà et al., 2017b; Cataldi et al., 2018). Such a short
photodissociation timescale implies a very recent collision, which Dent et al.
(2014) deemed less likely. Similarly, Cataldi et al. (2018) argued against a
recent collision based on the short timescale for atomic C production ($\sim
5000$ yrs) assuming production primarily via CO photodissociation. However,
Cataldi et al. (2018) pointed out that the C exhibits the same clumpy nature
as the CO, which should have spread out azimuthally within this $\sim 5000$
years.

Our dynamical model of the cat’s tail suggests a recent dust production event,
with a timescale of just $\sim 100$ yrs since production, giving credence to a
recent dust production event as the source of the continuum and gas emission
clumps. We emphasize that our quasi-steady-state episodic dust production
model does not require the entire secondary disk to be produced recently–only
the asymmetries/clumps. Our models also suggest abundant organic refractory
material in the secondary disk, which may enhance production of C through
destruction of carbon-bearing molecules (Cataldi et al., 2018), such that
production of C need not arise solely from CO photodissociation. While this
material does not necessarily imply a carbon-enriched progenitor, we find it
likely that a carbon-enriched progenitor would exhibit more of this material,
potentially consistent with the observed over-abundance of C gas (Roberge et
al., 2006). Alternatively, it is possible that the progenitor’s bulk
composition was not carbon-enriched while the ejected material is, perhaps due
to a collision stripping the outer layers from a differentiated body.

What is the probability of such a massive collisional event occurring in the
last 150 years? To estimate the probability, we distributed objects among 20k
logarithmically-spaced radius bins according to a Dohnanyi size distribution,
ranging from 10 microns (roughly consistent with the expected blowout size of
the main disk) to 1000 km. We normalized our distribution to produce a size-
integrated optical depth of $1.7\times 10^{-3}$ over a 10% wide belt centered
at 85 au, where the optical depth is taken to be equal to $L_{\rm
IR}/L{\star}$ (Kral et al., 2017). For each bin of radius $s$, we then
calculated the collision time $t_{\rm coll}\\!\left(s\right)=t_{\rm
orbit}/\left(4\pi\tau^{\prime}\\!\left(s^{\prime}\right)\right)$ (Wyatt,
1999), where $t_{\rm orbit}$ is the orbital period at 85 au,
$\tau^{\prime}\\!\left(s^{\prime}\right)$ is the optical depth of bodies
larger than than $s^{\prime}=s/2$, and $s^{\prime}$ is the smallest object
capable of fragmenting an object of radius $s$. The probability of collision
for any single body within the past 150 years can be approximated using the
Binomial theorem as $p\\!\left(s\right)=150\;{\rm yr}/t_{\rm
coll}\\!\left(s\right)$. The probability of at least one collision occurring
in the past 150 years among all bodies larger than the minimum progenitor is
approximately given by

$\sum_{s=s_{\rm
min}}^{\infty}1-\left(1-p\\!\left(s\right)\right)^{N\\!\left(s\right)},$ (1)

where $N\\!\left(s\right)$ is the number of objects in our bin of radius $s$.

In Section 5.2 we estimated the mass of the progenitor to be at least
$10^{19}$ kg and potentially as large as $10^{21}$ kg. In the solar system,
these two masses correspond roughly to radii of $\sim 100$ km and $\sim 500$
km, respectively. Adopting these radii for $s_{\rm min}$, we find a 65%
probability of at least one collision of a body larger than 100 km in radius
occurring in the past 150 years, and a 0.5% probability of at least one
collision of a body larger than 500 km in radius.

We stress that estimates of the probability of collisional events are wildly
uncertain. The estimate above must rely on extrapolations of a power law size
distribution over many orders of magnitude in size. Small changes in this
power law can greatly affect the estimated probability. Small changes in the
assumed blowout size that helps benchmark the optical depth can also greatly
affect the assumed probability. Precise estimates would require knowledge of
the absolute number and orbits of large planetesimals, for which we have no
knowledge. As a result, these estimated probabilities come with
unquantifiable, but large, uncertainties. We therefore conclude that the
calculated probabilities are consistent with our model, in that they do not
clearly rule out such a scenario.

$\beta$ Pic’s large number of comets—in fact, the largest of any exoplanetary
system studied thus far—might also tie in to our proposed dynamical scenario.
The exocomets were first detected as variable absorption features in the
spectra of $\beta$ Pic due to refractory materials (Ferlet et al., 1987;
Lagrange-Henri et al., 1989). The high frequency of exocomet events led to
several studies trying to find an origin (Beust et al., 1990, 1991), which
quickly linked the presence of a perturbing body in the disk (Beust &
Morbidelli, 2000; Thébault & Beust, 2001). The planet $\beta$ Pic b, detected
years later (Lagrange et al., 2009), seemed a natural match to the planetary
perturbation hypothesis. Through a study of hundreds of exocometary transits,
Kiefer et al. (2014) found that two families of exocomets can be distinguished
in orbit around $\beta$ Pic. They attribute the origin of one of the families,
named _S_ due to shallow transits, to the resonances caused with planet b.
However, family _D_ comets, with deeper transits and larger periastrons, have
a narrow distribution of periastron longitudes, and therefore very similar
orbits. This might indicate that they originate from the disruption of one or
several individual larger exocomets or parent bodies. A recent major dust
production event or events could plausibly send debris inward toward the star
as well as outwards. Whether family _D_ may be linked at all to a similar
scenario such as the one described in Sect. 5 will require a more in depth
analysis.

Our sporadic dust production model suggests a cat’s tail comprised of $\sim
100$ K dust with $\beta\sim 10$. Such a model leads to several testable
predictions. First, the spectra of the cat’s tail should reveal an SED similar
to the blue curves (for porous organic grains) shown in Figure 14, and the
variation of this spectrum over the length of the cat’s tail should suggest
smaller grains near the tip. Second, the high $\beta$ value suggests that dust
in the mid-to-outer regions of the cat’s tail, in spite of being located at
circumstellar distances approaching $\sim 1000$ au, should be fast moving. We
predict a $\sim 200$ mas shift in the location of the cat’s tail over three
years. Finally, our model predicts that the secondary disk should not be
perfectly edge-on, such that larger grains should track an ellipse with a
semi-minor axis $\sim 0.2$″, which should be resolvable with ALMA.

##  7 Conclusions

The JWST high contrast images of the $\beta$ Pic disk presented here have
proven to be a crucial addition to the long effort to characterize and
understand this complex system. The unprecedented sensitivity along with high
spatial resolution in near and mid-IR wavelengths allow for a better
understanding of the dust distribution along the disk, clearly revealing the
two disk components already inferred in previous works and detecting several
new components. In this work, we have reported:

  * •

The presence of two distinct disk components, with the secondary component
having a relative position angle to the main disk of 5.6$\pm$0.4 ∘. Moreover,
the main and secondary disk exhibit different colors in MIRI’s bandpasses,
suggesting the material in the secondary disk is hotter, and therefore
presumably a different composition than the main disk. We also identified a
“fork” that only appears on the NE side of the disk extending beyond 15″,
appearing in MIRI observations and attributed to this secondary disk. This
inner portion of this secondary component has been observed before in
scattered light at similar stellocentric distances, appearing as an apparent
warp in the disk. Composite color images of the MIRI observations confirm the
secondary disk is eccentric, with periastron to the SW, and is comprised of
hotter grains. NIRCam observations are consistent with previous scattered
light studies of the disk.

  * •

A curved asymmetric feature to the SW side of the disk, which we dub the
“cat’s tail.” This feature exhibits the same F1550C/F2300C color as the
secondary disk and seems to have its origin at the location of the CO clump
and mid-IR asymmetry. The cat’s tail is only detected in mid-IR, i.e. thermal
emission. This is consistent with porous organic refractory grains, which can
appear bright in the MIR and faint in scattered light. Dedicated spectroscopic
and/or imaging studies are needed at additional wavelengths in order to
further investigate the dust composition and size distribution.

  * •

Tendrilous nebulosity to the SW and W of the disk. These arc-like features
also exhibit the same F1550C/F2300C color as the secondary disk.

We presented a sporadic dust production model dynamically linking the
eccentric secondary disk to the cat’s tail and tendrilous nebulosity observed
to the SW. We posited that the eccentric secondary disk is due to ancient or
ongoing dust production, similar to the steady state collision point model of
Jackson et al. (2014), but the cat’s tail and nebulosity, as well as the
previously observed mm wavelength excess and CO clump (Dent et al., 2014), are
the outcome of individual recent dust production events occurring at locations
that deviate from the collision point. To reproduce the curvature of the cat’s
tail, we require the production of abundant sub-blowout grains produced within
the last $\sim 100-200$ years, which is on par with the CO photodissociation
timescale.

We find that porous organic refractory material could simultaneously explain
the shape of the cat’s tail, the color of the secondary disk in the MIR, the
flux along the extent of the cat’s tail in the MIR, the apparent absence of
these features in scattered light, and potentially the previously detected
clump of atomic C (Cataldi et al., 2018). This material has been previously
invoked to fit the spectral energy distribution of the $\beta$ Pic disk and
suggests a cometary origin (Li & Greenberg, 1998). JWST’s MIRI instrument may
provide the unique capability to characterize this population of dust. In-situ
investigations of small bodies in the Solar System, such as OSIRIS-REx, will
also provide new information about composition of asteroids and comets, and
serve as a baseline for research on debris disks composition.

If the material in the disk has been stirred, it could also explain the large
influx of small sized bodies observed in the system as exocomets. While
exocomets have been previously studied as originated due to resonances with
the planets, new dynamical studies are needed to better understand their
origin.

Future observations from both ground and space should aim to expand the
wavelength and time regime in which the cat’s tail is detected, in order to
characterize the material and better constrain the origin of the collisional
event(s). The proposed model makes testable predictions for the dust SED along
the cat’s tail region and for proper motion over time of the features seen
with MIRI.

The authors are grateful for the many efforts of the JWST observatory
operations teams, in particular for the flight operations team members who
responded to the spacecraft anomaly and payload safing event which occurred
during the first attempt at these NIRCam observations, diagnosed the problem,
and returned JWST to science operations in the midst of the holiday season and
shortly before the one year anniversary of its launch. The authors thank
Alycia Weinberger, Aki Roberge, Stefanie Milam, and Nathan Roth for insightful
conversations while attempting to model this disk. This paper reports work
carried out in the context of the JWST Telescope Scientist Team
(https://www.stsci.edu/~marel/jwsttelsciteam.html) (PI: M. Mountain). Funding
is provided to the team by NASA through grant 80NSSC20K0586. Based on
observations with the NASA/ESA/CSA JWST, obtained at the Space Telescope
Science Institute, which is operated by AURA, Inc., under NASA contract NAS
5-03127. The data presented in this article were obtained from the Mikulski
Archive for Space Telescopes (MAST) at the Space Telescope Science Institute.
The specific observations analyzed can be accessed via
https://www.doi.org/10.17909/8yq1-qv46 (catalog DOI). These observations are
associated with JWST program 1411 (PI Stark). I.R. is supported by grant
FJC2021-047860-I and PID2021-127289NB-I00 financed by MCIN/AEI
/10.13039/501100011033 and the European Union NextGenerationEU/PRTR.

## References

  * Apai et al. (2015) Apai, D., Schneider, G., Grady, C. A., et al. 2015, ApJ, 800, 136, doi: 10.1088/0004-637X/800/2/136
  * Arnold et al. (2019) Arnold, J. A., Weinberger, A. J., Videen, G., & Zubko, E. S. 2019, AJ, 157, 157, doi: 10.3847/1538-3881/ab095e
  * Artymowicz (1988) Artymowicz, P. 1988, ApJ, 335, L79, doi: 10.1086/185344
  * Augereau et al. (1999) Augereau, J. C., Lagrange, A. M., Mouillet, D., Papaloizou, J. C. B., & Grorod, P. A. 1999, A&A, 348, 557.  https://arxiv.org/abs/astro-ph/9906429
  * Ballering et al. (2016) Ballering, N. P., Su, K. Y. L., Rieke, G. H., & Gáspár, A. 2016, ApJ, 823, 108, doi: 10.3847/0004-637X/823/2/108
  * Beichman et al. (2010) Beichman, C. A., Krist, J., Trauger, J. T., et al. 2010, PASP, 122, 162, doi: 10.1086/651057
  * Betti et al. (2022) Betti, S. K., Follette, K., Jorquera, S., et al. 2022, AJ, 163, 145, doi: 10.3847/1538-3881/ac4d9b
  * Beust et al. (1990) Beust, H., Lagrange-Henri, A. M., Vidal-Madjar, A., & Ferlet, R. 1990, A&A, 236, 202
  * Beust & Morbidelli (2000) Beust, H., & Morbidelli, A. 2000, Icarus, 143, 170, doi: 10.1006/icar.1999.6238
  * Beust et al. (1991) Beust, H., Vidal-Madjar, A., & Ferlet, R. 1991, A&A, 247, 505
  * Boccaletti et al. (2015) Boccaletti, A., Lagage, P. O., Baudoz, P., et al. 2015, PASP, 127, 633, doi: 10.1086/682256
  * Boccaletti et al. (2022) Boccaletti, A., Cossou, C., Baudoz, P., et al. 2022, A&A, 667, A165, doi: 10.1051/0004-6361/202244578
  * Burrows et al. (1995) Burrows, C. J., Krist, J. E., Stapelfeldt, K. R., & WFPC2 Investigation Definition Team. 1995, in American Astronomical Society Meeting Abstracts, Vol. 187, American Astronomical Society Meeting Abstracts, 32.05
  * Bushouse et al. (2023) Bushouse, H., Eisenhamer, J., Dencheva, N., et al. 2023, JWST Calibration Pipeline, 1.12.1, Zenodo, doi: 10.5281/zenodo.8380331
  * Carter et al. (2023) Carter, A. L., Hinkley, S., Kammerer, J., et al. 2023, ApJ, 951, L20, doi: 10.3847/2041-8213/acd93e
  * Cataldi et al. (2018) Cataldi, G., Brandeker, A., Wu, Y., et al. 2018, ApJ, 861, 72, doi: 10.3847/1538-4357/aac5f3
  * Crotts et al. (2022) Crotts, K. A., Draper, Z. H., Matthews, B. C., et al. 2022, ApJ, 932, 23, doi: 10.3847/1538-4357/ac6c86
  * Dawson et al. (2011) Dawson, R. I., Murray-Clay, R. A., & Fabrycky, D. C. 2011, ApJ, 743, L17, doi: 10.1088/2041-8205/743/1/L17
  * Debes et al. (2009) Debes, J. H., Weinberger, A. J., & Kuchner, M. J. 2009, ApJ, 702, 318, doi: 10.1088/0004-637X/702/1/318
  * Dent et al. (2014) Dent, W. R. F., Wyatt, M. C., Roberge, A., et al. 2014, Science, 343, 1490, doi: 10.1126/science.1248726
  * Dohnanyi (1969) Dohnanyi, J. S. 1969, J. Geophys. Res., 74, 2531, doi: 10.1029/JB074i010p02531
  * Draine & Lee (1984) Draine, B. T., & Lee, H. M. 1984, ApJ, 285, 89, doi: 10.1086/162480
  * Esposito et al. (2020) Esposito, T. M., Kalas, P., Fitzgerald, M. P., et al. 2020, AJ, 160, 24, doi: 10.3847/1538-3881/ab9199
  * Ferlet et al. (1987) Ferlet, R., Hobbs, L. M., & Vidal-Madjar, A. 1987, A&A, 185, 267
  * Gardner et al. (2023) Gardner, J. P., Mather, J. C., Abbott, R., et al. 2023, PASP, 135, 068001, doi: 10.1088/1538-3873/acd1b5
  * Gáspár et al. (2023) Gáspár, A., Wolff, S. G., Rieke, G. H., et al. 2023, Nature Astronomy, 7, 790, doi: 10.1038/s41550-023-01962-6
  * Girard et al. (2022) Girard, J. H., Leisenring, J., Kammerer, J., et al. 2022, in Society of Photo-Optical Instrumentation Engineers (SPIE) Conference Series, Vol. 12180, Space Telescopes and Instrumentation 2022: Optical, Infrared, and Millimeter Wave, ed. L. E. Coyle, S. Matsuura, & M. D. Perrin, 121803Q, doi: 10.1117/12.2629636
  * Golimowski et al. (1993) Golimowski, D. A., Durrance, S. T., & Clampin, M. 1993, ApJ, 411, L41, doi: 10.1086/186907