© The Authors 2024

1 Introduction

Recent studies of young stellar objects (YSOs) have concluded that dust growth from small particles to planetesimals may occur very early in the lifetime of protoplanetary disks (Tychoniec et al. 2020; Drążkowska et al. 2023). To stop the radial drift of dust particles and allow growth to happen, dust evolution models require dust particles to be trapped within disk substructures. This requirement implies that the formation of these substructures should be well underway during the Class 0-I phase of disk evolution. While such substructures have been detected in the disks of some Class I objects (Sheehan 2020; Segura-Cox et al. 2020), these are still very limited in number. Moreover, they have not yet been observed in Class 0 protoplanetary disks. High-resolution studies of more evolved (Class II) protoplanetary disks, however, have revealed that substructures are common (ALMA Partnership 2015; Andrews et al. 2018; Isella et al. 2016; Long et al. 2018). These substructures, such as gaps, rings, arcs in cavities, and spiral arms (e.g., Casassus et al. 2013; van der Marel et al. 2013, Pérez et al. 2016; Huang et al. 2018), have proven crucial in reconciling the timescales for dust drift and planet formation, allowing larger dust particles to decelerate and grow further (Pinilla et al. 2012).

Many authors have suggested that disk substructures form due to interactions between the disk and forming planets (Dong et al. 2015; Zhang et al. 2018). Other processes, however, may also contribute to their presence (Flock et al. 2015; Zhang et al. 2015; Okuzumi et al. 2016; Takahashi & Muto 2018). The presence of substructures in evolved disks raises questions about when these substructures form and how are they linked to planet formation. The origin of these substructures, however, remains a subject of debate.

In addition to radial drift, vertical settling is another crucial process that significantly influences the evolution of dust particles in protoplanetary disks allowing them to grow and move into the mid-plane. This settling refers to the vertical motion of particles within the protoplanetary disk, driven by the balance between the gravitational force from the central star and the gas drag experienced by the particles. This settling process is influenced by several factors, including the turbulence of the disk and the sizes of the dust grains (Dullemond & Dominik 2004). Larger dust grains are generally more efficient in settling due to their greater inertia decoupling from the gas and moving into the mid-plane where they can grow further while smaller grains experience stronger gas drag and tend to stay more mixed with the gas in higher layers of the disk (Barrière-Fouchet et al. 2005). Vertical settling is expected to occur faster than radial drift and can be particularly pronounced in the inner regions of disks (Laibe et al. 2014). In evolved Class II disks, some studies have observed that the larger dust particles are already primarily located at the mid-plane, indicating that settling has already occurred (Pinte et al. 2016). On the other hand, the larger dust particles in younger disks may not have had enough time to settle completely, and the settling process may be ongoing in Class I disks with larger vertical extensions (Villenave et al. 2020, 2023). Furthermore, previous studies have provided enough evidence supporting the occurrence of grain growth in Class 0 young sources, as indicated by their low millimeter spectral indices. Notably, these studies have shown that the spectral index values are larger (ranging from 3.5 to 5) within the envelope at scales extending beyond 2000 au compared to the values (<3.5) at smaller scales (<200 au; Kwon et al. 2009; Jørgensen et al. 2007).

More recent findings suggest that young objects, particularly Class 0 YSOs, exhibit significant degrees of flaring and have considerable scale heights (Sheehan et al. 2022; Michel et al. 2022). This flaring and large scale height, irrespective of resolution or optical depth considerations, may conceal substructures in these systems or even prevent their formation. As a result, the settling of large particles is still ongoing, and much material remains in higher layers of the disk. More recently, a large study was performed in the eDisk survey (Ohashi et al. 2023) of Class 0-I objects with several new findings about their young disks. Understanding the dust properties, vertical structures, and evolution of substructures from the early disk stages is crucial for comprehending the onset and progression of planet formation, as dust evolution and grain growth play vital roles in that process.

In this study, we investigated the protobinary system NGC 1333-IRAS4A (IRAS4A), which is situated in the Perseus molecular cloud at a distance of 293 parsecs (pc; Zucker et al. 2018). The system consists of two Class 0 protostars, namely IRAS4A1 and IRAS4A2, which are separated by an angular distance of 1.8″ (Tobin et al. 2018). Both IRAS4A1 and IRAS4A2 are surrounded by an envelope with a total mass of approximately 8 M_⊙ and a total luminosity of around 5 L_⊙ (Maury et al. 2019). Both objects have very well-distinguished outflows (Santangelo et al. 2015). We present here the continuum emission of both A1 and A2. We aim to study the structure of IRAS4A1 and investigate the absence of substructures in this component using radiative models. Additionally, in a separate paper, we will discuss the line emission and the presence of complex organic molecules in the IRAS4A system, as well as the continuum analysis of IRAS4A2.

2 Observations

The observations used in this paper were obtained using the Atacama Large Millimeter-submillimeter Array (ALMA). Band 4 (1.2 mm) and Band 6 (2.1 mm) data were taken as part of the project code 2018.1.00510.S (PI: James Di Francesco). The calibration of the data was performed by the ALMA staff and was restored by the allegro team at Leiden University. For Band 4, the observations were carried out in five execution blocks spanning from October 16, 2018, to September 12, 2021. Band 6 data were acquired in four execution blocks from November 19, 2018, to September 30, 2019. The total observing time on source for Band 4 was 3.34 h, while Band 6 had a total observing time of 2.88 h. Table 1 provides additional information regarding the characteristics of the data utilized in this study.

The data reduction process was carried out using the Common Astronomy Software Applications (CASA; McMullin et al. 2007) version 5.7.0. The continuum spectral windows were separated from the line spectral windows and then averaged into eight channels for both data sets. Band 6 has 12 spectral windows centered at 264 GHz, 252 GHz, and 250 GHz with a total bandwidth of 2 GHz each. Band 4 has 15 spectral windows centered at 138 GHz, 150 GHz, and 152 GHz with a total bandwidth of 2 GHz each. The flux calibration errors are set to the nominal values of 5% at Bands 4 and 6.

