Revealing asymmetry on midplane of proto-planetary disc through modelling of axisymmetric emission: methodology

This study proposes an analytical framework for deriving the surface brightness profile and geometry of a geometrically-thin axisymmetric disc from interferometric observation of continuum emission. Such precise modelling facilitates the exploration of faint non-axisymmetric structures, such as spirals and circumplanetary discs. As a demonstration, we simulate interferometric observations of geometrically-thin axisymmetric discs. The proposed method can reasonably recover the injected axisymmetric structures, whereas Gaussian fitting of the same data yielded larger errors in disc orientation estimation. To further test the applicability of the method, it was applied to the mock data for $m=1,2$ spirals and a point source, which are embedded in a bright axisymmetric structure. The injected non-axisymmetric structures were reasonably recovered except for the innermost parts, and the disc geometric parameter estimations were better than Gasussian fitting. The method was then applied to the real data of Elias 20 and AS 209, and it adequately subtracted the axisymmetric component, notably in Elias 20, where substantial residuals remained without our method. We also applied our method to continuum data of PDS 70 to demonstrate the effectiveness of the method. We successfully recovered emission from PDS 70 c consistently with previous studies, and also tentatively discovered new substructures. The current formulation can be applied to any data for disc continuum emission, and aids in the search of spirals and circumplanetary discs, whose detection is still limited.


