Observability of Substructures in the Planet-forming Disk in the (Sub)centimeter Wavelength with SKA and ngVLA

Current imaging observations of protoplanetary disks using the Atacama Large Millimeter/submillimeter Array (ALMA) primarily focus on the submillimeter wavelength, leaving a gap in effective observational approaches for centimeter-sized dust, which is crucial to the issue of planet formation. The forthcoming Square Kilometre Array (SKA) and ngVLA may rectify this deficiency. In this paper, we employ multifluid hydrodynamic numerical simulations and radiative transfer calculations to investigate the potential of SKA1-Mid, ngVLA, and SKA2 for imaging protoplanetary disks at subcentimeter/centimeter wavelengths. We create mock images with ALMA/SKA/ngVLA at multiwavelengths based on the hydrodynamical simulation output and test different sensitivity and spatial resolutions. We discover that both SKA and ngVLA will serve as excellent supplements to the existing observational range of ALMA, and their high resolution enables them to image substructures in the disk’s inner region (∼5 au from the stellar). Our results indicate that SKA and ngVLA can be utilized for more extended monitoring programs in the centimeter wave band. While in the subcentimeter range, ngVLA possesses the capability to produce high-fidelity images within shorter observation times (∼1 hr on source time) than previous research, holding potential for future survey observations. We also discuss for the first time the potential of SKA2 for observing protoplanetary disks at a 0.7 cm wavelength.


Introduction
Planet formation and its intricate relation with protoplanetary disks have remained at the forefront of astrophysical research for decades.One of the most pivotal pieces of this puzzle lies in the substructures of these disks, especially in the form of gaps and rings.The Atacama Large Millimeter/submillimeter Array (ALMA) has been invaluable in observing these substructures in many systems (e.g., ALMA Partnership et al. 2015;Fedele et al. 2017;Andrews et al. 2018;Clarke et al. 2018;Dipierro et al. 2018;Long et al. 2019Long et al. , 2020Long et al. , 2022Long et al. , 2023)), thanks to its high-resolution capabilities.
These substructures can provide crucial information about the dynamic interplay between forming planets and their natal environments.While several mechanisms, such as changes in dust properties around ice lines (Zhang et al. 2015;Okuzumi et al. 2016) or inherent (magneto)hydrodynamic (MHD) instabilities in disk (Flock et al. 2015;Suriano et al. 2018;Riols & Lesur 2019;Riols et al. 2020;Elbakyan et al. 2022;Wu et al. 2023b), have been postulated to account for these patterns, a prevailing hypothesis is that they are a direct consequence of interactions between the disk and embedded, still-forming planets (see Kley & Nelson 2012;Baruteau et al. 2014; Andrews 2020 for reviews).
As a young massive planet orbits its host star, it exerts gravitational forces on the surrounding disk material, potentially carving out a gap in its path.Notably, localized continuum excesses (Keppler et al. 2018;Tsukagoshi et al. 2019;Nayakshin et al. 2020) and discernible kinematic deviations in the disk's molecular gas-attributed to the influence of an embedded planet (Teague et al. 2018;Pinte et al. 2019Pinte et al. , 2020Pinte et al. , 2023))-have further underscored this interpretation in recent years.
While there is an abundance of observational evidence pointing to planet-induced substructures in discs, the leap from millimeter-sized dust particles to kilometer-sized planetesimals is still a challenge.Collisions at the millimeter scale often result in fragmentation or a rebound effect, rather than the anticipated merging (see Blum & Wurm 2008 for review).Additionally, the solid material within the disk is subject to aerodynamic drag from surrounding gases, which saps angular momentum, propelling these materials swiftly toward the central star (Weidenschilling 1977).Although some theories suggest accelerated growth mechanisms-like the streaming instability (e.g., Youdin & Goodman 2005;Johansen et al. 2007;Wu et al. 2024b)-it is evident that the transition from millimetersized dust to centimeter-sized pebbles requires further elucidation to unlock the intricacies of planet formation.Notably, at millimeter wavelengths, it is primarily millimeter-sized dust grains that dominate dust opacity (Draine & Lee 1984;Semenov et al. 2003;Draine 2006).Hence, ALMA, spanning a wavelength range from 0.3 to 3 mm, is adept at capturing emissions from roughly millimeter-sized grains.Yet, its sensitivity falters when observing emissions from larger, centimeter-sized dust grains.
In any sense, observations at longer wavelengths are crucial to improve the chance of revolving the high optical depth region, and the only way to robustly decide the size distribution of pebbles in the corresponding wavelength range.Yet, currently, the Janksy Very Large Array (VLA) is the only facility capable of providing reasonable sensitivities and angular resolution at the centimeter wavelength, whose observation could only be archived in a handful of the brightest disks.
Upcoming facilities, like the Square Kilometer Array's midfrequency component (SKA1-Mid; Braun et al. 2015Braun et al. , 2019) ) and its subsequent phase (SKA2), along with the nextgeneration Very Large Array (ngVLA) (Murphy et al. 2018), offer significantly improved sensitivity and resolution expected at (sub)centimeter bands.Only with these, larger-size grain populations and their distribution in protoplanetary disks can be resolved in the planet-forming astronomical unit scales (Ricci et al. 2018;Harter et al. 2020;Ilee et al. 2020;Ricci et al. 2021;Ueda et al. 2022).
In this work, we simulate the potential of SKA1-Mid, SKA2, and ngVLA in observing the substructures, particularly gaps and rings, created by gap-opening planets in protoplanetary disks.By juxtaposing their capabilities with the existing benchmark set by ALMA Band 1, we endeavor to provide a clearer picture of how future observations can shed light on the age-old questions surrounding planet formation, particularly through the lens of dust grain distribution in protoplanetary disks.
Our paper is organized as follows: In Section 2, we introduce the numerical setup for our hydrodynamic simulations, radiative transfer calculations, and parameter choices for observational mock imaging.In Section 3, we described the mock images for ALMA band 1, SKA, and ngVLA at the (sub) centimeter wavelength.We present a summary of our findings and engage in a comprehensive discussion regarding the implications of our results in Section 4.