Self-calibration techniques were employed for each spectral window individually using solution intervals of inf, 60 s, and 30 s. Initially, we performed phase-only self-calibration to the short baseline data, which resulted in significant improvements for Band 4 data (signal-to-noise ratio, from 88 to 380). For Band 6 data, however, only an increase in the signal-to-noise ratio (S/N) from 61 to 69 was achieved. In any case, sufficient self-calibration solutions were found, leading to enhanced data quality. Amplitude self-calibration was also performed, but we stopped after a single step for most of the spectral windows as it did not yield substantial improvements in the signal-to-noise ratio. For the long baseline data, we also performed phase self-calibration, although the improvement observed was comparatively less significant than in the short baseline data (Band 4 S/N from 42 to 86 and Band 6 S/N from 16 to 19). The reason for the lesser improvement in long baselines compared to short baselines could be attributed to a higher frequency of returns to the phase calibrator source during the long-baseline observations. Additionally, since self-calibration was exclusively applied in the same configuration, the cleaning process for long-baseline data often struggles to model the largest angular scales, even though they are present. This limitation affects the visibility data, especially considering the substantial amount of large-scale emission present in these data. Amplitude self-calibration was only applied to a few specific spectral windows due to minimal enhancements in the S/N. The final data sets were obtained after concatenating all the spectral windows in which no alignment was needed for any of them. Moreover, Very Large Array (VLA) data of IRAS4A were obtained from the VLA Nascent Disk and Multiplicity (VANDAM) survey (Tobin et al. 2016) conducted in the Perseus molecular cloud. The observations took place on October 21, 2013, employing the B-array configuration. For the correlator setup, two basebands with a total bandwidth of 4 GHz were utilized. These basebands were centered at frequencies of 36.9 GHz and 28.5 GHz, respectively. The setup was then further divided into 32 spectral windows each having a bandwidth of 128 MHz. The VLA Ka-band data in B-configuration has a shortest baseline length of 210 m and an estimated amplitude calibration uncertainty of ~10%.

The final continuum images were created using task tclean in CASA. In addition, we used the MTMFS deconvolver (Rau & Cornwell 2011) with nterms=2, together with scales of 0, 10, 30, and 50 times the pixel size (0.003″ and 0.01″ for ALMA and VLA images, respectively). Briggs weighting was found to be optimal for the purpose of this project as it provided the best compromise between sensitivity and resolution, and several Robust values were explored when making the final images. Furthermore, for the Band 4 data the uv range was modified to decrease the resolution. A smooth tapering was applied by setting uvtaper to 0.058″. Both the Band 4 and Band 6 images were convolved to have the same 78 (milliarcsecond) beam. This common resolution allowed for a consistent analysis alongside the VLA data at 9.1 mm. Tables 2 and 3 provide an overview of the characteristics of the images for Band 4 and Band 6, along with the VLA image obtained from the VANDAM survey.

We acknowledge that the ALMA data for IRAS4A at Band 4, with its high resolution, time on source, and rms, can be favorably compared to the data obtained in the ALMA 2014 Long Baseline Campaign (LBC) Science Verification (SV) data of HL Tau at Band 6 (ALMA Partnership 2015). The data for HL Tau was specifically designed to search for substructures, a goal in common with the observation of IRAS4A1. The IRAS4A data set has a resolution of 47 mas, a time on source of 3.34 h, and an rms of 13 µJy, while the HL Tau data had a resolution of 35 mas, a time on source of 4.5 h, and an rms of 11 µJy. Given the numerous substructures identified in the HL Tau disk and the comparable nature of the data, one would expect these observations to be sufficient for detecting substructures in the IRAS4A1 disk.

3 Results

Figure 1 displays the continuum images obtained from the observations. IRASA1 and IRAS2A2 are well resolved in both the Band 4 and Band 6 images. The majority of the submillimeter emission originates from within each of these two sources. There is, however, additional faint emission observed between and surrounding A1 and A2, indicating some form of structure between the two sources. This structure is particularly evident at 2.1 mm and 1.2 mm wavelengths, but not at 9.1 mm, which might be related to the lower sensitivity to thermal dust emission at 9.1 mm (see Fig. A.1). The origin of this extended emission remains unknown, but the material could potentially be associated with the surrounding molecular cloud or with some diffuse envelope-core material at these scales. In contrast, the emission surrounding A1 or A2 is likely originating from the inner envelope or a very optically thick flared disk. Moreover, the brightness peak emission from A1 is lower than that from A2 at 1.2 mm, contrary to what is observed at 2.1 mm and 9.1 mm (see Sect. 3 for the flux values). One possible explanation for this discrepancy could be that both sources have different scale heights and different optical depths. Despite the objects’ similar age, Band 6 may be tracing different layers in A1 and A2, possibly not corresponding to the mid-plane.

Furthermore, our ALMA images have been thoroughly examined, and no additional compact objects, such as low-mass companions or distant galaxies, have been detected within the field of view (>3σ). Due to the high sensitivity and resolution of our data, it is highly improbable that any such objects have been missed. This suggests that A1 and A2 are unlikely to be part of a binary with a separation greater than 20 au, which is the long axis of the beam. However, is important to note that our data are not sensitive enough to detect a star lacking a circumstellar disk.

Lastly, a distinct asymmetry is observed for A1 in the 1.2 mm image, which is not apparent in the 2.1 mm and 9.1 mm images. The cause of this asymmetry warrants further investigation as it may provide valuable insights into the vertical structure of its Class 0 protoplanetary disk.