INTRODUCTION
High-resolution imaging of proto-planetary discs by the Atacama Large Millimeter/submillimeter Array (ALMA) has revolutionized our view of planet formation.The spatial resolutions, down to ∼ 0.01 ′′ , have enabled the identification of rings and gap structures in the disc continuum emission (ALMA Partnership et al. 2015;Andrews et al. 2018;Huang et al. 2018a).The mechanisms for invoking such annular structures have been under debate, and multiple explanations have been provided, such as a disc-planet interaction (Lin & Papaloizou 1979;Goldreich & Tremaine 1980;Dong et al. 2018), snowline (Zhang et al. 2015), sintering (Okuzumi et al. 2016), secular gravitational instability (Takahashi & Inutsuka 2014), and non-ideal MHD effects (Flock et al. 2015).Nevertheless, the disc-planet interaction is independently supported by recent discoveries of kinematic signals of embedded planets (Pinte et al. 2018(Pinte et al. , 2020(Pinte et al. , 2023)).Here, the signals, usually referred to as velocity "kinks", appear as localized deviations from the Keplerian motion of discs in velocity maps, and there have been dozens of such detections.
In contrast to multiple detections of rings, gaps, and kinks, the ★ E-mail: masataka.aizawa.sw35@vc.ibaraki.ac.jp direct detection of embedded planets and the circumplanetary discs within gaps remains limited.Currently, there are a few cases presenting the robust detection of embedded planets in gaps, such as, PDS 70 (Keppler et al. 2018;Wagner et al. 2018;Benisty et al. 2021), AB Aur (Currie et al. 2022), HD 169142 (Reggiani et al. 2014;Hammond et al. 2023), and MWC 758 (Wagner et al. 2023).These embedded planets have been somewhat selectively identified around (pre-)transitional discs with the ages greater than 4 Myr.In addition, Andrews et al. (2021) searched for circumplanetary discs in gaps of DSHARP discs (Andrews et al. 2018), whose typical age is 1 Myr; however, no compelling case remains.
Another supportive evidence for the planetary hypothesis can be obtained from a planetary spiral, which is excited by planetary gravity.However, despite the discovery of many spirals in discs (Muto et al. 2012;Grady et al. 2013;Huang et al. 2018b), only the spirals in systems with detected embedded planets are highly likely to be produced by the planetary gravity.The other spirals can also originate from non-planetary mechanisms, for example, gravity from stellar companions, gravitational instability or stellar flyby.Indeed, the spiral in HD 100453 A is a notable example, where its M-dwarf companion is thought to be the driver for the structure (Wagner et al. 2015;Dong et al. 2016).
Recently, Speedie & Dong (2022) searched for planetary spirals in 10 discs with kinematic signatures of embedded planets; however, they did not find any conclusive case.The limited detection of circumplanetary discs and spirals might indicate that the current sensitivities are still inadequate for most of the cases, or certain annular structures may have been caused via non-planetary mechanisms.
Nevertheless, an observational effort to search for circumplanetary discs and spirals is still essential for testing the planetary hypothesis.However, one obstacle to their detection is that these signals appear as small perturbations to the bright disc emission.This is similar to the situation of searching for a needle in a haystack.Thus, the precise modelling of the background disc emission is necessary to subtract its contribution.
One reasonable approximation for the background emission is an axisymmetric disc.However, the modelling of a vertically extended axisymmetric disc is difficult because of the geometric effect (e.g., shadows) and additional parameters to be solved (e.g., radial disc heights).Such geometric effects are significant at optical to infrared wavelengths, because of their sensitivity to vertically extended small dust grains.On the other hand, the disc emission at the radio band is well approximated by a geometrically thin disc model with no vertical structure owing to its sensitivity to large grains, which are more settled to the disc midplane.This approximation is to a certain extent justified by the observations of many clear concentric structures of discs in radio continuum emissions (e.g., ALMA Partnership et al. 2015;Andrews et al. 2018).
Recently, assuming a geometrically thin and axisymmetric disc, Jennings et al. (2020) proposed an algorithm for deriving a radial surface brightness profile of continuum emission in the radio band.Their method can resolve annular structures finer than the standard imaging technique, CLEAN, using simulated and real data (Jennings et al. 2020(Jennings et al. , 2022a,b),b).Using their method, Andrews et al. (2021) also searched for circumplanetary disc emission from DSHARP data by subtracting the axisymmetric structures of discs.
Although the method proposed in Jennings et al. (2020) can reasonably recover the brightness profile of the disc, there is still room for improvement.Specifically, in their modelling, the geometric parameters of the axisymmetric model, such as an orientation and a central position of a disc, and hyperparameters for the Gaussian Process model that is used to prevent overfitting must be fixed.This limitation can be problematic, because a slight change in geometric parameters of the disc model, for example, several mas in a central position, can introduce the apparent non-axisymmetric structures in the image with their axisymmetric structure being subtracted (Andrews et al. 2021).Such false structures can be degenerate with the real structures, and may disrupt the detection of spirals and circumplanetary discs.Till date, in the frameworks proposed for axisymmetric disc models, these geometric parameters or hyperparameters for regularization have been manually tuned (Jennings et al. 2020;Andrews et al. 2021), or estimated through parameterized models (Zhang et al. 2015;Tazzari et al. 2017;Jennings et al. 2020;Kanagawa et al. 2021) and image-based analysis (Huang et al. 2018a).
Therefore, this study proposed an analytical framework to derive all of the parameters for a geometrically-thin axisymmetric disc and hyperparameters for the Gaussian Process kernel from observations assuming a geometrically thin disc.This study employed the formulation for the inverse modelling presented in Kawahara & Masuda (2020), which recovered a planetary surface map and its spin and orbital geometry from planetary reflection light.Interestingly, the problem considered in that study is mathematically similar to the current problem.Thus, we can apply their methodology for the current problem.We demonstrate the feasibility of the current method via its application to mock and real data.In addition, we perform injection and recovery experiments for circumplanetary discs and spirals, and discuss the applicability and limitations of the method.
The remainder of this paper is organized as follows.In Sec 2, the method for estimating an axisymmetric structure from visibilities is formulated.In Sec 3, the methodology for studying nonaxisymmetric features in a residual image is discussed.In Sec 4, the feasibility of the proposed method for mock observations is investigated.In Sec 5, circumplanetary discs and spirals are injected to mock data, and the ability to recover the structures is examined.In Sec 6, our method was also applied to the real data for Elias 20 and AS 209.Sec 7 presents an analysis of the PDS 70 data, demonstrating the recovery of circumplanetary emission.Finally, the conclusions and future improvements are presented in Sec 8.

Model
We model a geometrically thin axisymmetric disc, whose parameters are composed of the brightness profile  and its geometry .Here,  = { (  )} for  = 1, 2, ...,  is a vector obtained by discretizing the radial surface brightness profile  (), and   is the -th collocation point of the Fourier-Bessel series: where  out is the outer boundary of  () with  ( >  out ) = 0, and  0 is the -th root of the zero-order Bessel function of the first kind  0 ();  0 (  0 ) = 0.The disc geometry is specified by the disc orientation and the central position:  = (Δ cen , Δ cen , cos , PA).Here, the disc centre (Δ cen , Δ cen ) is assumed to be shifted from a phase centre of the observation, and the disc orientation is specified by the position angle PA and inclination .PA is defined as the angle of the major axis measured counter-clockwise from the north direction, and  is defined as the angle between the line of sight and the axis normal to the disc plane.Further, the formulation adopts cos , which represents the aspect ratio for the disc ellipse.

Face-on disc
As the simplest problem, we consider the face-on disc with no positional offset from the phase centre;  = (Δ cen = 0 ′′ , Δ cen = 0 ′′ , cos  = 1, PA = 0).Visibilities for an azimuthally symmetric geometrically thin disc with a brightness profile  () are expressed as the Hankel transformation (Thompson et al. 2017;Jennings et al. 2020): where  (, ) is the complex visibility at a spatial frequency (, ),  cos =1 = √  2 +  2 is the deprojected spatial frequency for the faceon disc, and  0 is the zero-order Bessel function of the first kind.

Inclined disc
The computations of model visibilities can be simply extended to the case of the inclined disc.For the calculation, we prepare three different coordinates (, ), ( ′ ,  ′ ), and ( ′′ ,  ′′ ); Fig. 1 shows a schematic of the coordinates and disc.Here, (, ) is the observational coordinate system with the phase centre being the origin, ( ′ ,  ′ ) is the shifted coordinate system with the disc centre being the origin, and ( ′′ ,  ′′ ) is the deprojected coordinate system according to (PA, ).(, ) and ( ′ ,  ′ ) are related as (, ) = ( ′ + Δ cen ,  ′ + Δ cen ), and ( ′ ,  ′ ) and ( ′′ ,  ′′ ) are converted as follows: The brightness profiles for the coordinates  (, ),  ′ ( ′ ,  ′ ), and  ′′ ( ′′ ,  ′′ ) are assumed to be invariant to coordinate transformations: This assumption is a mathematical hypothesis for modeling purposes, and it is not necessarily consistent with the physical picture of what would actually occur if the disc were observed from the face-on view 1 .Assuming the brightness profile in a deprojected frame  ′′ ( ′′ ) with  ′′ ≡ √︁  ′′2 +  ′′2 , we aimed to compute the model visibilities.First, we transform equation (2) as follows: where we define ( ′ ,  ′ ) as follows: Using equation (3), equation ( 9) is reduced as follows: 1 Nevertheless, the assumption is physically consistent when the disc is a geometrically-thin, optically thick disc.
For  real and  imag , which are defined as real and imaginary parts of , respectively, we derive a linear equation as follows: where  ′ ,  real , and  imag are defined as follows: (16)

Inverse modelling
On the basis of the forward modelling in Sec 2.2, we consider inverse modelling for brightness profile  and geometric parameters  = (Δ cen , Δ cen , cos , PA) from visibilities.

Data with Gaussian noise
Let us assume that there are  visibilities  obs = { obs, (  ,   )} ( = 1, 2, ..., ), with the noise of each visibility data obtained via the standard deviation of  obs = { obs, }.We separate the real and imaginary parts of the visibilities as  real ≡ ℜ(  obs ) and  imag ≡ ℑ(  obs ), respectively.We assume that observational noises for  real and  imag obey the multivariate normal distribution.Here, the multivariate normal distribution for a -dimensional random vector  is expressed as follows: where  is the mean vector and  is the covariance matrix.For  real and  imag , we assume the covariance matrix to be  =  obs ≡ { ,  / 2  }.

Gaussian prior on brightness profile
In interferometric observations, a brightness profile  is generally under-constrained.If it is directly obtained by solving a linear equation, for example, equation (4), the solution can diverge, leading to overfitting.This can be avoided by introducing regularization in optimization.Specifically, we assume the Gaussian Process kernel with a set of hyperparameters  = (, ) on  as follows: where we adopt the radial basis function (RBF) kernel   (), which is expressed as; The disc is inclined with cos  and oriented by PA with respect to the -axis in the projected frame.The origin of ( , ) is set to be the phase centre, and (  ′ ,  ′ ) is shifted from ( , ) by (Δ cen , Δ cen ).(  ′′ ,  ′′ ) is the coordinate for the deprojected frame, and the brightness profile is assumed to be only radially dependent on  ′′ ( ′′ ).
where  determines the relative weight of the prior information, and  determines the spatial scale for regularization.The second term in Eq (19) is the small identity matrix, which stabilizes the computation of the inverse matrix Σ  () (Kawahara et al. 2022).We adopted  = 10 −8 in this paper.
This prior tends to promote smooth solutions, and can aid in preventing overfitting.A previous study by Jennings et al. (2020) adopted parameters for controlling the Gaussian Process prior on powers of the brightness profiles in the frequency domain.However, in this study, we adopted the parameters for regulating the profile in the spatial domain to render the discussion simpler.

Likelihood function for geometry and hyperparameters
Following Kawahara & Masuda (2020), we derive the likelihood function ( | , ).We briefly overview their formulation and apply it to the current problem.For further details, see equations ( 14)-( 21) and Appendix A.1-A.3 in the previous study.

Posterior distribution for geometry and hyperparameters
The posterior distribution for ( ,  | ) is given by the Bayes' theorem as follows: Based on ( ,  | ), we draw samples for ( , ) using an MCMC sampler.The prior distribution ( , ) is assumed as follows.We assume uniform distributions U (0, 1) and U (0, ) for cos  and PA, respectively.Further, we consider U (−1 ′′ , 1 ′′ ) for Δ cen and Δ cen .In addition, we assume the log-uniform prior from 10 −4 to 10 4 for  and a uniform prior U (0.01 ′′ , 0.15 ′′ ) for .We checked that the likelihood is exceedingly small when  and  fall outside their respective prior ranges, thus confirming the sufficient broadness of the prior.As the computation for ( | , ) can be time-consuming, we transform the equation to a more efficient form (further details in Appendix A).
Moreover, the computation also relies on the number of data points.In Appendix B, we introduce and discuss the use of data binning to reduce the computation.Specifically, we prepare two-dimensional linear and log grids, apply them to the simulated visibilities, and compare the binning errors.The binning with the log grid was found to yield overall low errors.A more detailed discussion is presented in Appendix B. In the following analysis, we simply adopt the log grid for the analysis.

Posterior distribution for all of parameters
We can draw samples for {, ,  } by using the conditional formula for the joint probability: (, ,  | ) = (| , , ) ( ,  |).Specifically, we can draw a sample (  †,  †) from ( ,  |) as done in Sec 2.3.4,and subsequently take a sample for  from (| ,  †,  †) using equation ( 28).We can iterate this procedure to make samples, and this sampling is indeed equivalent to drawing sample from the joint distribution (, ,  | ).Using the samples for {, ,  }, we can also compute statistics for parameters (e.g., mean and standard deviation for ).

Difference in formulation between current study and frank
We here highlight the key differences between our approach and that of frank.Note that the direct comparison of the recovered profiles is also given in Section 6.While both methods vary in their regularization strategies, the distinct difference is that our approach solve all parameters, including brightness profiles, geometric parameters, and hyperparameters, whereas frank is limited to solving brightness profiles.In this sense, our formulation can be seen as a natural extension of that of frank.
One advantage of our method is that it can optimize geometric parmaeters and hyperparameters directly from the data, and thus it is less susceptible to human biases caused by manual tuning.Moreover, while the current paper is limited to the simplest model with single frequency and single source, the methodology itself can be easily extended to more complex problems with multi-frequency data or multiple sources.In such cases, manual tuning of optimal parameters can be both challenging and less reliable, making our approach more effective.
On the other hand, the disadvantage of our approach is that it is more computationally demanding than frank.This is because we attempted to solve all parameters, including non-linear parameters, rather than deriving only brightness profiles like frank.Furthermore, our current formulation cannot support imposing a nonnegative condition on brightness profiles, a feature available in frank.This is because the analytical expression for the marginalization over  becomes inapplicable under the non-negative condition.Although we can implement the condition by considering the same problem as frank, the current study prioritizes simultaneous fitting of all parameters, so we do not implement the function.

EXTRACT ASYMMETRIC FEATURES
This section summarizes the process of extracting and quantifying the non-axisymmetric components from observations using the derived parameters in Sec 2.3.The method is applied to data analyses in Sec 4 and 5.

Making residual image
After subtracting the axisymmetric model from visibilities, we expect the residual visibilities to contain only non-axisymmetric components.Therefore, imaging using these residual visibilities can reveal the non-axisymmetric component of the disc (Jennings et al. 2020;Andrews et al. 2021).Based on this principle, we created images to study the non-axisymmetric components.
Drawing samples from a posterior distribution as described in Sec 2, we selected the maximum a posteriori probability (MAP) estimate for  and .Subsequently, we randomly drew a brightness profile using equation (28).For the chosen parameters and the brightness profile, we computed the model visibilities of an axisymmetric disc, and subtracted them from observed visibilities.Then, the phase centres of measurement sets were aligned with the disc centres using fixvis in CASA.
For the shifted and subtracted measurement set, we created images using the CLEAN algorithm, which is implemented as tclean in CASA.In this study, the pixel scale for the image was set as 6 mas, and the image size was 1000 × 1000 pixels.In the imaging process, Although it is typical to adopt 0 iteration for CLEAN in making the residual map (e.g., Jennings et al. 2020), we found that employing non-zero iterations for CLEAN result in slight or moderate improvement in the quality of the residual image, especially when there remain substantial residuals.Thus, in this paper, we implemented thresholds with nsigma= 3.5, which stops the iterations if the maximum residual on the image is below 3.5 times the rootmean-square error.Further, we adopted the Briggs weighting with a robust parameter of 0.5 and without UV-taper.Appendix C presents a discussion on the effects of changing robust parameters on the recovered residual images, and 0.5 was determined to be the balanced choice.The brightness of the output image from CASA was measured in Jy beam −1 , where 'beam' is the nominal area where the brightness is defined2 .The standard deviation of the residual image was also measured outside the disc emission.

Deprojection of residual image
In addition to the residual image in an observational frame, we created the deprojected residual image using the disc geometry.Spatial frequencies of the shifted measurement set were converted using the geometric parameter, (, ) → ( ′ ,  ′ ) in equation (10).Consequently, the image was created in the same manner as Sec 3.1.1.This operation aligns the major axis of the image with the -axis (north direction), and stretches the image along the  axis (east direction).This corresponds to the conversion in Figure 1 from a deprojected to a projected frame.
In the imaging process, we should be cautious that the transformation (, ) → ( ′ ,  ′ ) apparently decreases the brightness in Jy sr −1 by a factor of cos , while brightness in Jy beam −1 remains unchanged before and after deprojection.This is because the transformation increases the disc area (and beam) while the total flux is unchanged.Throughout the paper, we consistently presented the brightness of images in Jy beam −1 , which remains constant, and we did not apply any correction by a factor of cos .
Although widely used, the current simple deprojection methods provide the true image of the disc viewed face-on only in the case of an infinitesimally thin disc.The other geometric effects such as the vertical structures or shadows cannot be incorporated with our simple operations.

Mode decomposition of residual image in polar direction
Spirals are often characterized by dominant modes in a polar direction.We can decompose  (, ) in a polar coordinate (, ) as follows (e.g., Binney & Tremaine 2008): where   () is the amplitude at the -th mode, and   () is a phase shift for the -th component.Note that  0 () corresponds to the brightness profile.The amplitudes and phases can be derived from the residual images  (, ) as follows: where atan2 is the 2-argument arctangent, and (  ,   ) is defined as follows: The values of (  ,   ) can be computed for the real data.In this study, we computed them by interpolating a residual image.

Extraction of odd and even symmetric components of images
Odd-symmetric and even-symmetric components of an image can be decomposed by using either of each imaginary or real part of the data.
We start with decomposing an image into even-and odd-symmetric components.
where we define which satisfy  even (, ) =  even (−, −), and  odd (, ) = − odd (−, −).However, using equation ( 2), the real and imaginary parts are connected to  even (, ) and  odd (, ) as follows: Thus, we can create  even (, ) and  odd (, ) using either the real or imaginary part of the data.Notably, previous studies have already used the imaginary part of the data to extract an odd-symmetric component (Hashimoto et al. 2021;Kanagawa et al. 2021).

Injection and recovery test
As a test case, we applied our method to simulated data for an axisymmetric disc to examine the ability to recover input parameters.We used two radial profiles for our test.The continuum brightness profiles for AS 209 and WaOph 6 of DSHARP observations3 (Andrews et al. 2018) were used.AS 209 is chosen as the most structured disc in a radial direction, and WaOph 6 is among  = 2 spiral discs in the DSHARP sample.Subsequently, we assume the model discs have the geometric parameters of (Δ cen , Δ cen , cos , PA, ) = (0 ′′ , 0 ′′ , 0.75, 45 • ).We normalize the brightness profiles to render the total fluxes of the AS 209 and WaOph 6 models as 0.288 and 0.161 Jy, respectively.
To produce mock data with realistic UV-coverage and noise, we incorporated them using the actual data of DSHARP observations.We downloaded the self-calibrated continuum measurement set for AS 209 and WaOph 6.The data were averaged at each spectral window to produce a single channel data.Thereafter, time averaging was applied for 30 s.A new measurement set was created using the averaged data to obtain the UV-coverage and the weight at each spatial frequency.The signals of the model disc at each spatial frequency were calculated via Fourier transforms of the model surface brightness distribution.Consequently, the Gaussian noises were added to the model visibility according to the the weights.During the analysis, we used tb.getcol in CASA to export the measurement sets to handle them numerically in python4 .
The data weights recorded in the measurement sets were overestimated for the entire DSHARP data.Following Appendix E in Hashimoto et al. (2021), we deprojected the real and imaginary parts of visibilities, and compared the recorded standard deviations and the deviations of observed visibilities from the binned visibilities in the deprojected spatial frequencies.Consequently, we confirmed that the overestimation existed in both the real and imaginary parts of the visibilities, and the degree of overestimation was within ∼ 10% for both.Therefore, we manually reduced the recorded weights by 3.44 and 3.66 for AS 209 and WaOph 6, respectively.
To reduce the computational time, we binned the data using a logarithmically spaced grid with  bin = 500, where the number of grid points was (2 bin + 3) 2 .The boundaries of the grids are defined by the parameters ( min ,  max ), which are introduced in Appendix B. We set ( min ,  max ) to (10 2 , 10 7 ) in our analyses.This binning decreased the number of data points by an approximate factor of 10.Appendix B presents the details of data gridding using a log grid, and we show that the choice of  bin = 500 is acceptable.
Using equation ( 33), we drew samples from the posterior distribution for (, log 10 , Δ cen , Δ cen , cos , PA).In sampling, we used emcee, which implemented the affine-invariant ensemble sampler for Markov Chain Monte Carlo (Foreman-Mackey et al. 2013).The initial parameters for emcee were set to be close to the input parameters.Here, 16 walkers were prepared and then evolved for at least 8,000 steps.For the drawn samples, we discarded the initial 1,250 steps during the burn-in phase.The convergence is predominantly achieved within 2,000-4,000 steps according to the auto-correlation time analysis.Specifically, we have verified that all samples meet the condition  > 50, where  represents the step size and  represents the integrated autocorrelation time for the chain 5The total computational time for the single disc was approximately 12 hours using 8 cores.This represents a significantly higher computational demand in contrast to frank's approach.This stems from our strategy to retrieve all parameters in a simultaneous manner, as opposed to frank's approach of fixing parameters except for the brightness profiles.
Figs. 2 and 3 show the recovery and injection test results for the simulated data.The upper-left panels shows the posterior distribution of the parameters for the two simulated cases.All the geometric parameters were reasonably recovered by the current method.
The analyses resulted in two different length parameters;  ≃ 0.04 ′′ for the simulated case of AS 209, and  ≃ 0.1 ′′ for the case of WaOph 6.The differences in  were largely due to the differences in the injected brightness profiles; the injected profile for WaOph 6 was smoother than that of AS 209.The UV-coverage would also affect the length scale, but it remains unclear in the current two cases.The length scale parameter for WaOph 6 was larger than the beam size (0.091 ′′ , 0.037 ′′ ), whereas it was smaller than the beam size (0.076 ′′ , 0.040 ′′ ) for AS 209.Note that the beam sizes reflect the UV-coverage, but they vary with assumed parameters in CLEAN; for example, robust parameter and UV-taper.In Appendix D, we further study the variations in  by stretching the injected brightness profiles and UV-coverage.In conclusion, the optimal length scales  appear to be determined by multiple factors, including the brightness profiles and UV-coverage (beam size).
The upper-right panels in Figs. 2 and 3 compare the injected brightness profiles with 10 samples for  using the method presented in Sec 2.3.5.The radial profiles were also adequately recovered by the current method.However, the residuals for the brightness profiles were not perfectly consistent with the zero line for both cases, although an oscillating feature was observed in the radial direction.The oscillating length scales for the observed residuals indeed corresponded to the derived length scales , which determined the characteristic scales that could be resolved.Moreover, we also identified the large residuals at the very innermost parts owing to the steep changes in the fluxes in the radial direction.The oscillating features have been also reported in Jennings et al. (2020), which adopted regularization in the frequency domain, and such residuals might be inevitable in the Gaussian Process.
The lower-left panels in Figs. 2 and 3 show the simulated and model visibilities.The model visibilities were computed by choosing the MAP solution for the geometry and hyperparameters, and the brightness profile was randomly sampled.The residual visibilities are shown in the lower panels.The models were well consistent with the simulated visibilities from the low to high spatial frequencies.However, we observed that the powers of the model visibilities at high frequencies were forcibly suppressed.This cut-off scale was determined by the optimized length scale , which determined the scales below which the oscillating features appeared in the brightness profiles.At the low spatial frequencies, particularly at  < 0.3 M, we also found oscillating residuals with a length scale of ∼ 0.05 M.This length scale corresponded to the radius for the outer boundary  out = 2 ′′ , which limited the model's applicability to the large-scale structure.
In the lower-right panels in Figs. 2 and 3, we show the residual images, which are produced by the method presented in Sec 3.1.The residual images did not exhibit noticeable features, demonstrating that we reasonably estimated the injected parameters.

Residual images constructed from shifted geometric parameters
To demonstrate how incorrect determination of geometric parameters affects the resultant residual images, we shifted each geometric parameter before subtracting axisymmetric models from the visibilities.For the simulated case of AS 209, we shifted each geometric parameter from the MAP estimate as follows: Δ cen = 5 mas, Δ cen = 5 mas, Δ cos  = 0.01, ΔPA = 2 deg.For the simulated case of WaOph 6, we shifted the parameters as follows: Δ cen = 2 mas, Δ cen = 2 mas, Δ cos  = 0.03, ΔPA = 3 deg.To simulate the realistic situation as much as possible, we assume these incorrect geometric parameters with the hyperparamters from the MAP solution, and draw a brightness profile from the posterior distribution, which nevertheless gives the profile similar to our best-fit model.With the subtracted visibilities made with the new models, we created the residual images using the CLEAN following the method presented in Sec 3.1.We also created images from unsubtracted visibilities for comparison.For that, we adopted 0.04 mJy as the CLEAN threshold and used the multiscale CLEAN algorithm with scale parameters of [0,30,120,360,720,1440] mas (Rau & Cornwell 2011).Fig. 4 and 5 show the unsubtracted and residual images with their geometric parameters being shifted.They are all shown in observed frames.The residual images based on the best-fitting models did not exhibit any noticeable structure.This demonstrates the feasibility of the current method.However, the residual images with shifted geometric parameters exhibited coherent patterns.These residual images obtained from the current study were overall consistent with those presented in Andrews et al. (2021), wherein the dependence of residual images was studied by shifting the geometric parameters.As shown in the figures, the difference in the injected brightness profiles changed the residual images.We highlight two outer rings for AS 209 by black dashed ellipses in Fig. 4, and we find that the ring locations indeed corresponded to the boundaries, where the signs of the residuals were flipped.However, the brightness profile for WaOph 6 was more continuous than AS 209, and the residual images were less structured.
Fig. 6 shows  =1,2 () for the residual images in the simulated case for WaOph6.The shifts in the central positions introduced residuals with the  = 1 component, whereas those in cos  and PA introduced the residuals with  = 2 component.As discussed in the later sections, these residuals can be degenerate with real signals (e.g., spirals); thus they can introduce biases for estimating such non-axisymmetric structures in discs.
As a comparison with the proposed method, we also applied the Gaussian fitting to the visibilities to estimate the geometric parameters.Here, the Gaussian fitting of visibilities corresponds to the fitting of the Gaussian function on the image plane, and it is frequently used for estimating a geometry of a disc.The model was specified by six parameters, where four parameters, (Δ cen , Δ cen , cos , PA) are common to our model.We implemented a posterior sampling for the Gaussian model using emcee (Foreman-Mackey et al. 2013).
The first and second data from the left in each panel of Fig. 7 show the injected and recovered geometric parameters obtained from the current method and the Gaussian fitting.The central positions were well recovered in both of cases.However, the parameters cos  and PA obtained from the Gaussian fitting deviated from the injected values; Δ cos  ≃ 0.01 − 0.05 and ΔPA ≃ 1 − 5 deg.These deviations can be problematic when searching for faint signals (e.g., circumplanetary discs and spirals) because they introduce fake non-axisymmetric structures in the residual images as shown in Figs. 4 and 5.

INJECTION AND RECOVERY TEST FOR SPIRALS AND CIRCUMPLANETARY DISC EMBEDDED INTO AXISYMMETRIC EMISSION
As a simple extension to Sec 4, we considered an additional nonaxisymmetric perturbation to the symmetric disc.We injected spirals and a point source into an axisymmetric model and simulated the visibilities incorporating noises.Subsequently, we applied our method to extract the axisymmetric model and construct residual images, which were then compared against the injected structures.

Injection of spirals
We injected spiral structures into an axisymmetric structure assuming the same observational setup as the simulated case for WaOph6.The brightness profile was assumed to be the same as the case for WaOph6.In this simulation, we assume zero brightness outside  = 1.11 ′′ , beyond which the brightness in the literature can become negative.Three different models were considered; an oddsymmetric spiral with  = 1 component, an even-symmetric spiral with  = 2 component, and the combined image for both spirals.Specifically, we assumed perturbations in the form of Δ (, , ) = Δ  () cos(( −   ())), where we adopted an Archimedean spiral   () =  + .We assumed (, , ) = (1, 20, 1) for the odd-symmetric spiral and (, , ) = (2, −10, −/6) for the evensymmetric mode.The amplitudes   () were assumed as follows: where the beam area was assumed to be 0.00354 [arcsec] 2 .For the calculation of model visibilities, we used functions for the nonuniform fast Fourier transform from PyNUFFT (Lin & Chung 2017; Lin 2018).We adopted the size of the image grid   = 1024, the size of the oversampled Fourier grid   = 2048, and the size of the interpolator   = 6.Fig. 8 shows the amplitudes of the injected brightness profile and spirals.The amplitudes of the injected spirals were at most ∼ 30 per cent of that of the axisymmetric part of the surface brightness, and they were comparable to the amplitudes of the observed spirals with  = 2 mode (e.g., WaOph6 and IM Lup) such that the current simulation indeed considered the realistic situation.The left panels in Fig. 9 show the injected images for the three cases.

Comparison of injected and recovered models
After simulating the visibilities, we applied our method to estimate the axisymmetric structures.We selected the MAP estimate on geometric parameters, drew a brightness profile, and subtracted their contributions from the observations.The residual visibilities were then used to construct the residual images using CLEAN.
For the comparison, we also created the residual images by adopting geometric parameters estimated from the Gaussian fitting to the visibilities.Specifically, we assumed geometry with Gaussian fitting and hyperparameters from our method, and we drew a brightness profile.Subsequently, the brightness model was subtracted from the visibilities, and the residual visibilities were imaged with CLEAN.
The residual images are shown in Fig. 9.In the second column, we show the recovered residual images obtained via the current method.The last column shows the residual images created from the geometry given by the Gaussian fitting.The residual images from our method were reasonably consistent with the injected maps, whereas those based on the Gaussian fitting exhibited larger residuals.However, the innermost structures of the discs at  < 0.2 − 0.3 ′′ were not perfectly recovered even with our method.Specifically, the excess in the amplitude of the spiral was observed at  < 0.3 ′′ for the oddsymmetric spiral.Moreover, no clear spiral structure was observed at  < 0.2 ′′ for the even-symmetric spiral.These inconsistencies arise from the biases in the estimated geometric parameters, as discussed next.
The comparison of geometric parameters for three cases is shown in Fig. 7.For every parameter, larger deviations from the injected parameters were observed in the case of the Gaussian modelling than that using the proposed method.However, our method was still affected with biases; (Δ cen , Δ cen ) for the odd-symmetric spiral, and (PA, cos ) for the even-symmetric spiral.In the case of the combined spiral, the shifts in all the geometric parameters roughly corresponded to the summation of shifts in the odd and even-symmetric spirals.These shifts are reasonable because, as discussed in Sec 4, the shifts in (Δ cen , Δ cen ) introduce the residuals with  = 1 mode, which can be degenerate with the injected odd-symmetric spiral.Similarly, the shifts in (cos , PA) introduce the residuals with  = 2 mode, which can be degenerate with the even-symmetric spiral.These degeneracies produce the biases in estimating geometric parameters.Consequently, these biases impede recovery of injected non-axisymmetric structures, as shown in Fig. 9.
The left panels in Fig. 10 compare the residual image with the injected odd-symmetric spiral.The phase for the spiral was reasonably recovered in the residual image.The upper right and lower panels show the recovered radial phases  1 () and the amplitudes  1 ().The amplitudes and phases were reasonably recovered in the range of 0.3 ′′ <  < 0.5 ′′ , where the spiral amplitude was large.However, the amplitudes were clearly overestimated at inner radii  < 0.3 ′′ .This overestimation was due to the biased central position, and is consistent with the failure of the recovery of the spiral at the innermost part in Fig. 9.
The right panels in Fig. 10 show the comparison in the case of the even-symmetric spiral.The phases were reasonably recovered in most of the radial range; however, the deviations were large at the innermost part  < 0.1 ′′ .The amplitudes were well recovered at 0.3 ′′ <  < 0.6 ′′ , whereas they were underestimated in the inner regions.The deviations in the phases and amplitudes originated from the biases in the estimation of (cos , PA), which are associated with  = 2 feature.
Finally, Fig. 11 shows the residual images obtained with either real or imaginary parts of images.As evident, the imaging with only the real part successfully extracted the even-symmetric spiral, whereas that with only the imaginary part extracted the odd-symmetric spiral.This is consistent with the discussion in Sec 3.3.Notably, we successfully extracted odd-and even-symmetric spirals for the case of the combined image.Thus, in a realistic situation, where oddand even-symmetric components are mixed, the imaging with either the real or imaginary part is indeed useful for separating them.In addition, the noise level decreased by a factor of ∼ √ 2, or the signalto-noise ratio increased by a factor of ∼ √ 2 because of the assumed symmetry  (, ) =  (−, −) or  (, ) = − (−, −).Thus, this method would be helpful in the search for faint signals.

Axisymmetric structure + circumplanetary disc emission
We injected a point source to the simulated data assuming circumplanetary disc emission and tested the recovering ability.The flux for the point source was assumed to be 0.1 mJy, which is approximately 10 times more significant than the noise level but 10 times lower than the ambient emission at the inner disc.The planetary location was set to be at 0.5 ′′ in the deprojected frame ( ′′ ,  ′′ ) = (0.5 ′′ , 0), such that the source resided in a large gap of a disc.We considered the  same brightness profile and geometric configuration as the simulated case for AS209 in Sec 4 and computed the visibilities and added noises to the data assuming the same observational setup.
Using our method, we similarly recovered the brightness profile and geometric parameters from the simulated data.With the MAP estimate of geometric parameters, we randomly drew one brightness profile and subtracted the model from the simulated visibilities.Fig. 12 shows the residual images, which are generated from the residual visibilities.The position of the injected point source was reasonably recovered, and the flux density was well recovered as well.
Moreover, the central position of the disc was slightly biased, suggesting that the point source might introduce another bias.To check this possibility, we injected the brighter sources with fluxes of 1 and 2 mJy, and attempted to estimate the central positions of the discs by repeating the same analyses.Resultantly, the central positions were estimated to be (Δ cen , Δ cen ) = (0.58 +0.07  −0.07 mas, −0.50 +0.07 −0.07 mas) for the 1 mJy source and (Δ cen , Δ cen ) = (0.97 +0.07  −0.07 mas, −0.82 +0.07 −0.07 mas) for the 2 mJy source.However, the positions of the flux centre with respect to the disc centre were calculated as (0.9mas, −0.9mas) and (1.8mas, −1.8mas) for the 1 and 2 mJy sources, respectively, and they were comparable to the biases in the estimation.Thus, even with our method, the central positions of the disc can be biased by approximately half of the positional difference between the disc and flux centres.This holds true not only for the point source but also for localized emission; for example, a crescent-shaped emission.In the data analysis, in case of bright localized emission, it is recommended that such an additional emission be modelled or removed because it affects the residual image creation process.Specifically, in the case of the point source, we can directly include the point source model in the axisymmetric model.However, if the emission is more complicated, we can model them on the image plane and remove their visibilities from the observations, as presented in Andrews et al. (2021).

APPLICATION TO REAL DATA: ELIAS 20 AND AS 209
We applied the proposed method to the real DSHARP data of two discs, Elias 20 and AS 209, to demonstrate its feasibility.AS 209 is one of most structured discs in the DSHARP sample, and Elias 20 shows a non-axisymmetric feature in the residual map as shown later.More systematic studies on the DSHARP discs will be presented in future studies.
We downloaded the measurement sets of the DSHARP data for the two discs, and applied time and spectral averaging in the same manner as that of Sec 4.1.We manually reduced the recorded weights by 3.50 and 3.44 for Elias 20 and AS 209, respectively, to match the observations.The data were then binned with a log grid with  bin = 500, and they were analyzed with the current method.After sampling the posterior distribution for the parameters, we created the residual images following the method presented in Sec 3.1.1.To create residual images, we used two different geometries: one is obtained from our method, and the other obtained from Huang et al. (2018a), who estimated the geometric parameters through ellipse fitting of annular substructures on image planes.
Table 1 presents the derived geometries for Elias 20 and AS 209.For a reference, we also show the geometric parameters from Huang et al. (2018a).The length scale parameters  for Elias 20 and AS 209 are comparable to their beam sizes (0.076 ′′ , 0.040 ′′ ) and (0.048 ′′ , 0.028 ′′ ), respectively.This result is consistent with the discussion in Sec 4.1 and Appendix D;  was determined by the combination of the intrinsic brightness profile and the UV-coverage.
The derived geometries from our method were mostly consistent with the previous estimates from Huang et al. (2018a), although there are slight or moderate discrepancies.The central positional estimates for AS209 are offset by approximately 1 mas for both directions.Similarly, in the case of Elias 20, the positions also show a difference of 1-2 mas, and notably, there is a 0.07 difference in cos .
These small discrepancies are important for creating the residual images.Figs. 13 and 14 show the brightness profiles, visibilities, and residual images.For comparison, we also show the brightness profiles and model visibilities obtained by Jennings et al. (2022a), who systematically analysed the DSHARP data.Here, they used frank (Jennings et al. 2020), which reconstructs the radial brightness profile by fitting the real part of the deprojected visibilities.Note that Jennings et al. (2022a) adopted the non-negativity condition on the brightness profile in case of AS 209, and assumed additional pointsource emission with flux of 0.66 mJy at the disc centre to suppress the artificial oscillating features for Elias 20 (detailed discussion on point-source correction is presented in Appendix A in their paper).Their model visibilities, as shown in Fig. 13, for Elias 20 are thus converged to 0.66 mJy, which corresponds to the flux of this pointsource emission.
The brightness profiles and model visibilities of the proposed method and that of frank were mostly consistent; however, there existed slight or moderate differences.
In the case of Elias 20, the model visibilities are horizontably offset, and this is due to the difference in the cos .In addition, the models using the proposed method gradually converged to zero at high spatial frequencies, whereas those from frank sharply converged to the flux of the point source around 5 M.The locations of the tipping points for the convergence were determined by  in our method or hyperparameters in frank, and the differences in the way of convergence were due to the different choices for the regularization.
The notable difference was observed in the peak brightness values near  = 0 ′′ .The peak brightness is generally hard to estimate accurately because of the small flux at the small radii.This is mainly due to the point-source correction with 0.66 mJy adopted in Jennings et al. (2022a), which was however not included in their brightness profile.The innermost brightness  ( 1 ) in our model corresponding to the flux 0.66 mJy is ∼ 6 × 10 10 Jy sr −1 , which can largely explain the observed difference of about 10 × 10 10 Jy sr −1 between this study and Jennings et al. (2022a).Another potential explanation could be the difference in the strength of the regularization for smoothing, although identifying the superior model remains challenging at this stage.
In the case of AS 209, the brightness profiles and model visibilities were mostly consistent for our method and that of frank.The peak brightness values, on the other hand, showed the difference of about 2.5 × 10 10 Jy sr −1 , which is yet smaller than that of Elias 20.The discrepancy potentially arises from the difference in regularization strengths or the non-negative condition adopted in Jennings et al. (2020).Overall, our method exhibited the same capability at recovering the brightness profiles as that of frank.
The estimation of the residual images was also improved using our updated geometry.In case of Elias 20, the residual image derived with the geometry from Huang et al. (2018a) exhibits the  = 2 pattern, which is mainly attributed to the inclination offset, as shown in Figs. 4 and 5.This coherent pattern, also seen in Jennings et al. (2020), mostly disappears in our update image, suggesting that our method suppresses the artificial pattern owing to wrong geometric parameters.In case of AS 209, the previous literature identified significant residual patterns (Guzmán et al. 2018;Jennings et al. 2020).The residual image with the updated geometry is less noisy than the previous estimate; however, structured patterns were still present.We thus conclude that the residual pattern for AS 209 cannot be removed by solely optimizing the geometric for a geometrically-thin axisymmetric model.
There can be multiple possible explanations for the residual pattern in AS 209.It may be due to the real non-axisymmetry in the physical parameters, or because of the geometric effect.The individual rings may have different geometric parameters, that is, misalignment or positional offsets, where the latter case was reported for HL Tau (ALMA Partnership et al. 2015).Further, the vertical structure may also render the residual pattern complicated, although the observed residual image is not perfectly consistent with this hypothesis, which predicts the symmetric pattern with respect to the  axis (minor axis).To unveil the origin of the residual pattern for AS 209, we require a more complicated model with multiple geometric parameters or vertical structures; however this is beyond the scope of this paper.

Analysis with axisymmetric disc model
As a practical application to recovering a circumplanetary emission, we applied our method to the data of PDS 70 in Benisty et al. (2021).The same data were already analyzed by frank in Benisty et al. (2021).We used the combined dataset, including long, medium, and short baseline data, for continuum emission in Band 7 as used in Benisty et al. (2021).The data were then averaged in a similar manner to Sec 4 to reduce the data size.In the continuum emission, there is a notable crescent feature in the North-West direction.citebenisty2021 removed the asymmetric feature by following the method described in Andrews et al. (2021).Specifically, using the CLEAN model image, they defined the asymmetry model by isolating the emission of the crescent feature, and subtracted the mean radial profile outside the area from the model, leaving only the asymmetric contribution.The constructed model was then Fourier transformed, and the model visibilities are subtracted from the observed visibilities, which were analyzed using frank.As demonstrated in Andrews et al. (2021), the method is effective for extracting strong asymmetric features that hinder the detection of weak signals, such as emission from CPD (see Appendix B of their paper for details).
In this paper, we applied our method to the original data without any such prior subtraction, aiming to minimize the manual adjustment.A major concern with this simpler approach was that the localized crescent feature might bias the estimates of the geometric parameters, especially the position estimate (Δ cen , Δ cen ), similar to that in Sec.5.2.However, as will be shown later, the derived parameters agree well with those estimated by Benisty et al. (2021), with slight differences, such as 1 mas in Δ cen and 1 deg in PA.The effect of the positional difference in the residual image is negligible in the current analysis (see Fig. F1 in Appendix F).Therefore, we simply present the result of the analysis with our method, without implementing the prior subtraction of the asymmetric features.
We employed a logarithmically spaced grid with  bin = 500 and ( min ,  max ) = (10 2 , 1.1 × 10 7 ) to bin the data.We drew samples from the posterior of parameters using emcee with 16 walkers and 10,000 samples, and we ensured the convergence of MCMC.Table 1 lists the parameters from this study and Benisty et al. ( 2021) (specifi-cally, Appendix B in their paper).The two studies yielded consistent values, while there are slight differences: approximately 1 mas in Δ, 0.007 in cos , and 1 • in PA.
We constructed a visibility model using the MAP solution of parameters, and subtracted it from the observed visibilities.Subsequently, the residuals were processed with CLEAN, with a threshold of nsigma= 3.5.Note that we did not adopt JvM correction adopted in the previous literature (Czekala et al. 2021;Benisty et al. 2021;Balsalobre-Ruza et al. 2023) to prevent potential exaggeration of signal-to-noise ratios in an image (Casassus & Cárcamo 2022).In the imaging process, we experimented with various robust parameters, specifically setting them to be 0.0, 0.5, 1.0, and 1.5.We found that a setting of 1.0 offers a small standard deviation image while maintaining relatively high angular resolution.For comparison, we also generated a residual image using the geometric parameters in Benisty et al. ( 2021) while adopting the hyperparameters in the MAP solution from our modeling.
Figure 15 illustrates the brightness profiles and visibilities derived from our method.For the brightness profiles, we show 30 random samples, indicated by the light orange lines.The brightness profile was characterized by the presence of the inner disc as well as the outer disc with two local maxima, consistently with previous studies (Keppler et al. 2019;Benisty et al. 2021).At the outer disc, there is another shoulder at 0.85 ′′ , as also reported in Benisty et al. (2021).A slight positive excess was also observed in the cavity region, between the inner and outer discs.The flux would arise from emission around PDS 70 b&c and their Lagrange points (Benisty et al. 2021;Balsalobre-Ruza et al. 2023).
Figure 15 2022a) is also shown for comparison.(lower panels) Residual image made with the geometry from our method in the observational frame, that in the deprojected frame, and the image made with that of Huang et al. (2018a) in the deprojected frame.In the residual images, the reference lines at  = 0.5 ′′ are also shown.As treference, the ellipse and circles are shown.The synthesized beam sizes before and after deprojection are (0.048 ′′ , 0.028 ′′ ) and (0.079 ′′ , 0.029 ′′ ) are shown at the bottom left, respectively.the observational frame, the middle one is shown in the deprojected frame, and the right one is based on the geometry in Benisty et al. (2021).The residual images from both geometries are consistent, but there is a slight discrepancy due to difference in geometry.We investigated which parameter made the large difference, and found that the difference in PA largely accounts for the discrepancy (also see Appendix F).
The upper panels in Figure 16 shows the residual images in deprojected coordinate and the polar coordinate.Our methods successfully recovered the circumplanetary emission around PDS 70 c, in addition to the crescent feature.The emission of the circumplanetary disc remains unresolved, consistently with the analysis in Benisty et al. (2021).We measured the brightness of the circumplanetary emission around PDS 70 c.The estimated peak brightness of the emission were 95.2 ± 17.7 Jy beam −1 , 89.8 ± 13.1 Jy beam −1 , and 91.2 ± 11.9 Jye beam −1 for robust = 0, 0.5, 1.0, respectively.Our estimates for robust = 0 and 0.5 were consistent with the result from Benisty et al. (2021), which reported 86 ± 10 Jy beam −1 and 79 ± 7 Jy beam −1 .It should be noted that in our study, the minor/major beam sizes exhibit a slight deviation from those reported in Benisty et al. (2021) (see Table 4) with discrepancies ranging between 0.01 and 0.03 ′′ .This minor variation is likely due to the time and spectral averaging processes applied to our measurement sets, but the discrepancies give negligible impacts on the measurements of intensities for PDS 70 c.On the other hand, our estimate for robust = 1.0 significantly differed from 170 ± 4 Jy beam −1 reported in Benisty et al. (2021).The discrepancy likely arises from the contamination from the main disc emission, which is however mostly mitigated by our method.Ruza et al. 2023).We found the signal to be visually insignificant in the images.Indeed, the peak emission in this region showed significance of 2.0-2.7 for robust = (0, 0.5, 1), slightly lower than the significance around 3.3-3.4reported in Balsalobre-Ruza et al. (2023) for the robust parameters that we considered.The slight discrepancy might be due to the subtraction of the annular brightness in our study.On the other hand, Balsalobre-Ruza et al. ( 2023) observed higher significance at 5-6 if they employed robust ≥ 1.5 and JvM correction, which were not investigated in this study.The further observations and analyses of these marginal signals will be needed.