Numerical Setup
In this section, we introduce the methods that we used to simulate the high-resolution observation in radio wavelengths.

Hydrodynamic Simulations
In this study, we used grid-based code FARGO3D (Benítez-Llambay & Masset 2016) to do 2D {r, f} hydrodynamical planet-disk simulations.This code is the successor of the public FARGO code (Masset 2000); in the latest version, it added a multifluid module (Benítez-Llambay et al. 2019), and we can consider the dynamics of gas and dust in a disk with one or more embedded planet.This code replaced dust-particle size with a constant Stokes Number to describe the interaction with the gas.However, in this study, we used fixed dust sizes to describe the drag force between fluids.
We adopt the same solar-like system numerical simulation settings as the previous works (Ricci et al. 2018).In this study, we adopted the following initial radial profile of gas surface density: this equation is the self-similar solution for the evolution of a viscous accretion disk when the viscosity ν(r) ∝ r γ (Lynden-Bell & Pringle 1974).For Equation (1), r is the radius in the disk, r c is a truncation radius, and normalization constant 2 is a parameter used to describe the initial total gas mass in the disk.In our simulations, we set values of γ = 0.8, r c = 17 au, and the total disk mass M disk = 8 × 10 −3 M e .
For the initial dust surface density, we set where dust-to-gas mass ratio ò = 0.01.We distributed 11 different sizes of dust ranging from 7 μm to 1 cm in our simulations.For each size of dust particles, their mass follows grain-size (a) distribution n(a) ∝ a − q (e.g., Ricci et al. 2010), where q = 3.5.We set the internal density of dust to be 1.2 g cm −3 for all sizes.
The mass of the central star is set M å = 1 M e .We set the geometrical aspect ratio of the disk, H/R, to where h 0 = H/r = 0.05, f = 0.25 is a constant appropriate for the disk midplane temperature profile T(r) ∝ r −1/2 , and r 0 is the planet orbital radius.For the gas viscosity, we utilize the α prescription of Shakura & Sunyaev (1973) and adopt two values of the α parameter, 10 −3 and 10 −5 .
In this study, we focus on the disk structure of massive planets, and we only consider the same orbital radius, We let the orbital radii be kept fixed during each simulation.In our simulations, we put a Jupiter mass planet (planet-to-star mass ratio is q = 1 × 10 −3 ) for gap-opening.We run the simulations for 1500 orbital periods, and the planets' orbital radii are set to 5 au in our simulations; for our orbital radii, this corresponds to 16,700 yr.
Our numerical grids extend from 0.1 r 0 (0.5 au) to 12 r 0 (60 au) in the radial direction and 0 to 2π in the θ direction.We utilized 1536 grids in the radial direction and 1024 grids in the θ direction.

