A Warped Dust Ring in the CQ Tau Disk

Introduction

Exoplanet observations have shown that planetary systems exhibit a wide variety of orbital inclinations and misalignments [55, 1]. This diversity suggests that protoplanetary disks, the birthplaces of planets, should be understood not only as single planar structures but also as systems with three-dimensional geometries [49, 28]. High-angular-resolution observations with the Atacama Large Millimeter/submillimeter Array (ALMA) have revealed that protoplanetary disks commonly host a variety of substructures, such as rings, gaps, spiral arms, and azimuthal asymmetries (e.g., [2, 32, 61, 47]. More recently, warped disks, in which the disk rotation axis varies with radius, have been reported in several protoplanetary disks (e.g., [35, 12, 44, 40] (Bi et al. 2020). Bohn et al. [10] constrained orientations of compact inner disks in 20 transition disks using near-infrared (NIR) interferometric observations and found significant misalignments between the inner and outer disks in six systems. Such misaligned inner disks can partially block stellar irra-

Corresponding author: Shu Ishibashi

[email protected] diation and cast shadows on the outer disks, providing an observational signature of disk misalignment (e.g., Marino et al. [35]). Theoretically, disk warps and/or misalignments can be induced by gravitational perturbations from companions (Juhász & Facchini [26]), stellar flybys (Kimmig et al. [27]), misaligned infall (Kuffmeier et al. [29]), and magnetic star–disk interactions (Foucart & Lai [21]). Such variations in the rotation axis suggest that disk structure should be characterized not only in terms of localized substructures but also in terms of how the disk orientation changes with radius. However, most of these warps have been identified through gas kinematics or NIR observations. The geometry of the dust disk should also be important, because dust grains constitute the solid reservoir from which planetesimals and planetary cores form. Establishing an observational framework to quantify dust-disk warps as radial variations in disk orientation is therefore essential for connecting disk geometry with planet formation.

CQ Tau, a relatively evolved protoplanetary disk system with an estimated age of $\sim10$ Myr (Garufi et al. [23]; Wallack et al. [54]) at a distance of 149.4 pc (Gaia Collaboration et al. [22]), is a particularly promising target for investigating dust disk warps. This transition disk exhibits a variety of structures, including a bright dust ring, crescents, and arcs in the dust continuum, NIR shadows on the dust ring, as well as spiral arms seen in NIR and gas observations (Ubeira Gabellini et al. [51]; Wölfer et al. [57]; Hammond et al. [20]; Ren et al. [42]; Yamaguchi et al. [61]; Curone et al. [13]; Zagaria et al. [63]). In addition, a misaligned inner disk (with a tilt of $\sim40^\circ$) has been reported by Bohn et al. [10]. Winter et al. [56] showed that the non-axisymmetric CO line-of-sight velocity pattern of CQ Tau can be reproduced by the velocity field expected for a warped disk. The misaligned inner disk and the candidate gas warp suggest that CQ Tau may host a warp extending over a wide radial range. Gaia astrometry has revealed a proper motion anomaly that is difficult to explain with a single-star model, suggesting the presence of a companion (Blakely et al. [9]; Vioque et al. [52]).

We combine three archival ALMA Band 6 datasets of CQ Tau and reconstruct a 1.3 mm dust continuum image using the sparse-modeling-based two-dimensional super-resolution imaging method PRIISM (Nakazato & Ikeda [38]; Nakazato et al. [39]). PRIISM performs visibility fitting under the assumptions of image sparsity and smoothness, and reconstructs images consistent with the observed visibilities. In contrast to the standard CLEAN algorithm, PRIISM does not require deconvolution and convolution procedures, providing the potential to achieve higher image fidelity and super resolution. Using this image, we characterize the geometry of the dust disk as a function of radius and test whether this geometry can be explained by a warp.

This paper is organized as follows. In § 2, we describe the ALMA Band 6 data, data reduction, and CLEAN and PRIISM imaging procedures. In § 3, we present the 1.3 mm dust continuum morphology revealed by the CLEAN, PRIISM, and real-part PRIISM images. In § 4, we derive the global disk orientation and its radial variations through ellipse fitting and compare the results with those of coplanar-ring and warped-ring toy models. In § 5, we discuss the evidence for the dust warp and its possible physical origins. § 6 summarizes the main results and interpretations.

Data Reduction and Imaging

Data Reduction

We combined and analyzed three ALMA archival datasets of the 1.3 mm dust continuum emission of CQ Tau, which consist of Data 1 (2013.1.00498.S), Data 2 (2016.A.00026.S), and Data 3 (2017.1.01404.S). Data 3 consists of two execution blocks, hereafter referred to as Data 3-1 and Data 3-2. These data were observed in 2015 August, 2017 August, and 2017 November, respectively. The 1.3 mm dust continuum data from these observations have already been published (Ubeira Gabellini et al. [51]; Wölfer et al. [57]), and the configuration, baseline, and other observational properties of the datasets are summarized in Wölfer et al. [57]. Each dataset was calibrated with the standard ALMA calibration pipeline, with Common Astronomy Software Applications (CASA; McMullin et al. [36]) versions 4.3.1, 4.7.2, and 5.1.1-5 for Data 1, Data 2, and Data 3, respectively. The further calibration described below was performed with CASA version 6.1.0.

To combine Data 1, 2, 3-1, and 3-2, we first aligned their phase centers using the CASA tasks fixvis and fixplanets. In this step, we made PRIISM continuum images separately for Data 1, Data 2, and Data 3, following the procedure described in § 2.2. Here, Data 3-1 and Data 3-2 were imaged together because their observing dates were separated by only three days. These images were used to check whether the compact inner disk emission was detected and could serve as the reference position for phase alignment. For Data 1, since the compact inner emission was not detected, likely due to limited angular resolution, we adopted the ICRS coordinates R.A. = $05^\mathrm{h}35^\mathrm{m}58.4700^\mathrm{s}$ and Dec. = $+24^\circ44^\prime53.683^{\prime\prime}$, corresponding to the average of the two stellar positions expected for 2015 and 2016 from Gaia DR3 astrometry. Since Data 1 was obtained in 2015 August, the average stellar position obtained from the 2015 and 2016 Gaia observations, which corresponds to the position in the middle of 2015, should provide a sufficiently accurate stellar position. Assuming the linear proper motion, the offset between this averaged position and the stellar position at the actual observing epoch is at most $\sim5$ mas, much smaller than the effective spatial resolution of the PRIISM image for Data 1, $\sim0.12^{\prime\prime}$. For Data 2 and Data 3, the compact inner disk emission was detected in the PRIISM images, and its peak position was adopted as the reference position. The adopted coordinates were R.A. = $05^\mathrm{h}35^\mathrm{m}58.4697^\mathrm{s}$ and Dec. = $+24^\circ44^\prime53.615^{\prime\prime}$ for Data 2, and R.A. = $05^\mathrm{h}35^\mathrm{m}58.4710^\mathrm{s}$ and Dec. = $+24^\circ44^\prime53.618^{\prime\prime}$ for Data 3. The peak-based positions are close to the stellar position expected for 2018 from Gaia DR3 astrometry, with offsets of $0.011^{\prime\prime}$ for Data 2 and $0.0067^{\prime\prime}$ for Data 3. The position adopted for Data 3 was applied to both Data 3-1 and Data 3-2. Because fixplanets accepts directions only in the J2000 frame, the adopted ICRS positions were converted to J2000 before applying fixvis. The J2000 coordinate of the Data 3 position was then assigned as the common phase-center coordinate for all datasets using fixplanets. Accordingly, we use this common phase center as a proxy for the stellar position throughout this study.

After aligning the phase centers, we adjusted their flux scales. Specifically, we first deprojected the $uv$ coordinates assuming a position angle (PA) of $235^\circ$ and an inclination $i$ of $35^\circ$, following the disk geometry adopted by Ubeira Gabellini et al. [51], and compared the azimuthally averaged real visibility profiles within overlapping deprojected $uv$-ranges to derive the relative flux-scaling factors. The real- and imaginary-part profiles show that the imaginary part remains close to zero, and the null positions in the real part associated with the ring structure are consistent among the datasets. Here, Data 2 was used as the reference for the flux scaling because it provides the most extensive coverage at short $uv$-spacings. The flux discrepancies relative to Data 2 were 35%, 3%, and 8% for Data 1, Data 3-1, and Data 3-2, respectively. The corresponding antenna-based gain factors of 1.164, 1.014, and 1.040 were input to the CASA task gencal to generate the amplitude calibration tables, which were then applied to each dataset using applycal. The larger flux deviation of Data 1 may reflect an absolute flux-calibration offset in the oldest observations at Cycle 2.

Following the phase-center alignment and flux rescaling, all self-calibration steps described below were performed using the automated self-calibration pipeline AJISAI v0.1.0 ([58])¹. A detailed description of AJISAI is provided in Appendix § 10 of Yamaguchi et al. [62]. Self-calibration was first performed independently for each dataset. Two rounds of phase self-calibration, with a longer solution interval for the first round and a shorter one for the second (Data 1: inf and 121 s; Data 2: inf and 40 s; Data 3-1: inf and 40 s; Data 3-2: inf and 40 s), were applied to improve the signal-to-noise ratio (SNR). The SNR was evaluated as the ratio of the peak intensity to the rms noise measured in an emission-free circular annulus between $1.2^{\prime\prime}$ and $2.0^{\prime\prime}$ from the phase center. These self-calibration steps improved the SNRs of the individual datasets by factors of 2.5, 2.0, 2.0, and 1.8, respectively.

Finally, the combined data were self-calibrated with one round of phase self-calibration with a solution interval of inf. This final self-calibration improved the SNR of the combined dataset by a factor of 1.2 relative to its value before the final self-calibration.

Imaging with CLEAN and PRIISM

After obtaining the combined visibility data, we reconstructed 1.3 mm dust continuum images using both CLEAN and PRIISM. We first produced a CLEAN image with Briggs weighting and a default robust parameter of 0.5. The resulting CLEAN image has a synthesized beam of $0.074^{\prime\prime} \times0.058^{\prime\prime}$ with a PA of $0.6^\circ$. The rms noise was measured in an emission-free circular annulus between $1.2^{\prime\prime}$ and $2.0^{\prime\prime}$ from the phase center, yielding 2.78 mJy arcsec$^{-2}$ (13.6 $\mu$Jy beam$^{-1}$). The peak intensity is 480 mJy arcsec$^{-2}$, corresponding to a peak SNR of 173. The flux density measured from pixels above the $5\sigma$ threshold is 137.9 mJy.

To enable a direct comparison with the PRIISM image, where the effective spatial resolution is $0.052^{\prime\prime} \times0.042^{\prime\prime}$ as described below, we also produced a CLEAN image with a robust parameter of $-0.2$. The resulting CLEAN image has a synthesized beam size of $0.052^{\prime\prime} \times0.036^{\prime\prime}$ with a PA of $167.4^\circ$. The rms noise is 7.71 mJy arcsec$^{-2}$ (16.5 $\mu$Jy beam$^{-1}$), and the peak intensity is 515 mJy arcsec$^{-2}$, corresponding to a peak SNR of 67. The flux density measured from pixels above the $5\sigma$ threshold is 131.5 mJy.

For the PRIISM reconstruction, we used the Python package PRIISM version 0.12.0 ([38, 39]), which enables super-resolution imaging with improved image fidelity ([60, 59, 61, 46 b]). The hyperparameters in PRIISM, $\Lambda_l$ and $\Lambda_{\mathrm{tsv}}$, were selected us-

ing 10-fold cross-validation (CV) available in PRIISM. This procedure allows us to objectively identify the image with the highest fidelity. Following the CV procedure described in Yamaguchi et al. [60], the optimal PRIISM image was reconstructed using the hyper-parameter set $(\Lambda_l, \Lambda_{\mathrm{tsv}}) = (3 \times10^{5}, 4 \times10^{13})$. ² The hyper-parameters for the PRIISM images used for the phase-center alignment in § 2.1 were also determined using the same CV-based parameter-search procedure. The adopted hyper-parameter sets $(\Lambda_l, \Lambda_{\mathrm{tsv}})$ were $(10^{4}, 10^{9})$, $(10^{5}, 10^{12})$, and $(2 \times10^{5}, 7 \times10^{13})$ for Data 1, Data 2, and Data 3, respectively.

For the PRIISM image, the detection threshold is defined in a different way from the rms noise commonly used for CLEAN images. This is because, unlike CLEAN images, PRIISM images are reconstructed through regularized optimization under a non-negativity constraint, and the resulting pixel values are spatially correlated. Therefore, following previous PRIISM studies [59, 61], we define the detection threshold, $I_{\mathrm{DT}}$, as the maximum intensity in an emission-free region. This provides a conservative threshold for identifying emission above the reconstruction noise floor. The resulting $I_{\mathrm{DT}}$ for the optimal PRIISM image was 24.0 mJy arcsec$^{-2}$. The flux density measured from pixels above $I_{\mathrm{DT}}$ within the source region is 136.2 mJy, comparable to the 137.9 mJy and 131.5 mJy measured from the CLEAN images with robust parameters of 0.5 and −0.2, respectively.

Spatial resolution is not automatically defined in PRIISM images as it is in CLEAN, because PRIISM reconstructs images by minimizing a cost function without any beam-convolution process. We evaluated the effective spatial resolution of the PRIISM image, $\theta_{\mathrm{eff}}$, following the method of Yamaguchi et al. [59]. The method consists of three steps. First, we created four point-source-injected visibility datasets by adding an artificial point source to the observed visibility data at a distance of 1.2$^{\prime\prime}$ from the phase center toward the north, south, east, and west. The flux density of the injected point source was 7.0 mJy, corresponding to approximately 5% of the total flux of the PRIISM image. Second, we performed PRIISM imaging for each point-source-injected dataset using the optimal hyper-parameters, $(\Lambda_l, \Lambda_{\mathrm{tsv}}) = (3 \times10^{5}, 4 \times10^{13})$. Third, we fitted the injected point source in each resulting PRIISM image with a two-dimensional Gaussian function and defined $\theta_{\mathrm{eff}}$ from the mean FWHM of the fitted Gaussians. The four injection tests yielded consistent FWHM values, confirming that the estimated effective resolution is insensitive to the adopted injection positions. The resulting value was $\theta_{\mathrm{eff}} = 0.052^{\prime\prime} \times0.042^{\prime\prime}$ (or 7.8 × 6.3 au) with a PA of 177.6°.

Results

The CLEAN and PRIISM images

The CLEAN and PRIISM images are shown in Figure 1. These images show a central cavity and a bright ring at radii of ∼ 20 au to ∼ 70 au, hereafter referred to as the main ring. The main ring is not azimuthally uniform and shows two bright crescents located at PAs of ∼ 70° and ∼ 190°. Outside the main ring, a faint disk-skirt structure extending to a radius of ∼ 80 au is seen in all images. In addition, a compact emission component around the central star, defined as the inner disk, is seen in all the images, although in the CLEAN image with a robust parameter of −0.2 its peak intensity is below the 5$\sigma$ level (∼ 4.5$\sigma$). These structures are consistent with those reported in previous ALMA studies [51, 61, 13, 63].

Image not extracted

The compact inner disk is detected above the $I_{\mathrm{DT}}$ level in the PRIISM image. A two-dimensional Gaussian fit yields a size of $0.057^{\prime\prime} \times0.047^{\prime\prime}$ (PA = 172°). As this size is comparable to the PRIISM spatial resolution, the central disk remains unresolved with an upper limit of 8 au on its size. The flux density measured from the pixels above $I_{\mathrm{DT}}$ is 0.03 mJy. Assuming that the emission is optically thin, the dust mass of the inner disk is estimated as $M_{\mathrm{dust}} = F_{\nu}d^{2}/\{\kappa_{\nu}B_{\nu}(T_{\mathrm{dust}})\}$, where $\nu$ is the frequency, $F_{\nu}$ is the flux density, $d$ is the distance, $\kappa_{\nu}$ is the dust opacity, and $B_{\nu}(T_{\mathrm{dust}})$ is the Planck function at the dust temperature $T_{\mathrm{dust}}$. At a representative frequency of $\nu= 225$ GHz, we adopt $\kappa_{\nu} = 2.25$ cm$^{2}$ g$^{-1}$ [7]. Assuming $T_{\mathrm{dust}} = 50$–200 K, the dust mass is estimated to be ∼ $1.7 \times10^{-3}$–$7.3 \times10^{-3}$ $M_{\odot}$.

Even though the spatial resolution of the CLEAN image with the robust parameter of −0.2 is comparable to that of the PRIISM image, the image contrasts in the PRIISM image are sharper. The PRIISM image exhibits the crescent features at PA ∼ 70° and ∼ 190° more clearly. The radial intensity variations in the PRIISM image are also sharper than those in the CLEAN image, particularly evident at PA ∼ 120° and ∼ 300°. These sharp features may be related to the three arcs reported in the high-pass-filtered CLEAN image of Zagaria et al. [63], suggesting that PRIISM recovers these substructures more directly in the original image.

To assess the imaging fidelity of these different reconstructions, the PRIISM and CLEAN component images are Fourier-transformed and sampled onto the observed $uv$ coverage. The reduced chi-square values between these reconstructed visibilities and the original observed visibilities were evaluated in bins of 30 k$\lambda$ in deprojected $uv$ distance, adopting PA = 230.8$^\circ$ and $i = 32.9^\circ$, as determined in § 4.1.1. The reduced chi-square value of the PRIISM visibilities ($\sim1.2$) is indeed lower than that of the CLEAN components ($\sim4.3$ and $\sim5.3$ for robust parameters of 0.5 and $-0.2$, respectively), indicating higher fidelity for the PRIISM reconstruction.

In the PRIISM image, the brightness temperature $T_{\mathrm{B}}$ ranges from 3.8 K to 18.7 K over the radial extent of the main ring from 20 au to 70 au. On the other hand, the inferred $T_{\mathrm{dust}}$ of the main ring passively heated by the central source with $L_* = 9L_\odot$ ([54]) is estimated to range from 25.8 K to 48.3 K, using Equation (16) of Yamaguchi et al. [61] with the flaring angle $\phi= 0.02$. A pixel-by-pixel comparison of the corresponding maps gives $T_{\mathrm{dust}}/T_{\mathrm{B}}$ ranging from 1.7 to 12.7. The 1.3 mm dust emission is thus most likely optically thin over most of the radial extent of the main ring.

Real-part PRIISM Image

Because the real part of the visibilities traces the centro-symmetric component with respect to the phase center, we also reconstructed a PRIISM image using only the real part, allowing us to examine the underlying disk geometry centered on the stellar position while reducing the influence of non-centro-symmetric emission. For consistency with the PRIISM image reconstructed from the full visibilities, we used the same hyper-parameters, $(\Lambda_l,\Lambda_{\mathrm{tsv}}) = (3 \times10^5, 4 \times10^{13})$, and adopted the same effective spatial resolution and detection threshold. The peak intensity is $\sim460$ mJy arcsec$^{-2}$, corresponding to $\sim86\%$ of that of the full PRIISM image. The bulk of the emission is thus contained in the real part.

Figure 2(a) shows the real-part PRIISM image. Similar to the full PRIISM image, the real-part image shows the compact inner disk, central cavity, main ring, and faint extended disk-skirt structure. Near the inner edge of the main ring, a curved non-axisymmetric pattern of slightly reduced intensity is seen around PA $\sim130^\circ$. The same feature appears at PA $\sim310^\circ$ ($= 130^\circ+ 180^\circ$) by construction of real-part imaging, which produces centro-symmetric image. This pattern is further examined using high-pass filtering in § 4.2.3. The non-centro-symmetric component is highlighted in Figure 2(b), which was obtained by subtracting the real-part PRIISM image from the PRIISM image reconstructed from the full visibilities. Prominent positive residuals are seen to the south and east, spatially corresponding to the crescents in the full PRIISM image (Figure 1). Additional residuals are distributed near the inner and outer edges of the main ring.

Image not extracted

Analysis

Disk Geometry from Contour Ellipse Fitting

Global Disk Geometry

Throughout this paper, PA is measured east of north on the sky and ranges from $0^\circ$ to $360^\circ$. The inclination $i$ is defined as the angle between the disk plane and the plane of the sky, ranging from $0^\circ$ to $90^\circ$, with $i = 0^\circ$ corresponding to a face-on disk. We denote the global major-axis PA and $i$ of the CQ Tau dust disk by $\mathrm{PA}_0$ and $i_0$, respectively. $\mathrm{PA}_0$ is defined such that the far side of the disk lies along the minor axis at $\mathrm{PA}_0 + 90^\circ$.

Wölfer et al. [57] interpreted the slightly brighter CO peak intensity on the northwestern side as evidence that this side is the far side. They attributed this brightness asymmetry to a larger projected area of the far-side CO-emitting surface and less severe beam dilution, arising from the vertical structure of the emitting layer. This orientation is further supported by the consistency between the observed spiral arm morphology and the disk rotation direction ([57, 20]), following the geometric interpretation presented in Muto et al. [37]. We therefore adopt the northwestern side of the disk as the far side.

$\mathrm{PA}_0$ and $i_0$ of the CQ Tau dust disk were determined by performing ellipse fitting to a contour of the real-part PRIISM image. The ellipse center was fixed at the phase center. This choice reduces possible biases from non-centro-symmetric emission components. The free parameters are the lengths of the semimajor and semiminor axes ($a$, $b$) and the PA of the ellipse major axis. We assume that the continuum intensity approximately traces the dust density distribution since the 1.3 mm continuum emission is optically thin (§ 3.1). The inclination is estimated as $i = \arccos(b/a)$. The detailed ellipse-fitting procedure, including the calculation of the fitting uncertainties and the reduced chi-square, $\chi^2_{\mathrm{red}}$, is described in Appendix § 9.

An ellipse fit to the $3\,I_{\mathrm{DT}}$ contour tracing the outer edge of the main ring provides a reliable fit with $\chi^2_{\mathrm{red}} \simeq0.3$, yielding $\mathrm{PA} = 230.8^\circ\pm0.6^\circ$ and $i = 32.9^\circ\pm0.7^\circ$. These values are different from previous CO-based gas-disk orientations of $\mathrm{PA} \simeq235^\circ$ and $i \simeq35^\circ$–$36^\circ$ ([51, 25]). The CO-based values represent the global gas-disk orientation traced by line emission and kinematics, with the CO disk extending to a radius of $\sim150$ au (Ubeira Gabellini et al. [51]), approximately twice the radial extent of the main dust ring. In contrast, our values are derived from the outer edge of the dust ring. Since this paper focuses on the dust disk structure, we adopt the continuum-based values, $\mathrm{PA}_0 = 230.8^\circ$ and $i_0 = 32.9^\circ$, throughout this paper.

Radial Variations in Disk Geometry

The same ellipse-fitting analysis is applied over the ranges of contour levels tracing the inner and outer edges of the main ring, from $1\,I_{\mathrm{DT}}$ to $7\,I_{\mathrm{DT}}$ and from $3\,I_{\mathrm{DT}}$ to $13\,I_{\mathrm{DT}}$, respectively, in steps of $I_{\mathrm{DT}}$. The fitting results are shown in Figure 3, which summarizes the resulting PA and $i$ as functions of $a$. Here, $a$ is used as a proxy for the disk radius at each contour level. At both the inner and outer edges, the PA and $i$ vary with $a$ beyond the fitting uncertainties. At the inner edge, over $a \sim0.14''$–$0.17''$ (21–26 au), the PA and $i$ span approximately $214^\circ$–$233^\circ$ and $23^\circ$–$39^\circ$, respectively. At the outer edge, over $a \sim0.39''$–$0.45''$ (58–67 au), they span approximately $231^\circ$–$235^\circ$ and $33^\circ$–$36^\circ$, respectively. The gray shaded region in Figure 3 approximately corresponds to the radial range of the emission ridge of the main ring, where the contours are affected by azimuthal intensity variations and do not simply trace the underlying disk geometry. This region is thus excluded from the ellipse-fitting analysis, which aims to reveal radial variations in the disk orientation. At the inner edge of the main ring, the $8\,I_{\mathrm{DT}}$ contour is locally distorted by the curved reduced-intensity pattern shown in Figure 2(a). At the outer edge, the $14\,I_{\mathrm{DT}}$ contour traces only the brighter azimuthal portions of the ring rather than the overall outer edge.

Image not extracted

To examine whether the radial variations in PA and $i$ depend on the imaging method, the same ellipse fitting procedure is applied to the CLEAN image reconstructed from the real part of the visibilities with the robust parameter of $-0.2$ at the same intensity levels as for the real-part PRIISM image. The results from the $1\,I_{\mathrm{DT}}$ contour are excluded because its intensity is below the $5\sigma$ level in the CLEAN image ($\sim3\sigma$), and that of the $7\,I_{\mathrm{DT}}$ contour tracing the inner edge is discarded because of the fitting instability. The remaining data points show radial trends largely consistent with those obtained from the real-part PRIISM image at both the inner and outer edges (Figure 3), indicating that the

radial variations are not specific to the PRIISM reconstruction.

To assess a possible effect of image smearing due to the combination of multi-epoch data, the ellipse fitting is also applied to the real-part PRIISM image reconstructed from Data 3 alone, adopting $(\Lambda_l, \Lambda_{\mathrm{tsv}}) = (2 \times10^{5}, 7 \times10^{13})$ (§ 2.2). The fitting is conducted at the same intensity levels as that of the combined-data real-part PRIISM image, but the $7 I_{\mathrm{DT}}$ contour tracing the inner edge is excluded because the fitting is unstable. The overall trends are again broadly consistent with those obtained from the combined data (Figure 3), indicating that the radial variations in PA and $i$ are not primarily driven by multi-epoch smearing.

We further test a simple eccentric-orbit interpretation by applying center-free ellipse fitting to the full PRIISM image at the same contour levels as in the center-fixed analysis. In an eccentric orbit, the central star is expected to lie at one focus of the ellipse [16, 30]. However, the fitted ellipse centers are typically located within one third of the effective spatial resolution from the phase center, whereas, for every contour, the distance from the phase center to the nearer focus is more than three times the distance to the ellipse center. The observed geometry therefore disfavors a simple eccentric-orbit interpretation.

Toy Modeling

Coplanar Ring Model

The base disk structure is derived from the fitting to the real-part PRIISM image of the main ring with a coplanar ring model of multi-Gaussian functions. The model disk is assumed to be infinitesimally thin. Appendix B describes the detailed construction of the coplanar ring model. The model image was convolved to the same resolution as the PRIISM image and analyzed using the same center-fixed ellipse-fitting procedure and contour levels as for the observations.

Figure 4 and Figure 5 show the coplanar ring model images and their ellipse-fitting results. At both the inner and outer edges, the PA and $i$ profiles exhibit only modest radial variations and remain close to $\mathrm{PA}_0$ and $i_0$, failing to reproduce the observed variations, particularly around a \simeq0.17\arcsec. To assess the modest variations in PA and $i$, we repeated the ellipse fitting on a coplanar model convolved with a circular Gaussian whose FWHM was set to the geometric-mean resolution. The variations were largely reduced but did not disappear completely. Thus, the variations likely reflect not only the elliptical resolution but also uncertainties associated with extracting and fitting the pixelized contours.

Image not extracted

Image not extracted

Warped Ring Model

A warped ring model was constructed by applying a radial profile of the disk-normal direction to the coplanar ring model. The disk-normal direction at each radius (i.e., warp) is described by the polar angle $\beta(r)$ and azimuthal angle $\gamma(r)$, as illustrated in Figure 6. $\beta$ is defined over $0^\circ\le\beta\le90^\circ$. $\gamma$ is taken to lie in the range $-180^\circ< \gamma\le180^\circ$, with $\gamma= 0^\circ$ corresponding to the direction along the global disk major axis at $\mathrm{PA}_0$.

Image not extracted

The warp geometry is represented by a set of thin annuli. For each annulus, the emission is projected onto the sky plane according to the local values of $\beta(r)$ and $\gamma(r)$ while conserving its flux, as appropriate for the observed optically thin continuum emission. As a result, when an annulus is projected onto a smaller area on the sky plane, its intensity can be locally enhanced.

We modeled the radial profiles of $\beta$ and $\gamma$ as continuous piecewise-linear functions with a common break

radius $r_{\mathrm{ref}}$:

If the piecewise-linear expression for $\beta(r)$ becomes negative, we set $\beta(r)=0^\circ$. $r_{\mathrm{ref}}$ was fixed at $0.3^{\prime\prime}$, approximately the radial midpoint of the main ring. $\beta_0$ and $\gamma_0$ are the values at $r=0$, while $s_{\beta,\mathrm{inner}}$ and $s_{\beta,\mathrm{outer}}$ are the radial slopes of $\beta(r)$ inside and outside the break radius, respectively; $s_{\gamma,\mathrm{inner}}$ and $s_{\gamma,\mathrm{outer}}$ are defined analogously for $\gamma(r)$.

We manually explored the parameters defining the $\beta(r)$ and $\gamma(r)$ profiles. Each trial model was convolved to the same resolution as the PRIISM image and analyzed using the same center-fixed ellipse-fitting procedure and contour levels as for the observations. The fiducial model was selected so that its ellipse-fitting results broadly reproduce the observed radial variations in PA and $i$. The explored parameter ranges and fiducial values are listed in Table 1, and the resulting fiducial $\beta(r)$ and $\gamma(r)$ are shown in Figure 7. In the fiducial model, $\beta$ decreases outward and reaches $0^\circ$ at $r\simeq0.44^{\prime\prime}$, while $\gamma$ varies with radius. The resulting warp geometry is shown in Figure 8.

Parameter	Explored range	Fiducial value
$\beta_0$ [deg]	1 – 30	5
$s_{\beta,\mathrm{inner}}$ [deg arcsec$^{-1}$]	−90 – −2	−5
$s_{\beta,\mathrm{outer}}$ [deg arcsec$^{-1}$]	−40 – −5	−25
$\gamma_0$ [deg]	−180 – 180	90
$s_{\gamma,\mathrm{inner}}$ [deg arcsec$^{-1}$]	−2400 – −2000	−2200
$s_{\gamma,\mathrm{outer}}$ [deg arcsec$^{-1}$]	−1000 – −200	−350

Image not extracted

Image not extracted

Figure 5 and Figure 9 show the ellipse-fitting results and the warped ring model images, respectively. At the inner edge, the model PA decreases by $\sim20^\circ$, from $\sim232^\circ$ to $\sim213^\circ$. The model $i$ decreases by $\sim10^\circ$, from $\sim38^\circ$ to $\sim28^\circ$. At the outer edge, the model PA decreases by $\sim5^\circ$, from $\sim236^\circ$ to $\sim230^\circ$, while the model $i$ decreases by $\sim3^\circ$, from $\sim36^\circ$ to $\sim33^\circ$. The model (red markers in Figure 5) broadly reproduces the observed trends and amplitudes of the variations (black markers). At the inner edge, the differences between the model and observed PA and $i$ are comparable to the typical ellipse-fitting uncertainties ($\sim4^\circ$), indicating that the model reasonably reproduces the observed variations. At the outer edge, a systematic discrepancy of $\lesssim3^\circ$ remains between the model and observed PA profiles, possibly due to the simplified parametric warp profile or unmodeled intensity structures. By contrast, the coplanar model does not reproduce the observed variations (blue markers), suggesting that the warped ring model is a better representation of the observed image.

Image not extracted

Comparison of High-pass-filtered Images

Figure 10 compares the high-pass-filtered real-part PRIISM image with those of the coplanar and warped ring models. Each image was obtained by subtracting its Gaussian-smoothed version with an FWHM of $0.1^{\prime\prime}$ from the original image. The high-pass-filtered real-part PRIISM image shows multiple curved residual patterns distributed in a centro-symmetric manner. In particular, curved negative residuals along the inner edge and narrow positive residuals near the outermost radius are seen, delineated by the white and black dashed lines, respectively. The former correspond to the curved pattern of reduced intensity seen at PAs of $\sim130^\circ$ and $\sim310^\circ$ in the original real-part PRIISM image (Figure 2(a)). By contrast, the coplanar ring model shows only axisymmetric residuals associated with the multi-Gaussian ring components. In the warped model, annuli at different radii have different sky-plane projections, producing radial and azimuthal intensity variations. Such behavior appears consistent with that in the high-pass-filtered PRIISM image, while the detailed emission distributions do not perfectly match. The high-pass-filtered image of the warped model indeed reproduces both the curved negative residuals along the inner edge and the narrow positive residuals near the outermost radius (Figure 10), similar to those in the real-part PRIISM image. These qualitative trends persist for the different smoothing FWHMs examined in Appendix § 11.

Image not extracted

Discussion

Evidence for the Dust Warp

Our PRIISM imaging of the real-part visibilities enables us to investigate the detailed, systematic disk structures while reducing the influence of non-centro-symmetric components. While the coplanar ring model fails to reproduce the observed radial variations in both PA and $i$, the fiducial warped ring model, in which $\beta$ decreases from $\sim4^\circ$ to $1^\circ$ for $r = 20$ to 60 au, provides a better match to these variations. Although this toy model is based on simple parametric warp profiles and is not intended to provide a unique quantitative fit, the comparison suggests that the main dust ring in CQ Tau is warped. Further support for the dust warp comes from the center-free ellipse-fitting results, which disfavor a simple eccentric-orbit interpretation (§ 4.1.2), and from the agreement between the observed high-pass-filtered morphology and that of the warped ring model (§ 4.2.3). The latter comparison also suggests that high-pass filtering may provide a supplementary tool for identifying disk warps.

A warped disk geometry is intrinsically three-dimensional, because the disk height relative to the global disk midplane varies with radius and azimuth. In addition to the ring-scale warp that gradually weakens toward the outer edge, CQ Tau hosts a strongly misaligned inner disk with a tilt of $\beta\sim40^\circ$ ([10]). Furthermore, Winter et al. [56] showed that the non-axisymmetric CO line-of-sight velocity pattern over radii of $\sim45$–$155$ au can be reproduced by the velocity field expected for a warped disk with a tilt of a few degrees. The CQ Tau disk may therefore be interpreted as a hierarchical three-dimensional disk structure.

Physical Origin of the Warp

Companion

A companion on an inclined orbit exerts a secular gravitational torque on the disk and induces nodal precession of the disk angular momentum vector ([26, 17]). Since the precession rate generally depends on radius, different annuli precess at different rates if radial communication within the disk is inefficient, potentially leading to disk breaking. In the bending-wave regime ($\alpha< h/r$), where $\alpha$ is the viscosity parameter in the $\alpha$ prescription ([45]) and $h/r$ is the disk aspect ratio, such differential precession can be communicated across the disk by bending waves. The radial communication timescale, $t_{\mathrm{c}}$, and the global nodal precession timescale, $t_{\mathrm{n,global}}$, therefore characterize how the disk responds to the companion torque (e.g., [33, 48] (Bi et al. 2020)). When $t_{\mathrm{c}}/t_{\mathrm{n,global}} > 1$, differential precession proceeds faster than radial communication, and the disk can break into distinct regions. In contrast, when $t_{\mathrm{c}}/t_{\mathrm{n,global}} < 1$, radial communication is faster than global precession, so the warped disk responds as a single coherent unit.

Following Bi et al. (2020), we estimate $t_{\mathrm{c}}$ and $t_{\mathrm{n,global}}$ for the main ring of CQ Tau, assuming that the ring is in the bending-wave regime. For the disk aspect ratio, we adopt $h/r = 0.125(r/56\,\mathrm{au})^{0.05}$, based on the radiative-transfer modeling of Ubeira Gabellini et al. [51]. For a companion on a circular orbit, these timescales are given by

where $r_{\mathrm{in}}$ and $r_{\mathrm{out}}$ define the radial range of the main ring over which the global precession is evaluated. The quantity $(h/r)_{\mathrm{out}}$ is the disk aspect ratio at $r_{\mathrm{out}}$. The angle brackets denote the angular-momentum-weighted average over this radial range, where $\Sigma(r)$ and $\Omega(r)$ are the disk surface density and angular frequency at radius $r$, respectively. We assume $\Sigma(r) \propto r^{-3/2}$ in evaluating this average, following Bi et al. (2020). $M_{*}$ is the stellar mass, and $M_{\mathrm{comp}}$ and $a_{\mathrm{comp}}$ are the companion’s mass and orbital radius, respectively. $\Omega_{\mathrm{comp}}$ is the orbital angular frequency of the companion. Here, we adopt $r_{\mathrm{in}} = 20\,\mathrm{au}$, $r_{\mathrm{out}} = 70\,\mathrm{au}$, $M_{*} = 1.6\,M_{\odot}$ (Wallack et al. [54]), and $(h/r)_{\mathrm{out}} = 0.126$ (Ubeira Gabellini et al. [51]). Because the companion has not been directly detected, we adopt the companion solution inferred from the Gaia DR2–DR3 proper-motion-anomaly analysis of Vioque et al. [52], $M_{\mathrm{comp}} = 309\,M_{\mathrm{J}}$ and $a_{\mathrm{comp}} = 1\,\mathrm{au}$, and assume a circular orbit. We do not assign a specific value to the companion–disk mutual inclination, because the following timescale comparison does not explicitly depend on this angle. This gives a radial communication timescale of $t_{\mathrm{c}} \sim9 \times10^{-4}\,\mathrm{Myr}$ and a global precession timescale $t_{\mathrm{n,global}} \sim1.7\,\mathrm{Myr}$, corresponding to $t_{\mathrm{c}}/t_{\mathrm{n,global}} \sim5 \times10^{-4}$. In this case, radial communication is efficient compared to global precession, so the main ring is expected to respond as a single structure, consistent with the main ring responding coherently as a warped structure. Because the proper-motion-anomaly analysis is subject to a broad mass–orbit degeneracy [52, 9], we also explore an extended parameter range, with $M_{\mathrm{comp}}$ up to 1000 $M_{\mathrm{J}}$ and $a_{\mathrm{comp}}$ up to 7 au. The upper limit on $a_{\mathrm{comp}}$ corresponds to about one-third of the central cavity radius of CQ Tau, $r_{\mathrm{cav}} \sim20$ au [51]. This choice keeps the main ring at radii larger than several times the companion orbit, where the secular torque formulas provide a reasonable approximation [33, 48]. Across this parameter space, $t_{\mathrm{c}}/t_{\mathrm{n,global}}$ remains below unity, indicating that the main ring is expected to respond as a single structure over a wide range of companion parameters, consistent with the warped-ring interpretation in § 5.1.

This scenario requires that the inclined disk–companion configuration can be maintained over a sufficiently long timescale against viscous alignment. Following Bi et al. (2020), we estimate the alignment timescale, $t_{\mathrm{a}}$, over which the mutual inclination between the disk and companion orbit is damped by viscous dissipation, as

where $\alpha$ is the viscosity parameter, and $\omega_{\mathrm{n,global}} = 2\pi/t_{\mathrm{n,global}}$ is the global nodal precession angular frequency. For a conservative estimate of the alignment timescale, we use the disk aspect ratio evaluated at $r_{\mathrm{in}} = 20$ au, $(h/r)_{\mathrm{in}} = 0.119$ [51], and assume $\alpha= 0.01$ when estimating $t_{\mathrm{a}}$ over the parameter space of companion orbital semimajor axis and mass. The resulting map is shown in Figure 11.

Image not extracted

In most of the explored parameter space, $t_{\mathrm{a}}$ exceeds 10 Myr. Adopting a lower viscosity would further lengthen the alignment timescale. Since CQ Tau is of order 10 Myr old [23, 54], this suggests that a warp induced by an inclined companion could persist over a timescale comparable to or longer than the system age. Therefore, the companion-driven warp scenario is plausible as an origin of the warp in the dust disk of CQ Tau.

A companion can also produce other structures, such as a cavity [15]. CQ Tau is known to host a central cavity of $\sim20$ au [51], which is qualitatively consistent with the idea that a companion can carve the disk. At the same time, caution is required when linking the inferred companion directly to the full cavity size. Indeed, Vioque et al. [52] argued that a single companion is unlikely to account for the full cavity size in CQ Tau. This suggests that additional processes would be required to explain the observed cavity, although a companion can still contribute to shaping the overall disk architecture.

A companion can also produce spiral arms [14] and a misaligned inner disk [65]. The companion scenario is therefore qualitatively consistent with the multiple structures observed in the CQ Tau disk, including the NIR spiral arms and the misaligned inner disk [10, 20]. In addition, the shadows cast by such a misaligned inner disk may further induce warps in the outer disk ([64]). Overall, the companion scenario may provide a natural explanation for the ring-scale warp in CQ Tau.

Flyby

A stellar flyby can also induce warps with amplitudes of a few degrees in the outer disk ([27]). Such amplitudes are broadly comparable to the moderate warp inferred for the CQ Tau disk. However, such perturbations are not expected to persist indefinitely. In low-viscosity simulations ([27]), flyby-induced warps survive on a timescale comparable to the viscous damping timescale $t_{\mathrm{damp}} \sim1/(\alpha\Omega_{\mathrm{out}})$ expected from linear bending-wave theory ([34]), where $\alpha$ is the viscosity parameter and $\Omega_{\mathrm{out}}$ is the Keplerian angular frequency at the outer disk radius. Rescaling this damping timescale to the CQ Tau case at $r_{\mathrm{out}} = 70$ au yields values ranging from $\sim10^{4}$ to $\sim10^{6}$ yr over the representative range $\alpha= 10^{-4}$–$10^{-2}$ (Rosotti 2023). While CQ Tau does not show clear interaction with the surrounding environment ([23]), a recent flyby cannot be ruled out as a possible origin of the warp, particularly if the disk viscosity is sufficiently low.

Late Infall

Late infall of material with misaligned angular momentum can also tilt disks and change their angular momentum direction (e.g., [6, 29]). Indeed, several young systems showing warps have also been reported to exhibit asymmetric infall or streamer-like features ([44, 18, 50]). On the other hand, CQ Tau is a relatively evolved system with an age of order 10 Myr ([23, 54]) and shows no clear evidence for ongoing large-scale asymmetric infall. However, a contribution from a past infall event cannot be ruled out, as simulations suggest that infall-induced inner–outer disk misalignment can survive for $\sim0.1$ Myr or longer ([29]).

Magnetic Torque

Magnetic torque arising from a misalignment between the stellar magnetic field and the disk rotation axis has been proposed as one mechanism for producing warps ([31, 21]). These studies mainly describe how magnetic warping is induced in the disk inner region and communicated through a continuously connected disk. In contrast, the CQ Tau system contains a misaligned inner disk ([10]), a cavity, and the warped main ring. Because magnetic torques are expected to act most effectively in the innermost disk regions, they may be relevant to the origin of the misalignment of the sub-au-scale inner disk in CQ Tau, but they do not appear to provide a straightforward explanation for the warp in the main ring.

Summary

In this study, we combined three archival ALMA Band 6 datasets of CQ Tau and analyzed the dust disk by applying the super-resolution imaging method PRIISM to the 1.3 mm dust continuum data. The reconstructed image reveals a compact unresolved inner disk and a bright main ring extending from $\sim20$ to $\sim70$ au. The PRIISM image provides a higher imaging fidelity than that of the CLEAN images and achieves an effective spatial resolution of $0.052^{\prime\prime} \times0.042^{\prime\prime}$, comparable to that of the CLEAN image with a robust parameter of $-0.2$. To examine the underlying disk geometry while reducing the influence of non-centro-symmetric structures, we also reconstructed PRIISM and CLEAN (robust $= -0.2$) images using only the real part of the visibilities. The main scientific results are summarized as follows.

1. Ellipse fitting applied to multiple contours of the real-part PRIISM image with the ellipse center fixed at the stellar position reveals radial variations in PA and $i$ at the inner and outer edges of the main ring. At the inner edge ($\sim21$–26 au), PA and $i$ span $214^{\circ}$–$233^{\circ}$ and $23^{\circ}$–$39^{\circ}$, respectively, while they span $231^{\circ}$–$235^{\circ}$ and $33^{\circ}$–$36^{\circ}$

at the outer edge ($\sim58$–$67$ au). These results obtained from the PRIISM image reconstructed from the combined multi-epoch data are consistent with those obtained from the CLEAN image with robust = $-0.2$ and the PRIISM image using only the single observing epoch. These results support the robustness of the PA and $i$ variations inside and outside the main dust ring.

2. Toy-model tests show that the coplanar ring model cannot reproduce the observed variations in PA and $i$, whereas the fiducial warped ring model provides a better reproduction of the observed radial variations in both quantities. Ellipse fitting including the ellipse center as a free parameter still provides a similar result and the ellipse center is consistent with the stellar position. This suggests that the ring is unlikely to trace an eccentric orbit. These results identify the main dust ring of CQ Tau as warped.

3. The high-pass-filtered real-part PRIISM image shows a pair of centro-symmetric but non-axisymmetric residual patterns. Similar curved negative residuals appear in the high-pass-filtered warped ring model, whereas they are absent from the coplanar ring model. This comparison suggests that high-pass filtering may provide a supplementary tool for identifying disk warps.

4. Together with the previously reported strongly misaligned inner disk and candidate gas warp, the dust warp identified in this work suggests that CQ Tau has a hierarchical three-dimensional disk structure. The disk geometry changes radially from a strongly misaligned inner disk with a tilt of $\sim40^\circ$ to a more weakly warped main dust ring in which $\beta$ decreases from $\sim4^\circ$ at 20 au to $\sim1^\circ$ at 60 au, together with a candidate gas warp with a tilt of a few degrees at $r = 45$ to 155 au.

5. As possible origins of the warp, a stellar flyby and late infall cannot be ruled out depending on the parameters that control the survival timescales of the induced warps. On the other hand, an inclined companion provides a plausible explanation for the warp and is also consistent with several other structures in CQ Tau, including the misaligned inner disk, central cavity, and NIR spiral arms.

Taken together, these results establish CQ Tau as a clear example of a warped dust disk. Combined with the previously reported misaligned inner disk and candidate gas warp, CQ Tau provides an important laboratory for investigating the connection between inner-disk misalignment, ring-scale warps, and planet formation.

Acknowledgments

We thank the referee for the constructive comments and suggestions, which helped improve the manuscript. We are grateful to A. Shoshi for valuable comments. We also thank the ALMA Helpdesk staff for their assistance in providing the calibrated measurement sets. This work is part of the ASIAA Summer Student Program 2024, and S.I. acknowledges support from the Academia Sinica Institute of Astronomy and Astrophysics (ASIAA) in Taiwan. S.I. was supported by the ALMA Japan Research Grant of NAOJ ALMA Project, NAOJ-ALMA-406. M.Y. acknowledges support from JSPS KAKENHI grant Nos. JP26K17220 and JP26K00741 and from NAOJ ALMA Scientific Research grant Code 2022-22B. H.-W.Y. acknowledges support from the National Science and Technology Council (NSTC) in Taiwan through grant NSTC 115-2112-M-001-006- and from the Academia Sinica Career Development Award (AS-CDA-111-M03). S.T. is supported by JSPS KAKENHI grant Nos. JP21H04495 and JP24K00674 and by NAOJ ALMA Scientific Research grant No. 2022-20A. Data analysis was in part carried out on the Multi-wavelength Data Analysis System operated by the Astronomy Data Center (ADC), National Astronomical Observatory of Japan. This paper makes use of the following ALMA data: ADS/JAO.ALMA #2013.1.00498.S, ADS/JAO.ALMA#2016.A.00026.S, and ADS/JAO.ALMA#2017.1.01404.S. ALMA is a partnership of ESO (representing its member states), NSF (USA) and NINS (Japan), together with NRC (Canada), NSTC and ASIAA (Taiwan), and KASI (Republic of Korea), in cooperation with the Republic of Chile. The Joint ALMA Observatory is operated by ESO, AUI/NRAO and NAOJ.

Facility: ALMA

Software: AJISAI v0.1.0 [58], Astropy [3, 4, 5] CASA [11], emcee [19], matplotlib [24], Plotly [41], PRIISM [38, 39], SciPy [53]

Appendix

Procedure for Ellipse Fitting

In this section, we describe the ellipse-fitting procedure applied to the observed emission contours. In the center-fixed case, the ellipse center is fixed at the phase center by setting $(x_c, y_c) = (0, 0)$, where $(x_c, y_c)$ denotes the offset of the ellipse center from the phase center. The free parameters are therefore $(x_c, y_c, a, b, \mathrm{PA})$ for the center-free fit and $(a, b, \mathrm{PA})$ for the center-fixed fit, and the fitting is performed using MCMC with emcee ([19]).

For each contour, we first extract a set of coordinates $(x_i, y_i)$ along the contour from the corresponding continuum image using matplotlib ([24]). Because adjacent contour points are sampled more densely than the spatial resolution of the image, they are not treated as independent measurements. Instead, consecutive contour points are grouped along the contour within a cumulative distance corresponding to the spatial-resolution scale, taken as the major-axis FWHM. Once this distance is exceeded, a new grouping is triggered. This grouping also reduces the influence of noise-induced local fluctuations in the contours. However, as shown in Figure 3, the observed contours already show clear elliptical shapes without significant fluctuations. We therefore adopt the spatial resolution as the grouping scale. For each group, the representative data point $(x_i, y_i)$ is defined as the mean position of the contour points included in that group. To this data point, the positional uncertainties, $(\sigma_{x_i}, \sigma_{y_i})$, are defined as the standard deviations of the grouped coordinates in the $x$ and $y$ directions, respectively. These uncertainties are used as effective positional uncertainties that account for the finite spatial resolution and the local spread of contour points within each spatial-resolution-scale segment.

The derived data points $(x_i, y_i)$ are then converted to polar coordinates $(r_{\mathrm{data},i}, \mathrm{PA}_i)$ with respect to the ellipse center $(x_c, y_c)$. The observed radial distance $r_{\mathrm{data},i}$ is defined as

For the center-free fit, this coordinate conversion is recalculated at each MCMC sampling step. On the other hand, the radial distance of the model ellipse at a given position angle $\mathrm{PA}_i$, denoted as $r_{\mathrm{model},i}$, is derived analytically:

The residual for each derived data point is defined as the difference between the observed radius $r_{\mathrm{data},i}$ and the model radius $r_{\mathrm{model},i}$. Observational uncertainties in the data point positions, $(\sigma_{x_i}, \sigma_{y_i})$, are propagated to estimate the uncertainty in the radial direction, $\sigma_{r_i}$, based on the standard first-order error propagation formula. The log-likelihood function for the ellipse fitting is defined as

For the center-free fit, $\sigma_{r_i}$ varies with $(x_c,y_c)$ and is therefore updated at each MCMC sampling step. Based on this likelihood function, we estimate the posterior distributions of the free parameters and adopt the median values and the 16th–84th percentile ranges as the representative values and uncertainties. The inclination $i$ is derived from the sampled ellipse parameters, $a$ and $b$, as $i=\arccos(b/a)$.

To assess the goodness of fit, we also compute the reduced chi-square,

where $N_{\mathrm{data}}$ is the number of grouped data points and $N_{\mathrm{par}}$ is the number of free parameters. This quantity is evaluated for the median posterior parameter set.

Construction of the Coplanar Model

To construct an infinitesimally thin coplanar ring model, we first deprojected the real-part PRIISM image using $\mathrm{PA}_0$ and $i_0$ and then computed an azimuthally averaged radial intensity profile (Figure B1). The uncertainty at each radius was estimated as the standard deviation over all azimuthal angles.

Image not extracted

Next, we modeled the radial intensity profile as a sum of Gaussian components,

where $A_k$, $\mu_k$, and $\sigma_k$ are the peak intensity, central radius, and radial width of the $k$th Gaussian component, respectively, and $N_{\mathrm{G}}$ is the number of Gaussian components. This multi-Gaussian profile was convolved with

a one-dimensional Gaussian whose FWHM was set to the geometric mean of the effective spatial resolution of the PRIISM image, and the convolved profile was fitted to the observed radial profile using the curve_fit function in scipy.optimize (Figure B1). The fitting was performed in the radial range of $0.1^{\prime\prime}$–$0.6^{\prime\prime}$ to avoid contamination from the inner compact disk emission. The number of Gaussian components was determined as the minimum number for which the absolute residuals between the observed and convolved model profiles remained below $I_{\mathrm{DT}}$ throughout the fitting range. This criterion resulted in $N_{\mathrm{G}} = 4$, and the corresponding best-fit parameters are summarized in Table B1.

Component	$A_k$ [Jy arcsec$^{-2}$]	$\mu_k$ [arcsec]	$\sigma_k$ [arcsec]
Component 1	0.150	0.189	0.0225
Component 2	0.268	0.258	0.0369
Component 3	0.322	0.355	0.0478
Component 4	0.0194	0.535	0.0414

Then, the best-fit intrinsic profile was azimuthally rotated to construct a two-dimensional face-on axisymmetric ring model. The face-on model was subsequently projected onto the sky plane using $\mathrm{PA}_0$ and $i_0$ while conserving the total flux, as appropriate for the observed optically thin continuum emission.

Dependence of the High-pass-filtered Images on the Smoothing Scale

Figure C1 compares the high-pass-filtered images obtained with smoothing FWHMs of $0.05^{\prime\prime}$, $0.15^{\prime\prime}$, and $0.20^{\prime\prime}$, while Figure 10 shows the images obtained with the adopted FWHM of $0.1^{\prime\prime}$. Although the detailed appearance and residual intensity vary with the smoothing FWHM, the qualitative differences among the real-part PRIISM image, coplanar ring model, and warped ring model remain unchanged.

Image not extracted