The radial profiles from these images were obtained by averaging the emission in elliptical rings for both sources, and the central position of the radial profile of A2 was determined based on the peak emission in the respective images. Since there is an asymmetry in the A1 source, the central position of the radial profiles was determined by a Gaussian fit using imfit. Although a slight bias might remain in the Gaussian fit, it was considerably less pronounced than that seen when using the peak emission center. Consequently, the center from the Gaussian fit is likely much closer to the actual center of the source. For A1, the inclination and position angle values were set to 20 degrees and 96 degrees (NE direction, from the north axis moving towards the east), respectively, as reported for the outflow in Ching et al. (2016). On the other hand, for A2 we took the inclination to be 14 degrees (Ching et al. 2016), while the position angle was taken from measurements on the inner outflow (Chuang et al. 2021) (122 degrees NE). Finally, the brightness temperature values were calculated by applying the full Planck equation to the radial intensity profiles as indicated by Rybicki & Lightman (1979); the concept refers to the temperature of a blackbody with the same brightness at that specific frequency. Moreover, two additional radial profiles were generated, representing the spectral indices between 9.1 mm and 2.1 mm, as well as between 2.1 mm and 1.2 mm. Figure 2 displays both the brightness temperature profiles and the spectral indices for A1 and A2.

By examining Fig. 2, we can observe the behavior of the spectral indices for A1 and A2. For A1, the spectral indices are very low at the center of the source. As the radius increases, however, these indices gradually become larger. This continues until the noise of the 9.1 mm image starts to dominate the emission. Comparing the spectral index between 2.1 mm and 1.2 mm for A1 with that of A2, we find that the index for A1 remains consistently below two across most radii. On the other hand, the spectral indices for A2 are consistently above two throughout the range of radii considered. This discrepancy suggests that the emission from A1 is significantly more optically thick than that from A2 at these wavelengths. The difference in spectral indices between A1 and A2 implies differences in the physical properties of the two sources. For example, A1 may have a denser and more optically thick environment, which affects the observed spectral behavior. Furthermore, dust self-scattering might be affecting the inner regions of the IRAS4A1 disk. Additionally, some free-free emission might be increasing the brightness temperature in the inner regions of the A1 VLA image, affecting the spectral index from 9.1 mm to 2.1 mm.

In the recent study conducted by Galametz et al. (2019) using independent measurements from the CALYPSO sample at 1.3 and 3.2 mm, the authors reported discovering remarkably low values of spectral indices (<2.0) within the inner regions of the IRAS4A1 envelope, specifically at distances of less than 200 au. This is in agreement with our spectral index values of the IRAS4A1 inner regions. Additionally, Galametz et al. (2019) also observed higher spectral index values extending up to 2000 au, which was attributed to grain growth processes occurring within the envelope. It is crucial to acknowledge that our high-resolution image might be causing the extended component of the envelope to be resolved out, thereby making it difficult to measure the spectral index of this particular component. Jørgensen et al. (2007) previously pointed out that when extracting emission from the envelope, the spectral index of compact components would be flattened. Then, spectral index values below 3.5 at smaller radii could be indicative of the presence of another component, most likely a disk.

Fig. 1 Images of IRAS4A at 1.2 mm, 2.1 mm, and 9.1 mm imaged at 78 mas resolution. The central RA and Dec positions for Band 4 and Band 6 are 03:29:10.510 and +31.13.31.010, respectively. For the VLA image, the central Dec position is 03:29:10.502 and the central RA position remains the same, aligning with +31.13.31.010. Both sources are visible at a separation of 1.8″, with some surrounded leftover emission seen between them at 1.2 and 2.1 mm but not at 9.1 mm (see Fig. A.1, where the image color scale was changed to show the extended emission). Moreover, the peak emission of IRAS4A2 is larger than that of IRAS4A1 at 1.2 mm. The emission of both sources in the VLA image is more radially compact, though it is very faint for IRAS4A2.

Fig. 2 Brightness temperatures and spectral indices across the entire radius of A1 and A2, illustrating their emission patterns. Upper right and left panels: brightness temperatures at 1.2 mm, 2.1 mm, and 9.1 mm of A1 and A2, respectively. Lower right and left panels: spectral indices between 9.1–2.1 mm and 2.1–1.2 mm of A1 and A2, respectively. In the inner parts of the A1 source, some free-free emission might be present at 9.1 mm. Moreover, the brightness temperature in A1 at 1.2 mm is lower than that at 2.1 mm, which, as seen from the spectral index, might indicate very optically thick emission at those wavelengths and small dust particles (<1 mm). We only considered the statistical uncertainties for the brightness temperature and spectral index.

3.1 Multiwavelength analysis of a Class 0 young stellar object

In this study, we adopted the hypothesis that the emission detected in our high-resolution images originates primarily from a disk rather than the envelope. The reason behind this is that the high resolution of our imaging may result in the loss of most of the emission from the extended components (envelope) and that, as mentioned before, our findings of low values of the spectral indices in Fig. 2 further support the notion of a disk scenario. Of course, it needs to be noted that some emission coming from the inner envelope might still be contributing to the total emission. While there might still be some confusion within the envelope, we unfortunately did not account for a dynamical distinction. Separating the continuum from the envelope is challenging, and due to the optical depth, analyzing the lines becomes quite limited. Additionally, to be aligned with previous research, we also consider the disk to be flared, similarly to observations and results found in Class 0 young stellar objects (YSOs) using ALMA data (Sheehan et al. 2022; Michel et al. 2022) and as suggested by edge-on observations (Villenave et al. 2020, 2023) and, recently, the eDisk survey (Ohashi et al. 2023).

Protoplanetary disks are commonly expected to have millimeter- or even centimeter-sized dust particles. Because of such grain growth, the albedo of the dust can be high at millimeter wavelengths, indicating that scattering plays a significant role in the opacity of the dust emission. When scattering is a dominant factor, the spectral index of the dust emission can no longer be directly associated with a spectral index of the dust opacity (i.e., β; e.g., Sierra & Lizano 2020; Zhu et al. 2019). To analyze the spectral energy distribution (SED) of protoplanetary disks properly, it is crucial to consider both absorption and scattering effects in the dust opacity. To include the scattering effect, we can write the source function in the radiative transfer equation as $$S_v(T)={\omega }_v{J}_v+\left(1-{\omega }_v\right)B_v(T),$$(1)