Radiative Transfer Calculations
We utilized the publicly available RADMC-3D code (Dullemond et al. 2012) to conduct our dust radiative transfer calculations.Then, to process our FARGO3D simulations for use as inputs in RADMC-3D and to compute the final beamconvolved synthetic maps of dust continuum emission, we employed the publicly available Python code fargo2r-admc3d (Baruteau et al. 2019).
The configuration of the dust radiative calculations closely follows the methodology described in Baruteau et al. (2021); this method has been applied in multiple recent studies (e.g., Wu et al. 2023aWu et al. , 2023b)).Notably, the dust temperatures employed in RADMC-3D align with the temperatures utilized in the hydrodynamical simulations.Here, we only briefly list the main parameters used in our study.
To extend the 2D surface density in hydrodynamical simulation to 3D spherical grids for radiation transferring calculation, we consider vertical hydrostatic equilibrium and set the dust scale heights to The dust size distribution is described in Section 2.1.
The opacities for dust absorption and scattering are computed using Mie theory, and are produced by OpTool (Dominik et al. 2021); the grain composition is the same as DSHARP composition (Birnstiel et al. 2018;Zhang et al. 2018), which contains the total opacities of absorption and scattering for a given wavelength and grain size.The specific intensity of the dust continuum emission is computed with 10 8 photon packages used for ray tracing.
In Figure 1, we plot the relationship between the total optical depth (absorption + scattering) and the maximum size dust a max at different observation frequencies (wavelengths).This shows a huge potential of the observations in longer λ to constrain the dust size distribution.
We set the disk distance as 150 pc away from antennas, as this distance encompasses hundreds of disks in regions such as Taurus (Loinard et al. 2007), Ophiuchus (Ortiz-León et al. 2017), Chamaeleon (Whittet et al. 1997), and Lupus (Lombardi et al. 2008).We assume a stellar effective temperature of 6000 K with a star radius of 1.5 R e .Taking radiative transfer at the 0.7 cm wavelength as an example, our total disk flux is 0.26 mJy (for a disk model with α = 10 −3 ) and 0.2 mJy (for a disk model with α = 10 −5 ).

SIMOBSERVE
For each hydrodynamic simulation, we make synthetic maps to mock ALMA, SKA, or ngVLA images at varied wavelengths.We utilized the SIMOBSERVE task within the Common Astronomy Software Application (CASA 6.4.9) package to transform synthetic images into visibility data sets in the Fourier domain.
For configuring the ALMA and ngVLA observations, we employed the alma.out28.cfgand ngvla-main-revC.cfgfiles provided within the CASA package.The ALMA configuration includes baselines extending up to 16 km, while the ngVLA's reference design features 214 antennas, each 18 m in diameter, with baselines reaching approximately 1000 km.For the ngVLA, we utilized SETNOISE for noise processing.And for TCLEAN task, we set psfcutoff to 0.75 because our pointspread functions have a large skirt.
For SKA1-Mid, we utilized the beam size and weighted continuum sensitivity calculated by the SKAO Sensitivity Calculator, 5 and then created synthetic maps using fargo2radmc3d.For observations at 12.5 GHz on 1 hr source time, the beam size of SKA1-Mid is 0 025 × 0 022, with a sensitivity of 14.5 μJy beam −1 .Since SKA2 remains conceptual at this stage, we adopted two design values for the beam size (0 01 × 0 01 and 0 019 × 0 019) provided by Braun et al. (2019) and a reference weighted continuum sensitivity of 3.2 μJy beam −1 for 1 hr of observation time to create synthetic maps.
For this study, we conducted our SIMOBSERVE at two different frequencies for SKA.One frequency of 12.5 GHz (wavelength ∼2.4 cm) is positioned at the midpoint of the suggested SKA1-Mid Band 5b frequency range (8.3-15.3GHz).This choice is similar to Ilee et al. (2020).The other one is 43.3 GHz (wavelength ∼0.7 cm), which is a possible upgrade path of SKA2 (as part of the SKA Observatory development plan; Braun et al. 2019).For ngVLA, we took into account three distinct frequency bands, Band 3 (16 GHz, wavelength ∼1.87 cm), Band 4 (30 GHz, wavelength ∼1.0 cm), and Band 5 (43.3 GHz, wavelength ∼0.7 cm).As a control, we also made the synthetic maps of ALMA Band 1 (43.3GHz, wavelength ∼0.7 cm), and specifically, the antennas of ALMA Cycle 11 (C43-10).We summarize the array layout definitions in Table 1 for the reader's convenience.