Additional analysis on residual images by subtracting crescent model
Apart from the crescent feature and circumplanetary emission, residual features still persist in the image.Benisty et al. ( 2021) also iden-tified the residual features in the image after subtracting the axisymmetric component and asymmetric feature including the crescent feature (rightmost end of Fig. 8 in their paper), although the detailed discussion was not presented.Here, we revisited the residual feature in more detail using our updated residual image.
For the further investigation, the high contrast of the bright crescent feature hinders the search of faint features in the images.In addition, the bright crescent unnecessarily introduces the negative flux, because the residual image is supposed to have zero flux when averaged over the polar direction.To avoid the latter problem, Benisty et al. ( 2021) removed the asymmetric feature in the image-based analysis before the visibility analysis.To mitigate these effects, we instead attempted to subtract the bright crescent feature by employing a provisional analytical model.Specifically, we considered a model comprising of the super Gaussian function in the radial direction and the von Mises distribution in the polar direction as follows: where  ,modfied bessel is the modified Bessel function of the first kind,  crescent,amp determines the amplitude of the crescent model,  is the exponent for the super Gaussiaan function,  0 denotes the radial position for the crescent peak,   indicates the radial width,  specifies the azimuthal concentration, and  determines the azimuthal location for the peak.The integral of the equation ( 48) along  is zero, consistent with the construction of the residual image.The adopted model does not perfectly explain the crescent feature, but it is still useful for removing the bright component of the crescent.We fitted the model to the residual image in Sec 7.1, and the optimization was performed us-ing curve_fit in scipy, yielding ( crescent,amp , ,  0 ,   , , ) = (0.339 mJy/beam, 0.793, 0.623 ′′ , 0.0723 ′′ , 5.88, −0.960 rad).
The lower panels in Figure 16 shows the residual images with the crescent models being subtracted in deprojected coordinate and the polar coordinate.We confirmed the presence of emission from PDS 70 c, while no significant detection for point-source emission was observed from PDS 70 b and its L5 point.
The subtraction process facilitated the identification of faint features in the residual images.We here briefly summarize their characteristics in the following: (a) Crescent is trailing : As evident in the residual images with and without the subtraction in a polar coordinate, the observed crescent feature appears to exhibit a trailing pattern rather than circular shape.
Here the disc rotation is clockwise, as we deduced from the velocity gradient map and the emission surface of the gas (Keppler et al. 2019).One interesting possibility is that a planetary gravity perturbs the crescent, resulting in the trailing shape.
(b) Extended emission from crescent to PDS 70 c?: We also observed that the flux excess at the crescent appears to extend toward the vicinity of PDS 70 c.One interesting possibility is that this feature implies a dust accretion flow to PDS 70 c from the outer disc.
(c) Faint arm stemming from inner end of crescent: A faint leading arm is observed to stem from the inner end of the crescent.This feature is also faintly discernible in the residual images without the subtraction.
(d) Arm/Vortex-like structure around 0.75 ′′ stemming from crescent: Another arc or vortex stemming from the crescent was observed, and it is radially concentrated around  = 0.75 ′′ .This feature is more pronounced in our updated residual image than the residual image based on the geometry from Benisty et al. ( 2021), due to the difference in the adopted geometries, especially in PA (also see Appendix F).The base of the arm is overlapped with the crescent model in the lower panels, potentially resulting in over-subtraction of the flux of the arm near the region.
(e) Depletion near crescent: We observed a strong negative feature near the crescent feature.It is visible on both residual maps, with and without subtraction of crescent models.This might be attributed to the counter effect of the dust concentration at the crescent.
(f) Blobs in crescent: We observed blobs inside the crescent, suggesting that the crescent is not unimodal but rather multimodal.There are at least two visible blobs in the crescent feature, located at ( ′′ ,  ′′ ) = (0.53 ′′ , −0.24 ′′ ) and ( ′′ ,  ′′ ) = (0.58 ′′ , 0.02 ′′ ) in the deprojected frames.The blobs exhibited an excess of 3-5 with respect to the background level of the the crescent feature.
All of the features appear to be associated with the crescent feature.In Appendix F, we presented the difference in the residual images derived with geometries from our study and Benisty et al. (2021).Apart from the possible arm (d), all of the features were similarly identified in both of the residual images, reinforcing the coherence of the signals within the current data set.
Substructures in the PDS 70 disc have also been reported in different bands.The arm-like structure was found in the northwest direction in the near-infrared bands, near the crescent observed in the radio band (Müller et al. 2018;Juillard et al. 2022;Christiaens et al. 2024).Juillard et al. (2022) investigated the motion of the arm over six years.However, they did not find the anticipated rotation expected if the arm were comoving with PDS 70 c.This suggests that the arm might be a vortex rather than a spiral.Indeed, as demonstrated in Fig 4 of Marr & Dong (2022), such a circular arc/vortex may be misunderstood as one-armed spiral in the near-infrared band, where the disc thickness becomes significant.On the other hand, the crescent observed in the radio band in this study displays a trailing pattern, which is unlikely to result from a purely geometric effect.One alternative explanation for the arm might be the presence of an undetected companion in the outer disc.Christiaens et al. (2024) set an upper limit on the mass of such a potential planet, obtaining limits 1-4 Jupiter masses from the JWST observation in the near-infrared band.
Additionally, Christiaens et al. (2024) identified a spiral accretion stream to the vicinity of PDS 70 c in the near infrared band by the JWST observation.It might be related to the extended emission excess near PDS 70 in the radio band (feature b in this study), but the connection between the features is not clear at this point.
The present findings rely on the assumption that the disc structure is well approximated by the thin axisymmetric disc with a single set of geometric parameters.This underlying assumption may be inadequate if the inner and outer discs are misaligned, or the disc thickness is substantial.In addition, there is still room for improvement for the modeling of the crescent feature, and the model might be better to be incorporated into the framework of the axisymmetric modelling.A comprehensive analysis of these points is essential for the robust detection the features, and we leave this issues for the future study.