where J_ν is the local mean intensity and ω_ν is the albedo, defined by the scattering coefficient and the absorption coefficient as ${\omega }_v=\frac{{\sigma }_v}{{\kappa }_v+{\sigma }_v}$. We can approximate this to the analytical solution found in Miyake & Nakagawa (1993) assuming a disk as a vertically isothermal slab and with isotropic scattering: $$J_v=B_v(T)\left[1+f\left(t,{\tau }_v,{\omega }_v\right)\right],$$(2)

where $$f\left(t,{\tau }_v,{\omega }_v\right)=\frac{\exp \left(-\sqrt{3}{ϵ}_{v}t\right)+\exp \left(\sqrt{3}{ϵ}_v\left(t-{\tau }_v\right)\right)}{\exp \left(-\sqrt{3}{ϵ}_v{\tau }_v\right)\left({ϵ}_v-1\right)-\left({ϵ}_v+1\right)};$$(3)

t is the optical depth variable and τ_ν = Σ_dustχν; both are measured perpendicular to the disk mid-plane. Also, ${ϵ}_v=\sqrt{1-{\omega }_v}$. Considering inclination effects by correcting the optical depth by the inclination angle (i) of the disk, we reach the emergent specific intensity obtained by Sierra et al. 2019: $$I_v=B_v(T)\left[\left(1-\exp \left({\tau }_v/\mu \right)\right)+{\omega }_{v}F\left({\tau }_v,{\omega }_v\right)\right],$$(4)

where $$\begin{array}{c}F\left({\tau }_v,{\omega }_v\right)=\frac{1}{\exp \left(-\sqrt{3}{ϵ}_v{\tau }_v\right)\left({ϵ}_v-1\right)-\left({ϵ}_v+1\right)}\\ \text{X}\left[\frac{1-\exp \left(-\left(\sqrt{3}{ϵ}_v+1/\mu \right){\tau }_v\right)}{\sqrt{3}{ϵ}_v\mu +1}+\frac{\exp \left(-{\tau }_v/\mu \right)-\exp \left(-\sqrt{3}{ϵ}_v{\tau }_v\right)}{\sqrt{3}{ϵ}_v\mu -1}\right].\end{array}$$(5)

It is important to mention that isotropic scattering is assumed for these equations, which may be an incorrect approximation for 2πa ⩾ λ. To reduce the effect of the approximation, we replaced the scattering coefficient in all equations with an effective scattering coefficient in the following form (Ishimaru 1978; Birnstiel et al. 2018): $${\sigma }_v^{eff}=\left(1-g_v\right){\sigma }_v,$$(6)

where g_ν is the asymmetry parameter, that is, the expectation value of cos θ, where θ is the scattering angle (e.g., Ishimaru 1978; Birnstiel et al. 2018). The values of g_ν depend on the dielectric properties of the dust particles. For our calculations, the values obtained by Birnstiel et al. (2018) for ${\sigma }_v^{eff}$ were used.

In our analysis, the particle size distribution is assumed to follow a power law with a slope (n(a) ∝ a^−p), where p is commonly assumed to be 3.5 according to measurements of the ISM (Mathis et al. 1977). Also, the DSHARP opacity data (Birnstiel et al. 2018) were employed, which considers particles without porosity and a composition of 20% water fraction by mass, 32.91% astronomical silicates, 7.43% troilite, and 39.66% refractory organics.

Equation (5) then ultimately depends on only three free parameters: dust temperature (T_dust), surface density (I_dust), and particle size (a_max). With three or more observed wavelengths, it becomes possible to solve the equation and obtain estimates for the three free parameters (T_dust, Σ_dust, a_max).

It is important to note that this model assumes a single temperature at each radius within the disk. This assumption generally holds when most of the dust is settled in the disk’s mid-plane. In cases where the emission is originating from an envelope or a flared disk involving different layers, however, this assumption may not be valid. So, it is worth noting that the temperature structure within protoplanetary disks can be complex, particularly if there are significant vertical temperature gradients or if different layers of the disk are contributing to the observed emission at different wavelengths. In these situations, a more sophisticated modeling approach that considers the vertical structure and temperature gradients within the disk would be necessary to interpret the observed SED.

A multiwavelength analysis similar to ours here was previously performed on HL Tau using four images between 8 mm and 0.9 mm (Carrasco-González et al. 2019) by simplifying the spectral behavior of the extinction coefficient using a power law. After that, it has been used in several other works (e.g., Macías et al. 2021; Sierra et al. 2021; and Guidi et al. 2022), including the work presented in this paper, using the exact values of the dust opacity at each wavelength. This model is a first approach in determining the dust properties around a Class 0 YSO such as A1.

A Bayesian approach was employed to obtain the posterior probability distributions of the model parameters (a_max, T_dust, Σ_dust) at each radius. To achieve this, a standard log-normal likelihood function was used, which is defined as follows: $$\ln p(\overline{I}(r)\mid \text{Θ})=-0.5\underset{i}{{{\displaystyle \sum }}^{\text{}}}\left({\left(\frac{{\overline{I}}_i-I_{m,i}}{{\widehat{\sigma }}_{\overline{I},i}}\right)}^2+\ln \left(2\pi {\widehat{\sigma }}_{\overline{I},i}^2\right)\right),$$(7)

where Ī is the azimuthally averaged intensity at radius r and at frequency ν_i, I_m,i is the model intensity from different combinations of the three free parameters at a radius r, and Θ is the vector of the three free parameters. In addition, we assumed that the uncertainty ${\widehat{\sigma }}_{I,i}$ at radius r is $${\widehat{\sigma }}_{\text{¯},i}=\sqrt{{\sigma }_{\overline{I,i}}^2+{\left({\delta }_i{I}_i\right)}^2}\text{, }$$(8)

where ${\sigma }_{\overline{I},i}^2$ is the error of the mean, obtained from the azimuthally averaged intensity profiles (see Sect. 3), and δ_i is the flux calibration error at each frequency.