Results
In this section, we display the results of our multifluid hydrodynamic simulations for a disk with an embedded planet and the interferometric observations' potential with SKA, ngVLA, and ALMA.In Figure 2, we present the sky model for different models that we used for the next SIMOBSERVE.We noticed that the model with α = 10 −3 exhibits a very symmetrical and bright ring with a deep gap.The model with α = 10 −5 , on the other hand, generates a bright vortex in the horseshoe region and asymmetric crescent-like structures in the inner region, as well as wider gaps and narrower rings in the outer region.This is reasonable because its extremely low viscosity will make it easier to open gaps (e.g., Kanagawa et al. 2016;Elbakyan et al. 2022), and make a secondary inner gap (Dong et al. 2017;Huang et al. 2018).And due to Rossby wave instability (Lovelace et al. 1999;Li et al. 2001), it will generate a vortex accordingly.
We regard the successful observation of bright rings, dark gaps, and all asymmetric structures presented in the fluid model as critical factors in assessing the antenna's performance and potential.

Observation at the 2.4 cm Band with SKA1-Mid
Figure 3 shows the results of the SKA1-Mid (a) observations for models with different α viscosity α = 10 −3 for the left panels and α = 10 −5 for the right panels).SKA1-Mid (a) represents the high-resolution configuration design for SKA1-Mid, featuring an angular resolution of 0 025 × 0 022.As shown in Figure 1, at a wavelength of 2.4 cm, the opacities are predominantly attributed to the scattering contributions from grains of sizes 5 mm and 1 cm.
The top panels and the middle panels show the realistic observational outcomes, assuming an observation duration of 1 and 10 hr of source time; here the sensitivity is 14.5 μJy beam −1 (for 1 hr) and 4.6 μJy beam −1 (for 10 hr).Under this condition, utilizing SKA1-Mid (a) for observations, we are virtually unable to discern any structure within the disk.
The bottom row displays the mock images without any noise, i.e., the observation results under ideal conditions.Under the resolution and observational wavelengths of SKA1-Mid (a),  6 SKA2 (plan b) 43.3 0.7 0.019 × 0.019 3.2 Figure 6 ngVLA Band 5 43.3 0.7 0.0027 × 0.0023 0.32 Figure 6 ALMA Band 1 C43-10 (Cycle 11) 43.3 0.7 0.1 × 0.1 8.5 Note.These sensitivity values are specifically applicable to observations conducted with 1 hr of on-source time.ideally, we can resolve most substructures appearing in Figure 2, including the gap at 5 au, the bright ring out of the planet, and the vortex in the low-viscosity model.However, SKA1-Mid is incapable of resolving the secondary inner gap within the low-viscosity model, due to the beam size resulting from the resolution exceeding this area (less than 1 au).

Observation at the 1.87 cm Band with ngVLA
For ngVLA Band 3, we evaluated the observational capabilities at a wavelength of 1.87 cm in Figure 4. Referring to Figure 1, we discern that at this wavelength, opacities predominantly arise from the scattering contributions of grains sized at 2.3 mm, 5 mm, and 1 cm.The overall opacities are akin to those observed at a wavelength of 2.4 cm.
In comparison to SKA1-Mid, the ngVLA boasts superior resolution (0 0051 × 0 0034) and sensitivity (0.23 μJy beam −1 on 1 hr source time).Under long-time observational conditions (such as 10 hr on-source time, as shown in the bottom panels), ngVLA can resolve most substructures within the disk at a wavelength of 1.87 cm, except the secondary inner gap within the low-viscosity model.However, constrained by the low dust emission rates of centimeter/subcentimeter-sized dust in the centimeter band, ngVLA Band 3 is incapable of resolving the disk's substructures with only 1 hr of source time, except for the emission of the host stellar.