SUMMARY
This study proposed a scheme for estimating the geometry, hyperparameters, central position, and brightness profiles assuming a geometrically thin disc in radio interferometry.Our approach is less susceptible to human biases due to manual tuning of parameters in contrast to frank, where the non-linear parameters need to be determined a priori.
Simulating observations for an axisymmetric disc, we demonstrated that the proposed method can successfully retrieve geometric parameters in a more precise manner than that using Gaussian fitting.Additionally, we performed injection and recovery tests for the nonaxisymmetric perturbations to the simulated data, and showed that the proposed method can reasonably recover the injected structures.However, the estimated geometric parameters were slightly shifted from the assumed values.This is attributed to the degeneracy between the non-axisymmetric perturbations and the residuals caused by biases in the geometric parameters.
The model was then applied to the real data for Elias 20 and AS 209, and the ability of the method to determine the disc geometry and brightness profile was demonstrated especially for Elias 20.Moreover, the data for the continuum emission of the PDS 70 were reanalyzed with our method.Our methods successfully identified the circumplanetary emission from PDS 70 c and the crescent feature.Furthermore, we tentatively identified several new features, including trailing nature of crescent, extended emission near PDS 70 c, armlike structures, dust depletion near crescent, and blobs.The origins of these features are unclear, and a further modeling is needed.
The current methodology can be applied to any type of continuum data in radio interferometric observations, and future studies will analyse the DSHARP data (Andrews et al. 2018) to explore their non-axisymmetric structures.
One of the notable strengths of this research is its capacity to adapt to more complex problems involving a multitude of geometric or hyper parameters, which are difficult to adjust manually.Below, we present the potential directions for extending our study: • The model can be extended to cover multiple rings/gaps with different central positions, inclinations, or position angles as observed in GW Ori (Bi et al. 2020).This can be realized by considering multiple sets of geometric parameters, each of which is applied to one of the separate ranges.
• The current model can be extended to incorporate frequencydependent radial profiles for multi-frequency data.We can straightforwardly expand our formulation to include the frequency dependence in a linearized formulation by Taylor expansions in a manner similar to multi-frequency CLEAN (Rau & Cornwell 2011): where  (; ) is the radial brightness at the frequency ,  0 is the reference frequency, and   , 0 () is the coefficient of -th Taylor expansions.With this expression, we can analytically derive the model parameters in the same manner as that presented in the current study.
• The current study focuses on the axisymmetric component  (); however, we may be able to include non-axisymmetric components (  (),   ()) in the model, and directly solve them.Such modelling efforts might ease the degeneracy between the geometric parameters and non-axisymmetric structures.
• Although we created residual images with CLEAN, we can also use other imaging techniques, for example, sparse modelling, to produce images alternatively (Honma et al. 2014;Nakazato & Ikeda 2020;Yamaguchi et al. 2020;Aizawa et al. 2020).In particular, sparse modelling favours solutions with many zeros, and it can achieve better angular resolutions; thus it will be helpful for resolving or identifying the new substructures, including spirals and circumplanetary emission.
• A more comprehensive analyse for PDS 70 will be essential for the reliable detection of the discovered structures.Incorporating the crescent model into our model will allow us to separate the residual signal accurately.Exploring more complicated models that consider a disc thickness or misalignment between inner and outer disc planes will be rewarding.The multi-wavelength data will be also useful for assessing the consistency of the signals across different wavelengths.
These will be considered in future studies. H Σ| H Σ−1   + , (A4) where  corresponds to the terms independent on ( , ).In the case of  ≫ , this incurs a computational cost that is much lesser than that using equation (32).