Figure 3 shows the analysis we performed; a model grid of intensities was created using various dust parameters. To infer the physical parameters of the dust particles, we compared the observed intensity at each radius with the expected spectral energy distribution (SED) derived from different combinations of the three free parameters in Eqs. (5) and (6) (a_max from 0.001 to 10 cm, T_dust from 0.1 to 250 K, Σ_dust from 0.1 to 1000 g cm⁻²). In order to better match the observational data, the probability distribution of each parameter was plotted, along with the corresponding expected value (represented by the red curve in Fig. 3). The expected value of each parameter was obtained by $$E(X)=\frac{{\displaystyle {\sum }_i{X}_{i}P\left(X_i\mid \overline{I}(r)\right)}}{{\displaystyle {\sum }_{i}P\left(X_i\mid \overline{I}(r)\right)}},$$(9)

where X_i is each value in all the parameters inside our grid, and P(X_i|Ī(r)) is the marginalized posterior probability of each parameter in every single cell of the grid.

Multiple equally likely solutions close to each other were found for the A1 source, which equally explained the observed intensity. All the possible solutions have similar Σ_dust and dust temperatures, which is explained in Zhang et al. (2023) as there are no strong Mie interference patterns when 2πa < λ. Finally, the optical depth values were derived from the analysis to provide insights into the dust properties at different locations within the disk. Figure 3 shows the dust parameters, the optical depths at each wavelength, and the intensities of the observations compared with the ones obtained from the model.

From Fig. 3, it is evident that the A1 disk exhibits high optical thickness at the inner radii of A1, which poses challenges in fitting the dust parameters accurately. This observation suggests that the disk is highly unstable and contains very small dust particles (hundreds of microns in size) relatively to dust grain sizes in protoplanetary disks. Notably, the derived temperature from the multiwavelength analysis in A1 appears to be higher compared to other Class II disks analyzed using similar methods (Macías et al. 2021; Sierra et al. 2021; Guidi 2019; and Carrasco-González et al. 2019). This discrepancy may be attributed to the young age of the source and other processes occurring within the system, such as infalling material that can contribute to the elevated temperature of the dust particles, viscous heating or even back warming by the envelope (Natta 1993).

Furthermore, A1 displays a notably high dust surface density and mass in comparison to Class II disks. This result aligns with the notion that a significant portion of the material remains distributed as submillimeter particles surrounding the star rather than having settled and grown in the disk’s mid-plane, where it cannot be detected by our observations due to high optical depths. The high mass inferred is expected for a very young source such as A1, which is likely to be very gravitationally unstable as it still has a substantial circumprotostellar mass not yet accreted by the central star. We note that the particle sizes found in other disk studies often vary significantly depending on the presence of substructures, which are not detectable in the A1 source. Moreover, the particle sizes observed in other Class II disks tend to be larger (centimeter-sized particles) compared to the 0.1 mm particles found in A1. This disparity can be attributed to the different evolutionary stages of the disks, with the dust in the other disks having evolved and settled more in the mid-plane.

Figure 3 indicates that the material flowing in to form the disk already contains large dust particles (>10 µm) compared to the average ISM dust sizes. This suggests widespread grain growth across the entire disk radius. However, as one approaches the mid-plane and the central star, particles tend to become larger. These large dust particles compared to the ISM particles imply that grain growth is not limited to the mid-plane and that it also occurs in the flared regions of the disk where infall is likely the process that triggers this growth. Additionally, the increase in error at the outer region of the A1 disk is a result of the spatial sensitivity of the VLA image.

We compared the temperature profile of A1 with other Class 0 sources and found that A1’s temperature agrees with those derived from CO and H₂CO snow lines in IRAS 04302 (Class I) and L1527 (Class I-O) by van’t Hoff et al. (2020). Comparing with models from Yang et al. (2017) on the Class 0 proto-star BHR71, however, we note that the derived densities in A1 are at least an order of magnitude lower. This difference could potentially stem from the observations utilized by Yang et al. (2017) being of shorter wavelengths (Herschel) that are more sensitive to the cloud, the surrounding envelope, and smaller dust grains.

Concerning particle sizes, our analysis indicates that the inner disk of A1 comprises particles nearly 0.3 mm in size. This suggests that the dust size distribution in the disk is primarily characterized by larger particles when compared to typical interstellar medium (ISM) dust sizes. However, in comparison to pebbles found in more evolved disks, these particles are relatively small. This suggests that although some dust growth has already occurred, the process is still ongoing. Several studies focused on Class 0 objects have measured dust sizes in the envelope using low dust emissivity indices, revealing that grain growth might already happen in this Class 0 object, maybe even up to millimeter-sized particles (Jørgensen et al. 2009; Galametz et al. 2019). These findings diverge from the multi-wavelength analysis in our work that shows smaller dust particles in IRAS4A1.

The findings presented in Fig. 3 do not provide a definitive explanation for the spectral index below two in Fig. 2. The particle sizes around 0.1 mm align with what is expected for dust self-scattering, indicating the presence of low spectral indices (Liu 2019). Nevertheless, these observed values can also be rationalized by considering a highly optically thick disk within r < 60 AU, where the inner layers are warmer than the outer layers. This scenario not only aligns with the observations but also corresponds to the outcomes illustrated in Fig. 3, showcasing the high optical depth across all radii.

Fig. 3 Multiwavelength analysis results of A1. Left panels: probability distributions of dust parameters (amax, Tdust, Σdust) of A1; the red line in each plot is the expected value of each of the dust parameters. Middle panel: optical depths at each radius of the three different wavelengths used in this work. Right panels: comparison of radial intensities from the observations and the model at each wavelength. The plots show very high temperatures and surface densities along the radius of the disk. Particle sizes have increased and are large in comparison with the ISM but still very small for the mid-plane of protoplanetary disks, in agreement with the expectation that not much settling has occurred in very young protoplanetary disks. Lastly, A1 shows optically thick emission at Bands 6 and 4 at the inner radii.

3.2 Generic gap models with large scale heights