Observation at the 1.0 cm Band with ngVLA
In Figure 5, we assume the naturally weighted point-source sensitivity is 0.24 μJy beam −1 on 1 hr source time for ngVLA observations at the 1 cm wavelength, with an angular resolution of 0 0031 × 0 0025, higher than that presented by Ricci et al. (2018).We observe that even with only 1 hr of observation time, ngVLA Band 4 is capable of fully resolving the substructures within the disk at a 1 cm wavelength.Although the signal is quite faint and the signal-to-noise ratio relatively low for the low-viscosity model, it remains significantly above the noise level.At this point, dust emission predominantly originates from grains measuring 1.1 mm, 2.3 mm, 5 mm, and 1 cm.

Observation at the 0.7 cm Band
The designed frequency bands for ngVLA include 43.3 GHz (i.e., λ = 0.7 cm, rms = 0.32 μJy beam −1 on 1 hr source time), which is also among the planned frequency bands for SKA2 (rms = 3.2 μJy beam −1 on 1 hr source time).Coincidentally, the existing ALMA has just acquired the capability to conduct observations in the subcentimeter band, specifically ALMA Band 1, which similarly covers the 0.7 cm wavelength.In Figure 6, we undertake a comparison of the observational capabilities of these three distinct facilities.Within this context, ALMA Band 1 exhibits inferior resolution and sensitivity, standing at 100 mas and rms = 8.5 μJy beam −1 on 1 hr source time, respectively.SKA2 (plan a) boasts a resolution on par with ngVLA, at 0 01 × 0 01.Conversely, SKA2 (plan b) presents a slightly reduced resolution of 0 019 × 0 019, falling short of both its counterparts.
As illustrated in Figure 6, with an observation time of 1 hr, ALMA Band 1 is incapable of resolving any structures at the scale of the disk, a result that is unsurprising.In contrast, SKA2 is capable of vaguely identifying the general outlines of bright rings and gaps within the high-viscosity model.Interestingly, these substructures exhibit a signal-to-noise ratio in SKA2 (plan b) that is, despite its lower resolution, even better than the signal-to-noise ratio in the high-resolution SKA2 (plan a).However, SKA2 falls short in resolving structures within the low-viscosity model.In comparison to these three facilities, ngVLA Band 5 performs exceptionally well, capable of resolving any structures in both models.

Discussion
The multiwavelength SIMOBSERVE from 0.7 to 2.4 cm presented in Section 3 shows that the SKA and ngVLA will bridge the observational band-gap left by ALMA in probing substructures, such as gaps and rings, created by the interaction between protoplanetary disks and forming planets.
We employed a hydrodynamics setup akin to that of Ricci et al. (2018) and expanded the observation frequency band of ngVLA, discerning a similar observational prowess in the ngVLA as previously anticipated.In contrast to Ilee et al. (2020), our approach involved the utilization of a gap-opening planet, thereby testing the SKA's potential in scrutinizing disk  substructures.Furthermore, we conducted preliminary studies for the proposed SKA2, anticipating its future contributions to the field.
Moreover, for substructures induced by planets in close proximity to their host stars, such as those near 5 au, both the SKA and ngVLA possess the resolution capacity to delineate them.The ngVLA, in particular, has the potential to resolve vortices and the asymmetry generated by planet-disk interactions, as well as secondary gaps that emerge in more inner regions under conditions of low viscosity.One of the hottest topics of late-the MHD disk wind-is particularly active in the more inner regions of disks (e.g., Armitage et al. 2013;Bai & Stone 2013, 2014;Bai et al. 2016;Elbakyan et al. 2022;Aoyama & Bai 2023;Wu et al. 2023b).At the same time, compared to broader disk regions, certain effects may have a more significant impact in the inner disk region, such as the Hall effect and strong turbulence, which possess undeniable potential regarding planetesimal formation (Wu et al. 2024b) and planet chaotic migration (Wu et al. 2024a).The resolving capabilities of SKA and ngVLA in the inner disk region are instrumental in understanding these effects.
While enhancements in the signal-to-noise ratio for SKA and ngVLA observations at centimeter wavelengths remain a consideration, subcentimeter wavelength observations hold the promise of high-fidelity imaging within reasonable integration times for point observations with these facilities.Consequently, there is an anticipated expansion in the scope of multiwavelength observations utilizing a broader spectrum range with the SKA and ngVLA in the future.
Recent multiwavelength studies (e.g., Carrasco-González et al. 2019;Macías et al. 2019Macías et al. , 2021) ) employing ALMA and the VLA, as well as rare CO isotopologue surveys (Booth et al. 2019;Booth & Ilee 2020) targeting protoplanetary disks, suggest that disks may indeed be optically thick at submillimeter wavelengths.This condition could lead to misconceptions regarding fundamental characteristics of protoplanetary disks, such as the mass of dust and gas, along with the distribution of dust sizes.By operating at longer wavelengths with the SKA and ngVLA, it might be possible to observe disks that are optically thin, revealing structures that, while present, are less discernible at submillimeter wavelengths due to their concealment beneath optically thick regions (Wu et al. 2023b).The feasibility of this hypothesis has been verified in (Liu 2019), as the spectral index of less than 2.0 in the center of TW Hya can only be explained by optically thick dust.Recently, the morphological comparison of DM Tau at 1 mm and 7 mm wavelengths indicates that DM Tau is extremely optically thick at 1 mm (Liu et al. 2024).