APPENDIX B: REVISITING VISIBILITY BINNING WITH UNIFORM AND LOG GRIDDING
The computational cost strongly depends on the number of visibilities.Here, we discuss the data binning in an interferometric observation assuming linear and log grids.

B1 Formulation
We consider the data  obs with the data weights  obs .We binned them on a grid specified by bin edges {( i ,  j )} with  = −( bin + 1), − bin , ...,  bin ,  bin + 1, wherein  bin determines the number of grids.For a cell with  −1 <   <   and   −1 <   <   , we define the summing operation for a vector  with the same length as  as follows: In each cell, we computed a weighted average  bin , a total sum of weights  bin , standard deviations for the noise  bin,,  , and weighted average of spatial frequencies ( bin ,  bin ) as follows:  (125,250,500,1000,2000,4000),  min = 10 2 , and  max = 10 7 .At the binned spatial frequencies ( bin ,  bin ), we computed model visibilities  bin,model , and quantified a binning error relative to an observational error as follows: where  d,bin is the number of data after binning.We computed the binning errors for the linear and log grids by adopting  bin = (125, 250, 500, 1000, 2000, 4000).The left panel in Fig. B1 shows the relation between the mean squared binning errors  2 bin,mean and the number of binned data  d,bin , which determines the computational time.As  bin increases, the binning errors became small, and  d,bin became large as expected, thereby increasing the computational cost.On average, uniform and log grids yielded similar binning errors with  d,bin being fixed as shown in Fig. B1.
The right panel in Fig. B1 shows the binning errors |  bin | at deprojected spatial frequencies .To obtain a similar number of data points  d,bin ≃ 10 5 , we assumed  bin = 2000 for the uniform grid and  bin = 500 for the log grid (see left panel in Fig. B1).The binning errors in the uniform grid increased at low spatial frequencies, and they reduced at higher frequencies.With  bin = 2000 for the uniform grid, the binning error can reach 40 per cent at most, and such errors might be problematic to estimations on large-scale emissions.The binning errors increase at small spatial frequencies because the binning errors in the uniform grid are proportional to the first derivative of the radial visibility profile, which tends to be large at smaller scales.This is supported by the small binning errors at large spatial frequencies.In the case of the log grid, the binning errors are suppressed at small spatial frequencies, in contrast to the uniform grid.This accords with the narrow bin widths at small spatial frequencies in the log grid; the grid is sufficiently fine to resolve the visibility profile.At the middle to high spatial frequencies, the binning errors increased and decreased, with a peak at 10 6 M.At the higher spatial frequencies, the grids become coarse, while the first derivative of the radial visibility profile becomes small.These two different trends result in this dependency.
Comparing the two grids, the log grid tends to yield moderate binning errors, whereas the uniform grid yields errors in a broad range.Considering the robustness, we adopted the log grid in the analysis of the main text, because the large binning errors at small spatial frequencies for the uniform grid can degrade the estimated accuracy of the flux on a large scale.