To explain the absence of observed substructures and the observed asymmetry in the IRAS4A1 disk, we employed radiative transfer models using RADMC-3D (Dullemond et al. 2012). For these models we assumed that instead of observing a highly optically thick envelope with an embedded disk, IRAS4A1 is actually a flared disk with a significant scale height (the surrounding envelope has been resolved out in these high-res images; see Sect. 4). When considering a greater scale height and flaring in the disk, it is crucial to differentiate between flared disks, which represent an equilibrium configuration of orbiting material, and an infalling model. In this study, we modeled the flared disk solely from the perspective of the dust continuum, without incorporating a dynamic approach such as infalling or rotational motions. This flaring effect can create an asymmetry in the disk and make it challenging to detect substructures, if they exist. The combination of disk inclination, large scale heights, and optically thick emission, even at Band 6, contribute to this effect. Evidence supporting the presence of a highly flared disk instead of an envelope has been observed in the Class 0 Protostar L1527 IRS by Sheehan et al. (2022).

To test this assumption, a model of the dust continuum emission at 1.2 mm was constructed using RADMC-3D. Initially, “generic gap models” inspired by the disk of HD 163296 were made to investigate the disappearance of substructures with increasing scale height. Subsequently, we developed a specific model to the IRAS4A1 disk to reproduce the observed asymmetry at 1.2 mm in combination with the absence of substructures. These radiative transfer models allowed us to perform a detailed examination of the disk’s dust evolution and provide insights into its vertical structure.

For the generic gap models, we fixed certain parameters based on previous studies of HD163296. For the star, the parameters in Table 1 from (Andrews et al. 2018) were used: M_* = 2.04 M_⊙, L_* = 17 L_⊙, T_* = 9332 K, and a distance of 101 pc. The positions of the two most prominent gaps were taken to be 49 and 86 au with a fixed width of 10 and 8 au, respectively.

The disk model was taken to have an inclination of i = 46.7 Å, a position angle of 133.3 Å, and a dust mass of 0.039_⊙ from Dullemond et al. (2020). In addition, a size of 110 au for the disk was chosen. Inside RADMC-3D, a generic protoplanetary disk model was used, with the scale height varied in each model. To incorporate the DSHARP dust particle opacities, the optool software (Dominik et al. 2021) was utilized, allowing for their utilization within RADMC-3D. Finally for completeness, RADMC-3D calculated the dust temperature using the density distribution for the generic protoplanetary disk model as follows: $$\rho (r,z)=\frac{{{\displaystyle \sum }}^{\text{}}(r)}{H_p\sqrt{2\pi }}exp\left(-\frac{z^2}{2H_p^2}\right),$$(10)

where r is the distance to the star from the disk, Σ(r) is the dust surface density, and H_p is the scale height of the dust disk.

The scale height (H_p) in the generic protoplanetary disk model follows a power-law dependence on the radial distance as follows: $$H_{\text{p}}=H_{100}{\left(\frac{r}{100\text{AU}}\right)}^{1+\text{Ψ}},$$(11)

where Ψ is the flaring index, with a predefined value 0.14, and H₁₀₀, is the value of the scale height at a distance of 100 au from the central star.

The scale height parameter was increased in the generic gap models until the substructures disappeared due to shadowing, obscuration, and-or contrast effects. Figure 4 shows the images of these models together with a cut through their major and minor axes. In Fig. 4, it is evident that substructures present in young Class 0 disks are challenging to observe, if present, due to the large scale heights that these disks may exhibit. The cuts shown in Fig. 4 provide additional insights into the behavior of the disk at different scale heights. Along the major axis, even at a low scale height of 0.05, a strong flattening effect on the rings and gaps is observed. This effect is highly dependent on the inclination and position angle of the disk. The intensity variations along the minor axis reveal another interesting aspect: in the SW part of the disk, a lack of intensity is observed, an asymmetry caused by a large vertical structure in the disk, also seen in other sources, such as Lee et al. (2021) and Lin et al. (2023). The direction of this asymmetry is determined by the orientation of the modeled disk. Furthermore, as the scale height increases, both the depth of the gaps and the visibility of substructures begin to flatten along the minor axis too. Eventually, there is a point where substructures (≤ 10 au) can no longer be distinguished. This example demonstrates the impact of a highly flared disk on the visibility and discernibility of substructures, if any, in a young protoplanetary disk.

3.3 Large scale height and very flared disk models of IRAS4A1

To investigate the asymmetry observed in the 1.2 mm image of the IRAS4A1 disk, additional modeling was performed in RADMC-3D. The objective was to determine whether or not the observed asymmetry could be reproduced in a large-scale-height flared-disk scenario. To set up the RADMC-3D models, we fixed specific parameters. Due to the difficulty of determining the stellar properties directly from the literature for a highly embedded Class 0 object such as IRAS4A1, average values of stellar properties in a number of Class I systems were obtained from Tables 1 and 2 in Fiorellino et al. (2023). These average values include the stellar mass (1.55 M_⊙), radius (2.1 R_⊙), and effective temperature (3700 K). The inclination and position angle of the disk were fixed at 20° and 99 Å, respectively. The dust mass in the disk was taken from the multi-wavelength analysis, resulting in a value of ${0.11}_{-0.04}^{+0.08}{M}_{\odot }$. The scale height in the RADMC-3D models for the IRAS4A1 disk was initially set to H/R = 0.3, based on the appearance of asymmetry in the generic gap models. In addition to this base model, eight more models were created. three with a fixed scale height, three with a fixed high flaring index (Ψ = 1.3), and two models with reduced gap widths. This variety allowed us to explore different scale heights within the context of a consistently high flaring profile. Figure 5 shows the corresponding cuts through the major and minor axes in all eight models. The IRAS4A1 observation and the model that best reproduces its intensity along the major and minor axis are shown in Fig. 6. By examining the outcomes of these various models, we can observe the influence of a gap presence, large flaring index, and large scale heights on the observed asymmetry and young Class 0 sources like the IRAS4A1 disk.

From the radiative transfer models of IRAS4A1, it is evident that an asymmetry is formed on the north (compared to south) part of the disk at large scale heights. The inclination and position angle in the models greatly influence the resulting asymmetry, emphasizing the uncertainties in these results. Furthermore, the simplicity of the model employed in this study may limit its ability to accurately reproduce the complexities of a Class 0 young stellar object such as IRAS4A1. Nevertheless, the intensity profiles along the major and minor axes suggest the presence of “substructures” or other unknown processes occurring in the actual observations, as most models appear flat unless a gap is included. In the upcoming paragraphs, we will speculate about the substructure scenario in the IRAS4A1 disk, although it is possible that something else is shaping the intensity profiles along the major and minor axis.