Conclusions
We present the hydrodynamic simulation results from FARGO3D, encompassing disks with embedded planets and the intricate interactions between disk gas and dust.The primary aim of this study is to explore the capabilities of current (ALMA Band 1) and forthcoming (SKA and ngVLA) telescopic facilities in probing planet formation at subcentimeter/centimeter wavelengths.Our main findings are as follows: 1. SKA1-Mid, ngVLA, and SKA2 all demonstrate the capacity to resolve substructures in protoplanetary disks as well as their inner regions at subcentimeter/centimeter wavelengths.2. At centimeter wavelengths, SKA1-Mid and ngVLA Band 3 require extended integration times and high resolution.
For studies of protoplanetary disks, future enhancements should prioritize boosting the signal-to-noise ratio to facilitate research within more reasonable observational timeframes.Alternatively, emphasis could be placed on more extended monitoring programs.3. The ngVLA exhibits exceptional performance in the subcentimeter bands (Band 4 and Band 5), capable of producing high-fidelity images of disk substructures due to gap-opening planets within feasible integration times (less than 1 hr on-source time).It holds promise for future high-resolution surveys.4. While SKA2 shares a similar observational band with ngVLA Band 5, its ability to detect planet signatures appears somewhat less strong in comparison.However, it will serve as a valuable complement to SKA1-Mid.
In summary, both SKA and ngVLA will serve as significant complements to the existing (sub)millimeter facility, e.g., ALMA.Our research indicates that SKA and ngVLA will play a crucial role in characterizing the distribution and properties of pebbles within protoplanetary disks.In subsequent studies, we plan to further explore how to employ SKA and ngVLA in conjunction with ALMA for multiwavelength collaborative observations using spectral energy distribution methods, aiming to better constrain the distribution of dust sizes.

Figure 1 .
Figure 1.The total (absorption + scattering) optical depth (κ abs + κ sca ) vs. the maximum dust are used in our radiation transfer calculations.Different colors indicate different observation wavelengths.

Figure 2 .
Figure 2. The sky models in units of μJy beam −1 that we used for this study.In each panel, the value of α viscosity of each model is labeled in the lower right corner of each panel.The white "+" symbol indicates the position of the planet.

Figure 3 .
Figure3.Simulated SKA1-Mid (a) continuum images of the planet-disk interaction models with a Jupiter mass planet at 5 au (marked by the white "+") from the central stellar.The band wavelength is 12.5 GHz (marked at the upper left corner).The beam size (0 025 × 0 022) is marked at the lower left corner.From top to bottom, there are images with a noise level equivalent to 1 hr onsource time (marked at the upper right corner), 10 hr on-source time, and a without noise image.From left to right, the α viscosity (marked at the lower right corner) in each model is 10 −3 and 10 −5 , respectively.

Table 1
The Label, Array Names, Frequencies, Wavelength, Angular Resolutions, and Sensitivities in the Computation of the Synthetic Observations