B3 Dependence of geometric parameters and hyperparameters on number of grids
In the main text, we adopted the log grid with  bin = 500.Here, we show that the choice of  bin does not affect the estimation on non-linear parameters.Simulating the same observational setup and the brightness profile as in the simulated case for AS 209 in Sec 4, we derived the geometric parameters by varying the number of grids,  bin = 125, 250, 500, and 1000.We injected (Δ cen , Δ cen , cos , PA) = (0 ′′ , 0 ′′ , 0.75, 45 • ).Fig. B2 shows the results for recovered geometric parameters and hyperparameters.The results appear to be unaffected by the choice of  bin .

Figure 1 .
Figure1.Schematic of a disc in observed (a) and deprojected frames (b).The disc is inclined with cos  and oriented by PA with respect to the -axis in the projected frame.The origin of ( , ) is set to be the phase centre, and (  ′ ,  ′ ) is shifted from ( , ) by (Δ cen , Δ cen ).(  ′′ ,  ′′ ) is the coordinate for the deprojected frame, and the brightness profile is assumed to be only radially dependent on  ′′ ( ′′ ).

Figure 2 .
Figure 2. Recovery and injection test for mock data of AS 209.(upper-left) Posterior distribution for (Δ cen , Δ cen , cos , PA, , ).Dotted lines indicate the input values.(upper-right) Injected brightness profile, denoted by the black line, and 10 samples drawn from the posterior distribution of brightness profiles, denoted by red lines.At the bottom, the residuals between injected and recovered models are shown.(lower-left) Simulated visibilities denoted by blue points and model visibilities denoted by an orange line.The visibilities are binned with a logarithmic grid with  = 2000.In the lower panels, the residuals indicating the difference between the simulated and model visibilities are shown.(lower-right) Residual image produced with the MAP estimate in the observed frame.The synthesized beam size (0.076 ′′ , 0.040 ′′ ) is shown at the bottom left.