The models with a very small gap exhibit intensity profiles that more closely resemble the observed profile at 1.2 mm in both the major and minor axes. This difference may indicate that the gaps at these early stages are still forming and that we will need still higher resolution to see them. Regardless of whether IRAS4A1 is indeed a flared disk, a lower limit on the scale height for generating an asymmetry can be established (H/R > 0.3). We note that the intensity profiles in the figure are normalized, as the primary goal of this study is not to replicate the flux of the IRAS4A1 source precisely, but rather to provide insights into the earliest stages of disk and planet formation. Nevertheless, if our observations capture emission from higher layers in the disk, and if the emission remains highly optically thick, it may be challenging to detect substructures with ALMA at the available resolution.

Fig. 4 Generic gaps models using RADMC-3D, with cross-sections made along both the major and minor axes. Top panels: generic gap models in RADMC-3D increasing the scale height up to H/R 0.3. Bottom panels: normalized intensities of major axis and minor axis cuts of the generic gap models. An asymmetry in the minor axis of the disk becomes more prominent as the scale height increases. In both the minor and major axes, substructures start to flatten with increasing scale height, and after 0.3 they become barely visible.

Fig. 5 Normalized intensities of major axis and minor axis cuts of the IRAS4A1 RADMC-3D models. Top panels: scale height of each model was fixed to H/R = 0.3. Middle panels: Ψ flaring index values were fixed to 1.3 in the RADMC-3D models. Bottom panels: scale height was fixed to 0.3 and flaring index to 1.3, but the gap width was reduced in each model. The “shoulder” shape and the asymmetry of the disk are better reproduced using a large scale height, a large flaring index, and a small gap width. Unfortunately, an excess of emission is found in the models and the small size of the gap needed cannot be observed at the current resolution.

Fig. 6 Comparison between the best-fit RADMC-3D model image and the ALMA observation at Band 6 of IRAS4A1. Left: RADMC-3D model 8 convolved with a 0.078″ × 0.031″ beam. Considering both the major and minor axis cuts, the substructures disappear in the major axis at a scale height of 0.3. The gap radial location is set at 40 au with a small width of 4 au. Right: ALMA continuum unconvolved emission image at 1.2 mm with a zoomed-in view of the A1 source. The respective intensities of each images is shown in the upper color bars.

4 Discussion

The inferred large scale height (H/R > 0.3) in IRAS4A1 has significant implications for planet formation. Despite the fast settling expected during the disk’s lifetime, the optical thickness and asymmetry observed at Band 6 indicate the presence of material with varying grain sizes in higher layers of the disk. This result implies that settling is still ongoing for millimeter-sized particles. Indeed, this state is expected considering the settling timescales (<1 Myr; Dullemond & Dominik 2004) and the estimated dynamical age of the outflows in IRAS4A (around 0.01 Myr; Taquet et al. 2020). Furthermore, settling and radial drift are likely acting together during these early stages of dust evolution and growth in the disk.

The large scale height of the disk may also obscure young substructures, as suggested by models, particularly when combined with very narrow substructures measuring less than 4 au in size. While we are unable to directly resolve substructures (i.e., gaps) in the disk, our models suggest that some must be present to explain the observed bumps in the radial profile of IRAS4A1. These small-scale features are challenging to observe directly with current resolution capabilities (no substructures observed in IRAS4A1), but their presence at these early stages could indicate two possibilities. Firstly, if these substructures are caused by planet-disk interactions, it suggests that planets formed nearly instantaneously after the collapse of the molecular cloud. Furthermore, given sufficient time, substructures within protoplanetary disks are expected to widen, resembling those observed in other systems. This widening occurs as planets grow in size by accreting and carving out material from their surroundings.

We note that substructures can potentially arise from mechanisms other than planet-disk interaction. If this is the case, it introduces an intriguing possibility. In this scenario, the substructures initially form early on and are narrow, as indicated by the narrow gap width in the models. Given settling-induced growth of dust particles and other processes occurring within these narrow substructures, substructure formation through alternative mechanisms may itself trigger planet formation within such gaps. Hence, planet formation could take place very early in the disk’s evolution, following an evolutionary path similar to the first scenario.

An additional crucial aspect to consider is the large flaring index observed in IRAS4A1. If the “shoulder” position around the 20–40 AU range – which served as the basis for defining the gap position for the models – is pointing to a substructure, this substructure could be relatively close to the disk’s mid-plane.

For example, smaller scale heights may occur in the center of the disk, meaning that the closer to the center a substructure forms, the easier it would be to detect it in a very flared disk. Consequently, any mechanism responsible for carving out these substructures probably starts in the mid-plane and is unable to reach large scale heights, as seen in other protoplanetary disks when scattered light observations and submillimeter observations with ALMA are compared. It is important to note, however, that the combination of large scale heights and a large flaring index could still hide further substructures in the outer radii of the IRAS4A1 disk. On the other hand, if there are no substructures in the flared disk of IRAS4A1, planet formation may then occur only at a later stage when larger particles have already settled in the disk’s mid-plane, taking into account the timescales required for settling.

Recently, similar results including the “shoulders”, asymmetries, and large scale heights were found in YSOs in studies by the eDisk survey team (Ohashi et al. 2023). Regardless of the specific dynamics within the IRAS4A1 disk, it is becoming evident that Class 0 YSOs exhibit flared disks with significant scale heights, providing valuable insights into the planet formation process.

5 Summary and conclusions

We have shown high-resolution ALMA images (78 mas) of the IRAS4A binary system in Bands 4 and 6. The key findings of this paper can be outlined as follows:

• No substructures were detected in either A1 or A2 at the current resolution;

• Analysis of spectral indexes and brightness temperatures indicated that A1 is significantly more optically thick than A2;