Figure 3 .
Figure 3. Recovery and injection test for mock data of WaOph 6 in Sec 4. The format of the figure is same as that of Figure 2. The synthesized beam size (0.091 ′′ , 0.037 ′′ ) for the residual image is shown at the bottom left.

Figure 4 .Figure 5 .
Figure 4. Unsubtracted image and residual images produced with incorrect geometric parameters in simulated case for AS 209.The standard deviation for the brightness  is computed in the region outside the disc.(upper-left panel) The image produced with CLEAN for unsubtracted data.(upper-middle panel) The residual image produced with residual visibilities produced from the MAP estimate.(other panels) The residual images with each of the geometric parameters being shifted.The synthesized beam size (0.076 ′′ , 0.040 ′′ ) is shown at the bottom left.

Figure 6 .Figure 7 .Figure 8 .
Figure 6. =1,2 ( ) computed for residual images in Fig. 5.The left and right panels show the results for the shifted central position and the shifted disc orientation, respectively.

Figure 9 .
Figure 9. Injected and recovered residual images for spirals in deprojected frames; (top rows) Even-symmetric spiral (middle rows) Odd-symmetric spiral (bottom rows) Combined spirals.(left columns) Injected perturbations (middle columns) Recovered residual images based on the proposed method (right columns) Recovered residual images based on the geometry obtained with Gaussian fitting.The synthesized beam size for the deprojected image (0.11 ′′ , 0.041 ′′ ) is shown at the bottom left, and the beam size for the image before deprojection is (0.091 ′′ , 0.037 ′′ ).

Figure 10 .
Figure 10.Phases and amplitudes for injected and recovered odd-(left) and even-(right) symmetric spirals in Sec 5. (upper panels) The residual images in deprojected frame overplotted with the locations of injected (black) and recovered phases (blue).(middle panels) Input and recovered phases in the radial direction.(lower panels) Input and recovered amplitudes of spirals.

Figure 11 .
Figure 11.Demonstration of extraction of odd-and even-symmetric spirals with real and imaginary parts of the visibilities; (top rows) Residual images for even-symmetric spiral in deprojected frames (middle rows) Residual images for odd-symmetric spiral (bottom rows) Residual images for combined spirals (left columns) Injected spirals (middle column) Image produced with only the real part of the data (right column) Image produced with only the imaginary part of the data.

Figure 12 .
Figure 12.Residual images for the simulated data that include a simulated circumplanetary disc emission, whose position is indicated by a red dashed circle.The central and right panels show the maps in an observed and deprojected frames, respectively.The synthesized beam sizes before and after deprojection are (0.076 ′′ , 0.040 ′′ ) and (0.089 ′′ , 0.045 ′′ ) are shown at the bottom left, respectively.

Figure 13 .
Figure 13.Analysis of DSHARP data of Elias 20. (upper-left panels) Brightness profiles recovered by this study and Jennings et al. (2022a).The profiles are shown in linear and log scales in the upper and lower panels, respectively.(upper-right panels) Model and observed real visibilities and deprojected spatial frequencies.The model by Jennings et al. (2022a) is also shown for comparison.(lower panels) Residual image made with the geometry from our method in the observational frame, that in the deprojected frame, and the image made with that ofHuang et al. (2018a) in the deprojected frame.In the residual images, the reference lines at  = 0.5 ′′ are also shown.As treference, the ellipse and circles are shown.The synthesized beam sizes before and after deprojection are (0.048 ′′ , 0.028 ′′ ) and (0.079 ′′ , 0.029 ′′ ) are shown at the bottom left, respectively.

Figure 14 .
Figure 14.Analysis of DSHARP data of AS 209.The format of the figure is the same as that of Fig. 13.The reference lines corresponding to gaps and rings are shown in brightness profiles and the residual images.The synthesized beam sizes before and after deprojection are (0.076 ′′ , 0.040 ′′ ) and (0.076 ′′ , 0.049 ′′ ), respectively.

Figure 15 .
Figure 15.Analysis for PDS 70.The format of the figure is the same as that of Fig. 13 except that the zoom-up view of the radial profile is shown in the linear scale.The reference lines are shown at  = (0.4 ′′ , 0.6 ′′ , 0.8 ′′ ) in the panels of the brightness profile and the residual maps.The synthesized beam sizes before and after deprojection are (0.063 ′′ , 0.053 ′′ ) and (0.098 ′′ , 0.053 ′′ ), respectively.

Figure 16 .
Figure 16.The residual images for PDS 70 in the deprojected coordinate (left panels) and in the polar coordinate (right panels).The upper and lower panels show the residual images without and with the subtraction of the crescent models, as discussed in Sec 7.2.Contours of the crescent models are also overlaid as black lines with levels of (50, 100, 200) Jy beam −1 .Annotations in the upper panels show the locations of the protoplanets and the L 5 point of PDS 70 b, as provided by Balsalobre-Ruza et al. (2023).Annotations in the lower panels specify the possible structures identified in this system, and the further explanations for the features (a)-(f) are provided in Sec 7.2.

Table 1 .
Hyperparameters for Gaussian Process and geometric parameters of Elias 20, AS 209, and PDS 70.Geometric parameters assumed in previous studies are also shown for comparison.