• A multiwavelength image analysis was carried out showing the dust parameters in A1. The expected values of the dust parameters inferred high temperatures (>50 K), high surface densities (>10 g cm⁻²), and large dust size particles (>30 µm) at all radii (< 60 au) in the IRAS4A1 disk. In addition, the analysis showed high optical depth in the inner disk in Band 6 and Band 4;

• Radiative transfer models using RADMC-3D have shown that a minimum scale height of H/R > 0.3 is adequate to render the substructures invisible and produce an asymmetry in the disks. Moreover, the models that incorporated a narrow gap around 34–50 au and increased flaring index, provided better matches to the observed intensity profiles, suggesting the presence of potential hidden substructures within a very flared disk in the IRAS4A1 system, even in these early stages of disk formation.

Observations with high resolution and sensitivity at centimeter wavelengths with the ngVLA can help unveil any substructure that might exist in IRAS4A1 and give us a better understanding of substructure formation in young protoplanetary disks.

Acknowledgements

We thank the referee for the very constructive comments. We also thank Dominique M. Segura-Cox for the useful discussion. We acknowledge assistance from Allegro, the European ALMA Regional Centre node in The Netherlands. This paper makes use of the following ALMA data: ADS/JAO.ALMA#2018.1.00510.S. ALMA is a partnership of ESO (representing its member states), NSF (USA) and NINS (Japan), together with NRC (Canada), MOST and ASIAA (Taiwan), and KASI (Republic of Korea), in cooperation with the Republic of Chile. The Joint ALMA Observatory is operated by ESO, AUI/NRAO and NAOJ. E.G.C. acknowledges support from the National Science Foundation through the NSF MPS-Ascend Fellowship Grant number 2213275. L.W.L. acknowledges support from NSF AST-1910364 and NSF AST-2108794.

Appendix A Faint emission between A1 and A2

Fig. A.1 Images of IRAS4A at 1.2 mm, 2.1 mm, and 9.1 mm imaged at 78 mas resolution. The color scale was changed to better observe the faint emission between the two targets at B6 and B4.

References

All Tables

All Figures

Fig. 1 Images of IRAS4A at 1.2 mm, 2.1 mm, and 9.1 mm imaged at 78 mas resolution. The central RA and Dec positions for Band 4 and Band 6 are 03:29:10.510 and +31.13.31.010, respectively. For the VLA image, the central Dec position is 03:29:10.502 and the central RA position remains the same, aligning with +31.13.31.010. Both sources are visible at a separation of 1.8″, with some surrounded leftover emission seen between them at 1.2 and 2.1 mm but not at 9.1 mm (see Fig. A.1, where the image color scale was changed to show the extended emission). Moreover, the peak emission of IRAS4A2 is larger than that of IRAS4A1 at 1.2 mm. The emission of both sources in the VLA image is more radially compact, though it is very faint for IRAS4A2.

In the text

Fig. 2 Brightness temperatures and spectral indices across the entire radius of A1 and A2, illustrating their emission patterns. Upper right and left panels: brightness temperatures at 1.2 mm, 2.1 mm, and 9.1 mm of A1 and A2, respectively. Lower right and left panels: spectral indices between 9.1–2.1 mm and 2.1–1.2 mm of A1 and A2, respectively. In the inner parts of the A1 source, some free-free emission might be present at 9.1 mm. Moreover, the brightness temperature in A1 at 1.2 mm is lower than that at 2.1 mm, which, as seen from the spectral index, might indicate very optically thick emission at those wavelengths and small dust particles (<1 mm). We only considered the statistical uncertainties for the brightness temperature and spectral index.

In the text

Fig. 3 Multiwavelength analysis results of A1. Left panels: probability distributions of dust parameters (amax, Tdust, Σdust) of A1; the red line in each plot is the expected value of each of the dust parameters. Middle panel: optical depths at each radius of the three different wavelengths used in this work. Right panels: comparison of radial intensities from the observations and the model at each wavelength. The plots show very high temperatures and surface densities along the radius of the disk. Particle sizes have increased and are large in comparison with the ISM but still very small for the mid-plane of protoplanetary disks, in agreement with the expectation that not much settling has occurred in very young protoplanetary disks. Lastly, A1 shows optically thick emission at Bands 6 and 4 at the inner radii.

In the text

Fig. 4 Generic gaps models using RADMC-3D, with cross-sections made along both the major and minor axes. Top panels: generic gap models in RADMC-3D increasing the scale height up to H/R 0.3. Bottom panels: normalized intensities of major axis and minor axis cuts of the generic gap models. An asymmetry in the minor axis of the disk becomes more prominent as the scale height increases. In both the minor and major axes, substructures start to flatten with increasing scale height, and after 0.3 they become barely visible.

In the text

Fig. 5 Normalized intensities of major axis and minor axis cuts of the IRAS4A1 RADMC-3D models. Top panels: scale height of each model was fixed to H/R = 0.3. Middle panels: Ψ flaring index values were fixed to 1.3 in the RADMC-3D models. Bottom panels: scale height was fixed to 0.3 and flaring index to 1.3, but the gap width was reduced in each model. The “shoulder” shape and the asymmetry of the disk are better reproduced using a large scale height, a large flaring index, and a small gap width. Unfortunately, an excess of emission is found in the models and the small size of the gap needed cannot be observed at the current resolution.

In the text

Fig. 6 Comparison between the best-fit RADMC-3D model image and the ALMA observation at Band 6 of IRAS4A1. Left: RADMC-3D model 8 convolved with a 0.078″ × 0.031″ beam. Considering both the major and minor axis cuts, the substructures disappear in the major axis at a scale height of 0.3. The gap radial location is set at 40 au with a small width of 4 au. Right: ALMA continuum unconvolved emission image at 1.2 mm with a zoomed-in view of the A1 source. The respective intensities of each images is shown in the upper color bars.

In the text

	Fig. A.1 Images of IRAS4A at 1.2 mm, 2.1 mm, and 9.1 mm imaged at 78 mas resolution. The color scale was changed to better observe the faint emission between the two targets at B6 and B4.
In the text