© The Authors 2024

1 Introduction

The formation of exoplanets heavily depends on the amount of material present in a planet-forming disk. Hence, efforts have been made to determine the total disk mass and link this to planet formation occurring inside them. The most common disk mass tracer used is the millimeter dust continuum of the disk. Assuming the canonical gas-to-dust mass ratio of 100 (Bohlin et al. 1978), the total disk mass can be determined after making a few assumptions on the grain properties, most notably the grain emissivity (Beckwith et al. 1990). Many population studies on dust masses have been done with the Atacama Large Millimeter/submillimeter Array (ALMA; see, e.g., Ansdell et al. 2016; Barenfeld et al. 2016; Pascucci et al. 2016; Eisner et al. 2018; Cazzoletti et al. 2019; Anderson et al. 2022; van Terwisga et al. 2022). These studies have shown clear trends regarding the disk mass with stellar mass (Ansdell et al. 2017; Manara et al. 2023), disk radius (Hendler et al. 2020), and stellar accretion rate (Testi et al. 2022). However, it is still unknown if the dust continuum does indeed trace the total disk mass, that is, gas and dust, directly. Recent works have shown that the dust can be optically thick, underestimating the total disk mass by an order of magnitude for the most massive disks (e.g., Liu et al. 2022; Kaeufer et al. 2023), and may even be optically thick at 3 mm (Xin et al. 2023).

An understanding of the total (gas) mass of the disk is therefore much needed (see Miotello et al. 2023 for a recent overview). The most abundant molecule in a planet-forming disk is H₂. However, its faint emission at the typically low temperatures (~20 K) present in a planet-forming disk do not make it a viable disk mass tracer. Hence, other more indirect tracers are necessary. One promising tracer is the H₂ isotopologue hydrogen deuteride (HD). Using thermo-chemical models, the HD emission can be used to determine the bulk mass of a disk (see e.g., Bergin et al. 2013; Schwarz et al. 2016; Trapman et al. 2017; Sturm et al. 2023). HD J = 1−0 observed with the Herschel Space Observatory has been detected in three disks resulting in gas mass measurements (Bergin et al. 2013; McClure et al. 2016), and a dozen non-detections in Herbig disks resulting in gas mass upper limits (Kama et al. 2016, 2020). However, after the end of the Herschel mission, there are no current facilities to observe the HD J = 1−0 line, so a different tracer is needed.

Carbon monoxide (CO) has been used extensively as a tracer of both the total gas mass and size of planet-forming disks. Especially its less common isotopologues ¹³CO and C¹⁸O (e.g., Miotello et al. 2017; Long et al. 2017; Loomis et al. 2018), or even the rarer C¹⁷O, ¹³C¹⁷O and ¹³C¹⁸O isotopologues (Booth et al. 2019; Booth & Ilee 2020; Loomis et al. 2020; Zhang et al. 2020b, 2021; Temmink et al. 2023), have been used as tracers of the bulk mass. Generally, two main methods are used to determine the gas mass from CO: simple scaling relations assuming an excitation temperature and optically thin emission (e.g., Loomis et al. 2018), and using thermo-chemical models (e.g., Miotello et al. 2014, 2016). It has been shown that in T Tauri disks the CO emission is much weaker than expected, which could be due to carbon- and oxygen-rich volatiles being locked up as ice resulting in chemical conversion into more complex ices (Bosman et al. 2018; Agúndez et al. 2018) and radial transport (Krijt et al. 2018), both of which lead to low CO fluxes and making CO less ideal to trace the bulk mass (see e.g., Ansdell et al. 2016, 2017; Pascucci et al. 2016; Miotello et al. 2017). Moreover, if CO is not abundant enough, it cannot self-shield, resulting in a drop in CO abundance due to UV photons dissociating the molecule (Visser et al. 2009). Alternatively, Miotello et al. (2023) show that the lack of CO emission can also be explained by compact gas disks for those disks that remain unresolved in ALMA observations. Because of the more luminous star, Herbig disks are expected to be warmer and thus less CO conversion should occur (Bosman et al. 2018). HD upper limits combined with dust based masses from Kama et al. (2020) do indeed imply that this is likely the case. Moreover, using literature values of gas-to-dust ratios Miotello et al. (2023) found that the CO depletion in Herbig disks is lower than for T Tauri disks. However, a general survey of Herbig gas masses is still lacking.

In addition to the gas mass, CO has also been used to trace the radius of the disk (Ansdell et al. 2017; Trapman et al. 2019). The outer disk radius can be used to distinguish between the two main evolutionary scenarios proposed for disks: viscous evolution and wind driven evolution (see for a recent overview Manara et al. 2023). The former results in an increasing outer radius with evolution, while the latter decreases the outer radius of the disk. Modeling has shown that this is also visible in CO emission (Trapman et al. 2022). In the Herbig disk HD 163296, a disk wind has been detected in ¹²CO and ¹³CO (Klaassen et al. 2013; Booth et al. 2021), suggesting a possible wind driven scenario operating in that Herbig disk. However, clear conclusions on which of the two scenarios is the main driver of evolution in disks is still missing. A survey is lacking here as well.

This work builds upon the previous work by Stapper et al. (2022) who obtained the dust masses of a comprehensive survey of Herbig disks observed with ALMA. We use their sample in addition to five observations done with the Northern Extended Millimeter Array (NOEMA) to obtain information on all Herbig disks with millimeter CO isotopologue observations. These CO observations are used to determine the gas masses in the same manner as Miotello et al. (2016), as well as gas disk radii. The obtained masses and radii are used to answer important questions in the field. The gas masses give us insights on whether there is less CO conversion in Herbig disks compared to T Tauri disks, and if CO can be used as an effective gas mass tracer in Herbig disks. The gas and dust radii can be used to investigate the efficiency of radial drift in Herbig disks, and the CO observations can be used to determine if there are any more wind signatures present in the available Herbig disk ALMA data. Also, comparisons to T Tauri and debris disks can be made.

This work is structured as follows. Section 2 details the selection of the sample and accompanying data reduction. Section 3 describes the DALI model grid parameters. The results are presented in Sect. 4, where the integrated-intensity maps of ¹²CO, ¹³CO and C¹⁸O are presented in Sect. 4.1, the ¹²CO radii are compared to the dust radii in Sect. 4.2, the luminosities of the observations and models are compared in Sect. 4.3 and the resulting gas mass estimates are shown in Sect. 4.4. We discuss these results in Sect. 5, where in Sect. 5.1.1 we investigate how to unravel the different masses of the disks using rare isotopologues, and introduce a new technique to compute the full mass of a disk using CO isotopologues in Sect. 5.1.2. In Sect. 5.2 we compare the obtained masses to those of T Tauri disks in Sect. 5.2.1, between the group I and group II Herbig disks in Sect. 5.2.2, and to those of debris disks in Sect. 5.2.3. Lastly, in Sect. 5.3 we present results into finding disk winds in the data, and discuss the wind driven versus viscously driven evolution of Herbig disks using the obtained disk radii. Section 6 summarizes our conclusions.

2 Sample and data reduction

The sample used in this work is compiled by Stapper et al. (2022), who selected all Herbig disks from Vioque et al. (2018) with ALMA data within 450 pc. In this work, we update this sample by using the Gaia DR3 updated stellar parameters and distances from Guzmán-Díaz et al. (2021). We include HD 34282, KK Oph, V892 Tau and Z CMa from Vioque et al. (2018) as well, which are not included in Guzmán-Díaz et al. (2021).

A selection of the available data on the ALMA archive¹ was made which cover the ¹²CO, ¹³CO, and/or C¹⁸O lines. In general, the J = 2−1 transition was chosen, except in cases where only the J = 3−2 transition was available. The data were chosen such that the largest angular scale is at least as large as the continuum disk size, as to not filter out larger scales. In cases where high resolution data were used (either due to a better sensitivity or no low resolution data were available), a uv-taper was applied, to taper the resolution to ~0.1″ to easier be able to determine the radius and integrated flux.

These selection criteria resulted in 30 disks for which at least one ¹²CO, ¹³CO or C¹⁸O transition is available. R CrA is not included because it only has ACA (Atacama Compact Array) data available, which is prone to having contamination from nearby sources due to the much lower spatial resolution. Out of the 30 datasets, two (VV Ser and HD 245185) only have ¹²CO observations available.

In addition to the ALMA archival data, we present observations made with the Northern Extended Millimeter Array (NOEMA) toward five Herbig disks which lie too far north for ALMA. These disks (BH Cep, BO Cep, HD 200775, SV Cep and XY Per) all lie at distances similar to Orion, ranging from 324 pc to 419 pc (Guzmán-Díaz et al. 2021). These data were calibrated using the IRAM facilities. Details on both the NOEMA and ALMA data and their project codes can be found in Table A.1.

The data were imaged using the Common Astronomy Software Applications (CASA) version 5.8.0 (McMullin et al. 2007). The data were first self-calibrated using the continuum, which was made by combining all spectral windows and flagging the channels containing line emission. For the ALMA data, up to six rounds of phase-only self-calibration were computed, starting with a solution interval equal to the duration of a single scan. The solution intervals were shortened by a factor of two each round. For the NOEMA data, only one phase calibration round was used. After phase-only self-calibration, a single round of amplitude self-calibration was performed. Both the phase and amplitude calibrations were only applied if an increase of more than 2% in peak signal-to-noise ratio (S/N) was found, for which the rms was determined by excluding the disk emission and using the full field-of-view of the observation. Typically a factor of 1.8 was gained in peak S/N. Afterwards, the resulting calibration table was applied to the line spectral windows by using the applycal task in CASA. See Richards et al. (2022) for more information on self-calibration. After self-calibration each dataset was continuum subtracted using the uvcontsub task in CASA using a linear fit. For masking during the CLEAN process, hand-drawn masks were aided by the auto-multithres² option in tclean to automatically generate masks. For the resolved data, the data were cleaned using the multiscale algorithm in tclean. The used scales were 0 (point source), 1, 2, 5, 10 and 15 times the size of the beam in pixels (~5 pixels). The last three scales were only used if the disk morphology allowed for it. All unresolved data (which is, in addition to the NOEMA data, some of the ALMA data as well) were cleaned using the hogbom algorithm. For the ALMA data, every dataset was imaged using a Briggs robust weighting of 2.0, optimizing the S/N of the resulting image. For the NOEMA data, a Briggs robust weighting of 0.5 was used, due to its lower spatial resolution, without impacting the S/N. Finally, the resulting image cubes were primary beam corrected. All resulting image parameters can be found in Table A.1 which also lists the root-mean-square noise of the data. Additionally, the NOEMA continuum data are presented in Appendix B.

The velocity integrated (moment 0) maps of the image cubes were obtained by using a Keplerian mask. We use the same implementation as Trapman et al. (2020; also see Salinas et al. 2017; Ansdell et al. 2018). The used inclinations, position angles, system velocities and internal velocities (i.e., a range in velocity to take into account that emission is not geometrically thin and comes from different layers in the disk, extending the velocity range over which the mask is generated, see Appendix A in Trapman et al. 2020) can be found in Table 1. The internal velocity is for the majority of the disks set at 1.5 km s⁻¹, and generally provides an extra factor to improve the coverage of all emission by the Keplerian mask. The distances and stellar masses are taken from Guzmán-Díaz et al. (2021) or Vioque et al. (2018), obtained via evolutionary tracks on the HR-diagram, except for V892 Tau, for which we use the mass estimate from Long et al. (2021), because both Guzmán-Díaz et al. (2021) and Vioque et al. (2018) do not include a mass estimate. We note that other stars have dynamically estimates of their mass as well (e.g., AK Sco, Czekala et al. 2015; HD 169142, Yu et al. 2021; HD 163296, Teague et al. 2021; HD 34282 Law et al. 2023), but we do not include these to keep the stellar masses homogeneously derived for the complete sample. As all emission is included in the Keplerian masking, the stellar masses are not expected to impact our results. The references for the inclinations and position angles as obtained from continuum observations are given in Table 1. The position angles have been altered by hand by small amounts such that the Keplerian masks fit the channel maps well. By applying the Keplerian mask the moment 0 maps are made, see Fig. 1. The integrated fluxes were determined using a curve-of-growth method, see Stapper et al. (2022) for more details. The final luminosity of each line was computed by multiplying the integrated flux by a factor 4πd² where d is the distance to the source in parsecs.

For the non-detections, the upper limits are estimated by using three times the noise in a single channel and multiplying this by the width of the channel and the square root of the number of independent measurements, assuming that the emission is coming from a range of 10 km s⁻¹ in velocity and within a single beam. We note that this was also done for the binary XY Per (see Appendix B). Similarly, the error on the integrated fluxes was estimated by multiplying the noise in an empty channel by the width of the channel and the square root of the number of independent measurements in the Keplerian mask and the area over which the integrated fluxes were computed.

In addition to the integrated flux, the curve-of-growth method also gives the size of the disk. We compute the radii encircling 68% and 90% of the flux. We use a minimum error on the radius of 1/5 the size of the beam (i.e., the pixel size). Similarly, the errors on the dust radius from Stapper et al. (2022) are taken as 1/5 the size of the beam. Upper limits on the size of the disk are given if, after fitting a Gaussian to it with the CASA task imfit, the major or minor axes of the Gaussian are smaller than two times the beam. We note that for BO Cep, HD 104237 and HD 245185, while marginally resolved with only two beams along the minor axis, we estimate an inclination and position angle to cover all emission with the Keplerian mask. See the resulting values in Table 1. For the other unresolved disks, we either use the inclination and position angle of previous works (HD 9672, HD 142666, HD 139614) for the Keplerian mask or we use no Keplerian mask (V718 Sco). We include the inferred radii for these seven disks as upper limits in Table 2. Due to the many different projects imaged, the sample has a rather inhomogeneous distribution of spatial scales and sensitivities (see Table A.1). All ¹³CO and C¹⁸O observations are from the same datasets, which is also the case for most ¹²CO observations. The exceptions are: AB Aur, HD 135344B, HD 141569, HD 142666 and HD 290764, which have a different dataset for their ¹²CO observations.

The fluxes and radii together with their errors are listed in Table 2. A summary of Table 2 of all luminosities and radii is shown in Fig. 2.

3 Model setup

In this work a grid of DALI models (Dust And LInes, Bruderer et al. 2012; Bruderer 2013) is run to determine the gas masses of the Herbig disks based on the ¹³CO and C¹⁸O luminosities. DALI is a thermo-chemical code which solves for the gas and dust thermal structure of the disk by taking heating, cooling, and chemical processes into account. We use the CO isotopologue chemistry network by Miotello et al. (2016) evolved to 1 Myr, which is a simplified version of the network by Miotello et al. (2014). This network includes the ¹²C, ¹³C, ¹⁶O, ¹⁸O, and ¹⁷O isotopologues, and processes such as isotope-selective photodissociation, fractionation reactions, self-shielding, and freeze-out.

DALI uses the simple parametric prescription proposed by Andrews et al. (2011) for the density structure. This density structure is motivated by the viscous accretion disk model, for which the solution follows an exponentially tapered power law (Lynden-Bell & Pringle 1974; Hartmann et al. 1998) given by $${\text{Σ}}_{\text{gas}}={\text{Σ}}_c{\left(\frac{R}{R_c}\right)}^{-\gamma }\exp \left[-{\left(\frac{R}{R_c}\right)}^{2-\gamma }\right],$$(1)

where Σ_c and R_c are respectively the surface density and critical radius and γ the power law index. The vertical distribution of the gas is given by a Gaussian distribution. The scale height angle of the gas is given by $$h=h_c{\left(\frac{R}{R_c}\right)}^{\psi }\text{,}$$(2)

where ψ is the flaring index and h_c is the scale height at distance R_c. The physical scale height H is then equal to hR.

The models include small and large dust populations, for which the small dust follows the gas distribution and the large dust is settled vertically. The dust settling is set by the settling parameter χ. The gas-to-dust mass ratio in the disk is given by ∆_g/d and is set to the ISM value of 100.

In this work we use a range of parameters to run a grid of DALI models, see Table 3. All of these parameters combined make a total of 3600 models. Besides varying the gas mass of the models from 10⁻⁵ M_⊙ to 10^−0.5 M_⊙ in steps of 0.5 dex, we also vary the parameters related to the vertical and radial mass distribution of the gas following Eqs. (1) and (2). We expand the range of values used in Miotello et al. (2016) by also including flat and compact disks, as these are likely present in the group II Herbig disks (Stapper et al. 2023), through h_c, ψ, and R_c. Therefore, the critical radius ranges from 5 to 200 au in steps of ~0.4 dex. The power-law index γ is set to 0.4, 0.8, and 1.5, which changes the steepness of the turnover in the surface density at large radii. Higher values of γ allow for the bulk of the mass to be distributed further out in the disk, effectively changing the size of the disk.

For the vertical distribution of the gas, we alter both the flaring index ψ and the scale height at the critical radius h_c.

The range in ψ is kept the same as Miotello et al. (2016) used, see Table 3. For the scale height h_c we use a larger range from a very flat disk with h_c = 0.05 rad, as very flat disks have been found to exist in the Herbig population (Law et al. 2021, 2022; Stapper et al. 2023), up to 0.2 rad and 0.4 rad for thicker disks. Flaring of the disk changes the (vertical) temperature structure of the disk and hence the emitting layer of the CO, which in turn changes the observed luminosity of the line.

To better accommodate for the range in stellar properties, we vary the stellar luminosity from 5 L_⊙ to 50 L_⊙ in steps of ~0.3 dex to include a similar range in luminosities present in our Herbig sample (not including the most luminous stars), as determined by Guzmán-Díaz et al. (2021). We keep the effective temperature T_eff at 10⁴ K as this matches well with the effective temperatures of the sample, and we use the sublimation radius R_subl corresponding to a stellar luminosity of 10 L_⊙ (using the scaling relation by Dullemond et al. 2001) for all models, as the sublimation radius is much smaller than the scales we probe.

We kept the remaining parameters the same as Miotello et al. (2016). The PAH abundance is set to 1% of the ISM, with the ISM value being a PAH-to-dust mass ratio of 5% (Draine & Li 2007). The cosmic-ray ionization rate (ζ_cr) is set to 5.0 × 10⁻¹⁷ s⁻¹, and the X-ray luminosity (L_X) is set to 10³⁰ erg s⁻¹. As noted in both Bruderer et al. (2012) and Miotello et al. (2016), L_X only has a minor influence on the CO pure rotational lines, which are the focus in this work.

The models have been run with 40 vertical cells and 42 log-spaced radial cells out to 1000 au. The highest mass models were radially sampled in 62 cells out to 1500 au. This sampling was found to be sufficient to yield line fluxes accurate to <5%.

The grid extends vertically to ten times the scale height, as given by Eq. (2).

For the built-in DALI ray-tracer, a distance of 100 pc is used and each model is ray-traced at inclinations of 10°, 40° and 70°, encompassing all inclinations in our sample (see Table 1). A total of six CO isotopologues are ray-traced: ¹²CO, ¹³CO, C¹⁸O, C¹⁷O, ¹³C¹⁸O, and ¹³C¹⁷O. The J = 2−1 and J = 3−2 transitions are ray-traced for each molecule (Yang et al. 2010; Schöier et al. 2005).

Fig. 1 Keplerian masked velocity integrated-intensity maps of the 23 Herbig disks with at least one detection of 12CO, 13CO, or C18O. For each disk, 12CO, 13CO, and C18O are given in the first, second, and third panel from left to right. If nothing is shown, the molecule is not covered. The bar in every first panel from the left indicates a size of 100 au with the corresponding angular size below. The size of the beam is indicated in the bottom left of each panel. A sinh stretch is used to make the fainter parts of the maps better visible. The minimum value is set to zero, and the maximum value in mJy beam−1 km s−1 is indicted on each panel in the top left corner. The disks with no detections in all three isotopologues are not shown here.

4 Results

4.1 Integrated-intensity maps

Figure 1 presents the integrated intensity maps of the disks in which any of the three CO isotopologues is detected: 20 out of 27 for ¹²CO, 22 out of 33 for ¹³CO, and 21 out of 33 for C¹⁸O for our 35 sources. If nothing is shown, the molecule is not covered. As C¹⁸O is covered but not detected, the integrated-intensity map of HD 141569 is presented as well, which is obtained by integrating over the same velocity range as done for the ¹²CO and ¹³CO emission. All other non-detections are not shown. For the majority of the disks, the transition used is J = 2−1, except for all lines in the HD 135344B and HD 290764 disks, and ¹²CO in AB Aur and HD 141569, for which the J = 3−2 transition is imaged. For the detections, azimuthally averaged radial profiles are presented in Fig. 3 which have been made using the inclinations and position angles listed in Table 1. Figure 3 also presents the radii listed in Table 2.

In general we find that the size of the disk decreases with the abundance of the isotopologue. For the resolved disks the median ¹²CO R_68% (R₉₀%) radius is a factor of 1.3 (1.3) larger compared to ¹³CO, and a factor of 1.8 (1.4) larger compared to C¹⁸O. Hence, for C¹⁸O the bulk of the emission is generally compact, but it can have more faint extended emission around the compact emission compared to ¹²CO and ¹³CO, resulting in a larger difference between the R_68% and R_90% radii. The decrease with abundance of the isotopologues observed is caused by the smaller column density of the rarer CO isotopologue and isotope selective photodissociation, which lead to optically thin millimeter emission at larger radii reducing the amount of emission. Regarding the luminosity of the disk, we find a clearly decreasing trend with rarity of the isotopologue. Most ¹²CO disks have a luminosity of at least 5 × 10⁶ Jy km s⁻¹ pc², while the ¹³CO and C¹⁸O disk luminosities are 2 × 10⁶ Jy km s⁻¹ pc² and 4 × 10⁵ Jy km s⁻¹ pc² respectively. As all three isotopologues are likely optically thick in our sample when detected, this decrease in luminosity is a result of a decrease in size of the detected disk, a reduction in temperature due to the emission from rarer isotopologues originating from deeper in the disk, and/or truncation by photodissociation.

The sizes of the ¹²CO disks are seen to vary significantly. The smallest detected gas disks range from an R_90% of less than 64 au (HD 104237) out to 82 au (HD 100453). The HD 104237 disk is particularly small both in CO and continuum: we measure a continuum disk extent of less than 21 au, which is still unresolved (we use a higher spatial resolution dataset than Stapper et al. 2022, who had an upper limit of 139 au). Most of the ¹²CO disks have a R_90% between 100–400 au (65%, 13/20; see Fig. 2). The median disk size lies close to the middle of this range at 335 au. The R_68% radius has a median of 226 au. The largest disk in ¹²CO is HD 142527 which has a R_90% of 793 au, by far the largest disk in the sample. AB Aur is likely a very large disk in ¹²CO as well, because of its ¹³CO radius of 1118 au. But due to a small maximum recoverable scale a large part of the ¹²CO disk is filtered out. A number of disks show large extended emission especially in ¹²CO (but also for ¹³CO and C¹⁸O) where there is an initial peak and a much flatter emission “shoulder” toward larger radii. This part contributes to the integrated flux only slightly, and mainly increases the inferred radius of the disk, especially for the 90% radius. Hence, R_90% is the main tracer of the outer radius, and we primarily use this over the 68% radius, except when explicitly stated otherwise.

The ¹³CO and C¹⁸O sizes are generally smaller than the ¹²CO sizes. As Fig. 2 shows, most of the ¹³CO and C¹⁸O emission is within 300 au in size. The range of the inferred radii is larger than for ¹²CO, ranging from the smallest resolved ¹³CO disk in the sample of 61 au, out to 1118 au in size for the largest disk. The most stringent upper limit obtained is similar as the smallest resolved disk, at 59 au, for HD 104237. There are also more upper limits on the ¹³CO and C¹⁸O radii compared to ¹²CO due to four relatively low resolution datasets used not containing ¹²CO.

Some disks suffer from foreground cloud absorption, such as HD 97048 and V892 Tau, as can be seen in Fig. 1 where emission along the minor axes of these disks seems to be absent. Especially the ¹²CO emission of AB Aur and some of HD 97048 suffers from this, resulting in a smaller outer radius found for ¹²CO compared to the other isotopologue(s). For ¹³CO and C¹⁸O the cloud absorption leaves much less of an imprint on the moment map, and these do not affect the obtained integrated fluxes and the radii significantly. Hence, we mainly use the ¹³CO radius as a measure of the size of the disk. Several of the non-detections show foreground cloud emission: HD 176386, MWC 297, TY CrA and Z CMa. For Z CMa emission from the disk is present at V_LSRK ~ 10 km s⁻¹, but the foreground cloud emission makes it difficult to extract any information from this, and thus it is considered as an upper limit. For the other disks with foreground cloud emission no signature of a disk is visible in any of the isotopologues, possibly due to the cloud emission.

Due to the different datasets, their spatial resolution is different resulting in unresolved ¹³CO and C¹⁸O disks for HD 142666, and different sensitivities resulting in upper limits larger than the ¹²CO disk. This is especially clear from the azimuthal aver ages in Fig. 3, where one can see that the inferred radii of ¹²CO are indeed smaller than for the rarer isotopologues in these two disks. Lastly, for BO Cep, HD 104237, HD 142666, and V718 Sco we do find very similar outer radii for ¹³CO and C¹⁸O because the disk is unresolved. The other two disks with upper limits on their size, HD 139614 and HD 9672, are marginally resolved, especially along the major axes, resulting in a more significant difference between the ¹³CO and C¹⁸O radii.

Fig. 2 Summary of the luminosities and 90% radii for the 12CO, 13CO, and C18O observations as presented in Table 2.

Fig. 3 Azimuthally averaged radial profiles of the three CO isotopologues in all 23 Herbig disks in which at least one isotopologue is detected. The 1σ uncertainty interval is indicated by the shaded region in the same color as the profile. The vertical solid lines indicate the derived 90% radii for each line. An upper limit on the radius is indicated with a left facing arrow. The beam size is illustrated by the horizontal lines in the top right of each panel.

Fig. 4 Gas radii versus dust radii of the Herbig disks with detected 12CO emission, both R90%. An upper limit on the gas radius or dust radius, or both, is marked with a triangle pointing in the corresponding direction(s). The colors indicate the Meeus et al. (2001) group of the disk: red is group I, blue group II. Lines for different Rdust to Rgas ratios are shown as well. The black line shows for the resolved disks the relation of the R90% radii. Using R68% does not alter the relation significantly. The yellow region indicates the region where the difference between the dust and gas radii cannot solely be explained by optical depth effects, and radial drift is necessary (Trapman et al. 2019). The red line corresponds to the relation found for the disks in Lupus (Ansdell et al. 2018).

4.2 ¹²CO radius versus dust radius

The ¹²CO observations can be used to compare the gas radii of the Herbig disks with the dust radii from Stapper et al. (2022). Figure 4 presents the dust radius compared to the gas radius. In addition to HD 104237, a higher resolution dataset of HD 141569 has been used as well compared to Stapper et al. (2022), which further constrains the dust radius to less than 38 au for the 68% radius and 54 au for the 90% radius. In addition to the region in the plot which is radial drift dominated, the relation found for the disks in Lupus is also indicated (Ansdell et al. 2018). If the gas radius is larger, the disk is thought to be radial drift dominated (Trapman et al. 2019). Anything below this value can still be caused by optical depth effects.

All disks but one (HD 9672) have larger gas radii compared to their dust radii. Most disks lie along the same relation as found for Lupus, where the gas is two times the dust radius. However, a larger spread in ratios is present in the Herbig sample compared to the disks in Lupus. The black line in Fig. 4 indicates the relation for the resolved disks between the dust and gas radii for the R_90% radii. A simple fitting routine curve_fit is used from the scipy package to fit the scaling between the dust and gas radii. We find that for the R_90% radii, the ratio between the dust and gas radii is a factor of 2.7. The R_68% are consistent with this value. HD 9672 is the outlier in this sample, showing an equal size in both the dust and gas. The gas radius is slightly unresolved however, and this could make the gas radius even smaller. A comparison with debris disks is shown in Sect. 5.2.3.

These ratios are high compared to Lupus due to a few disks which are close to or within the regime where radial drift is necessary to explain the differences seen between the dust and gas radii (Trapman et al. 2019). This is the case for HD 100546, HD 142666, HD 163296, HD 31648 (MWC 480) and V892 Tau, all of which have gas disks a factor of more than four larger than their dust radius. This is consistent with other works of for instance MAPS (Zhang et al. 2021). A few disks have gas radii only slightly smaller than four times the dust radius, HD 142527 and HD 245185. We note that HD 100546 has a very faint outer dust ring at 190 au (Walsh et al. 2014; Fedele et al. 2021), lowering the ratio to 1.7 instead of almost a factor of 8. Similarly, HD 163296 has a faint outer ring in both the DSHARP and MAPS data as well (Huang et al. 2018; Sierra et al. 2021). Hence, Fig. 4 could indicate around which disks a faint outer ring can be found, possibly linked to giant exoplanet formation happening in these disks.

4.3 ¹³CO and C¹⁸O luminosities

Figure 5 presents the resulting ¹³CO and C¹⁸O J = 2−1 line luminosities of the models in colors and observations in black (the same figure for the J = 3−2 transition can be found in Appendix C). The size of the markers are scaled by the ¹³CO 90% radius of the models and observations. The disks with upper limits on both ¹³CO and C¹⁸O are presented as the diagonal triangles. HD 141569, which has a C¹⁸O upper limit, is shown as a downward pointing triangle. The top and right hand-side panels indicate the distribution of the models for the ¹³CO and C¹⁸O luminosity respectively.

The two main deciding factors governing how luminous the disk is are the optical depth of the millimeter lines as seen from the observer and the self-shielding capacity of the CO isotopologues against UV (e.g., Visser et al. 2009). Additionally, these two factors depend on the distribution of the mass and the size of the disk (Trapman et al. 2019, 2020). If the disk is very compact a large range of masses can be ‘hidden’ behind the region where the gas becomes optically thick (Miotello et al. 2021). This results in a large degeneracy in possible luminosities for a single mass, as can be seen in Fig. 5. The only way to change the observed luminosity is to increase the emitting area as for optically thick emission, the intensity per unit area is constant, and thus increasing the radius of the disk. So a clear trend is present in the models where the smallest disks are at low ¹³CO and C¹⁸O luminosities, while the largest disks are on the opposite side with high ¹³CO and C¹⁸O luminosities due to an increase in radius.

Figure 6 presents this in another way: for the most massive disks an increase in radius results in an increase in luminosity. In addition, Appendix D shows this for the other values of R_c and γ. This trend is also evident in the observations, where the largest disks such as AB Aur and HD 97048 are in the top right of Fig. 5 and the smallest disks such as HD 100453 and HD 104237 are in the bottom left. We can therefore conclude that most disks, apart from HD 9672 and HD 141569, are likely to be optically thick in both ¹³CO and C¹⁸O as for these disks both luminosities scale with the size of the disk.

The models with higher mass disks can also make larger disks due to the disk surface density being higher at larger radii. This results in a maximum ¹³CO and C¹⁸O luminosity for both the 10⁻¹ M_⊙ and 10⁻² M_⊙ models. The distributions presented in the top and right panels of Fig. 5 show that the 10⁻¹ M_⊙ models can have luminosities up to ~4 × 10⁷ Jy km s⁻¹ pc² for ¹³CO and ~10⁷ Jy km s⁻¹ pc² for C¹⁸O. Similarly, the 10⁻² M_⊙ models have luminosities up to ~2 × 10⁷ and ~2 × 10⁶ Jy km s⁻¹ pc² for ¹³CO and C¹⁸O, resulting in a region where the mass of the disks needs to be at least 0.1 M_⊙ to explain the observed luminosities.

For the lower mass disks, an interplay between the optical depth and self-shielding of the CO molecules sets the observed luminosities. Due to its larger abundance, ¹³CO becomes optically thin further out in the disk than C¹⁸O. For the 10⁻² M_⊙ models, first an increase in radius results in an increase in luminosity in both isotopologues, but eventually C¹⁸O becomes optically thin. Consequently, an increase in radius does not increase the luminosity of C¹⁸O, but it still does so for ¹³CO, and the model moves horizontally with an increase in radius, also see the bottom left panel of Fig. 6. When the C¹⁸O is spread out even more, its self-shielding ability decreases and the C¹⁸O starts to photodissociate, decreasing its abundance and thus its luminosity. This results in the model moving down with an increase in radius in Fig. 6. Consequently, a maximum luminosity for C¹⁸O is set by photodissociation for disk masses ≲ 10⁻³ M_⊙. ¹³CO is on the other hand still optically thick, and increasing the radius still increases its luminosity, until ¹³CO starts dissociating as well, causing the models to curve toward lower luminosities in both isotopologues with increasing radius for M_disk ≲ 10⁻⁴ M_⊙.

Photodissociation of C¹⁸O is already the main driver for the luminosity of the lower mass models of almost all sizes, hence almost all models move horizontal or diagonally down to lower ¹³CO luminosities. For the smallest models, while still being optically thin (or marginally optically thick) in C¹⁸O, an increase in size can still result in an increase in luminosity if the photodissociation of the C¹⁸O does not yet dominate. Also in the M_disk ≲ 10⁻⁴ M_⊙ case the ¹³CO becomes optically thin for the observer and photodissocation starts to lower its luminosity with an increase in radius. Going to the lowest mass models of 10⁻⁵ M_⊙ both ¹³CO and C¹⁸O are (marginally) optically thin and photodissociation of both isotopologues reduce their luminosity with increasing radius.

Evidently, the mass and the radius of the disk have the biggest impact on the observables of the disk. The next most impactful parameter affecting the CO luminosity of the disk is γ in Eq. (1), which affects the overall distribution of the mass of the disk. A larger γ results in a higher surface density at the inner region and decreases the effect of the exponential taper in the outer region, which increases the size of the disk. γ especially impacts the lower mass models, see Fig. 6 and the extra figures in Appendix D. A larger γ makes the region where the C¹⁸O emission is coming from relatively constant in size and self-shielding occurs less than for smaller γ even for the largest R_c, while the ¹³CO emission increases with increasing radius because of a larger emitting surface. Hence, the model moves horizontally in the bottom left panel of Fig. D.6. For smaller γ at C¹⁸O can dissociate at large radii due to the lower surface density closer in and the luminosity decreases rapidly, as can be seen as the relatively large jump in luminosity from R_c = 60 au to 200 au for the 10⁻³ M_⊙ models in the bottom left panel of Fig. 6. This jump also causes the bimodal distribution present in the right panel of Fig. 5. This could have been alleviated by adding a value of γ between 0.8 and 1.5. For even lower masses an increase in γ results in less change in both ¹³CO and C¹⁸O with an increase in radius due to more efficient self-shielding.

All remaining parameters affect the CO luminosity in the same direction as R_c, see Figs. 6, and D.1 to D.6. Changing h_c does not affect the more massive disks, because these are optically thick already. For the lower mass models increasing h_c has the main affect of increasing the volume of the disk, making the gas optically thin and reducing the abundance of both ¹³CO and C¹⁸O due to photodissociation. Increasing the stellar luminosity decreases the abundance of both ¹³CO and C¹⁸O, hence the downward left movement of the models in Fig. 6. For the higher mass disks the CO isotopologues are already shielded and the luminosity barely affects the overall abundance. Lastly, both ψ and the inclination do not affect the observed integrated ¹³CO and C¹⁸O luminosities considerably.

Fig. 5 C18O luminosity versus 13CO luminosity for the DALI models (colored circles) and the observations (black circles). The colors indicate the gas mass of the disk model. The size of the markers indicate the 90% radius of the 13CO emission. Probability density curves of the models for each gas mass are shown along the x-axis in the top panel and along the y-axis in the right panel, i.e., most models of a particular gas mass reside around the peak of each curve.

Fig. 6 Overview of the different parameters explored in the DALI models and their effect on the luminosity of 13CO and C18O. The colors indicate the mass of the model. The arrows indicate in which direction the parameter increases in value (see Table 3). The black outlined circles are the same models for each mass in all panels. The black outlined models have the following parameters: γ = 0.8, ψ = 0.2, hc = 0.2 rad, Rc = 200 au, i = 10° and L★ = 10 L⊙. In Appendix D the same plot is shown, but with γ = 0.4 and γ = 1.5, and with Rc = 5 au, Rc = 10 au, Rc = 30 au, and Rc = 60 au, and for the J = 3−2 transition.

4.4 Obtaining the mass of the disk

As discussed in the previous section, Fig. 5 shows that most observations are consistent with the highest possible masses based on their ¹³CO and C¹⁸O luminosities. However, at least for the observations detecting ¹²CO, ¹³CO, and C¹⁸O, the mass range can be further constrained by using the observed size of the disk in the selected isotopologues.

4.4.1 Masses from detected emission

To obtain a measure of the mass of the disk, we select all models within a region around the observations in Fig. 5 based on the confidence intervals of the observations. There are two main sources of uncertainty that we take into account: a systematic uncertainty set by the calibration accuracy, and a statistical uncertainty set by the noise in the moment 0 map. We select all models which fall within 3σ from a line given by the 10% systematical uncertainty. For some of the more luminous targets the obtained uncertainty is relatively small and either no or very few models were found within the region confined by the two sources of uncertainty. If less than 200 models were found to fall within the given region of a disk, the closest 200 models were included, this occurred for 11 disks. Typically, the selected models deviated around 10% from the observed luminosity. For HD 141569 we take all models below the C¹⁸O upper limit and within a 10% systematic uncertainty on the ¹³CO luminosity because it dominates over the noise.

After obtaining this set of models, selections based on the size of the disk and its inclination, and on the stellar luminosity can be made (see Tables 1 and 2). Given the lack of ¹²CO observations for some of these disks, together with possible cloud contamination, we use the ¹³CO R₉₀% as a disk size tracer. We select all models within a factor of 1.4 in size compared to the observations, as this is the factor by which the models are on average smaller than the observations in ¹³CO, we discuss this further at the and of this section. Based on the radius, inclination, and stellar luminosity, the constraints on the disk mass can be tightened. Figure 7 presents the resulting range in possible disk mass values and the (logarithmic) mean mass of the models within that range, as shown by the light gray lines and darkblue diamonds respectively. The spread in mass of the selected models is given by the black lines. For the high-mass disks (on the right), that correspond to bright ¹³CO emission, the masses are well constrained to be above 10⁻¹⁵ M_⊙. These stringent lower limits are mainly due to their large and well-known sizes. Our obtained disk masses can be found in Appendix E.

Figure 8 shows that the C¹⁸O luminosity is sensitive to the gas mass when selecting models with the appropriate size based on our grid of models, each panel showing this for differently sized models in ¹²CO. We note that this is the R₉₀% radius, which should not be confused with R_c. This clearly shows that the size of the disk is an indicator of its mass as well. One cannot make disks larger than ~500 au without the disk mass being higher than ~10⁻² M_⊙. On the other hand, for the lower mass disks the tenuous ¹²CO gas cannot self-shield in the outer regions anymore for the largest disks, reducing the measured size of the disk. Consequently, one obtains a better constraint on the disk mass of the largest disks in the sample.

For the smaller disks, more masses are compatible with the observed ¹³CO and C¹⁸O line luminosity due to the higher optical depth. This results in low lower limits on the gas mass from our models. For example AK Sco and HD 104237 have lower limits of ~10^−3.5 M_⊙ and ~10⁻³ M_⊙ respectively. This is also evident from Fig. 8, clearly showing that the different gas masses give similar C¹⁸O luminosities for the smallest disks. For the smallest radius bin size a lower limit on the size of the disk can be given as well. Disks which are smaller than R_90% = 20 au in ¹²CO and have a mass of more than 10⁻⁴ M_⊙ do not occur in our models.

To constrain the disk masses even more, one can include other parameters as well. For example, for AK Sco and HD 142666 the vertical extent of the disk from Stapper et al. (2023) can be exploited. We implement this in a simple way, where for the very flat disks AK Sco and HD 142666 we choose models with h_c = 0.05 rad. While this does not improve the overall range in possible disk masses (as h_c does not considerably change the overall luminosity, see Sect. 4.3), for both flat disks we can rule out either lower mass models for AK Sco or higher mass models for HD 142666. This difference is due to their relative position in Fig. 5.

For two low mass disks, HD 9672 and HD 141569, we can determine their disk mass within an order of magnitude. Based on their position in Fig. 5 (the two disks with the lowest C¹⁸O luminosity or upper limits thereon), we can already see that the gas mass should be around 10⁻³–10⁻⁴ M_⊙. No other models fall within the computed confidence intervals. We note that for HD 9672 we do not reproduce the size of the ¹³CO disk well with the models. We comment more on this in Sect. 5.2.3. This results in a well determined disk mass of 10⁻⁴ M_⊙ for HD 9672 and 10⁻³ M_⊙ for HD 141569.

Based on the stellar mass of the Herbig stars, and the disk mass and disk outer radius of ¹²CO from the models, we can determine if a specific model would be gravitationally unstable around the Herbig star using the relation $$\frac{M_d}{M_{\star }}>0.06\left(\frac{f}{1}\right){\left(\frac{T}{10 \text{K}}\right)}^{1/2}{\left(\frac{r}{100\text{\hspace{0.17em}\hspace{0.17em}au}}\right)}^{1/2}{\left(\frac{M_{\odot }}{M_{\star }}\right)}^{1/2}$$(3)

from Kratter & Lodato (2016), where T is the temperature of the disk, r the outer radius of the disk (given by the measured ¹³CO R_90% radius) and M_★ the mass of the star, we set the pre-factor f equal to 1. The temperature T is determined with the luminosity of the Herbig star via the relation by Andrews et al. (2013). For seven disks the radius of the disk combined with Eq. (3) is constraining enough to result in a lower mean disk gas mass by ruling out the potentially gravitationally unstable models. For the disks with radius upper limits especially, the compact high mass disks can be ruled out if assuming that these observed disks are not gravitationally unstable.

For eight of the Herbig disks, upper limits on the gas masses have been estimated from HD (Kama et al. 2020), see the downward facing triangles in Fig. 7. In general, the gas masses we find are consistent with the HD upper limits. For HD 163296 and HD 36112, the HD upper limits are roughly equal to our gas mass estimates, suggesting that the observations might have been close to detecting HD. For HD 100453, HD 169142 and HD 97048 the mean gas masses of the models are higher than the HD upper limits, implying that these disks are not as massive as their size and ¹³CO and C¹⁸O luminosities suggest. A combination of modeling the CO emission and HD upper limits could add additional constraints to the obtained gas masses.

Figure 7 also shows 100× the dust mass from Stapper et al. (2022). For the five NOEMA targets we compute the dust masses in the same way as Stapper et al. (2022), see Appendix B. The middle panel in Fig. 7 shows the resulting gas-to-dust ratio from the mean gas mass as the dark red diamonds and from the range in gas masses as the vertical gray line. The horizontal gray line indicates a gas-to-dust ratio of 100. For most disks we find that the total disk mass derived from the dust mass is consistent with a gas-to-dust ratio of 100. Some disks suggest a depletion of dust compared to the interstellar medium. Primarily for AB Aur the 100× dust mass falls well below the gas mass range, giving a gas-to-dust ratio that is two orders of magnitude higher than the canonical value. For some of the other higher mass disks, such as HD 169142, the dust mass also suggests a depletion of dust compared to gas of a factor of a few. This apparent depletion of dust might be related to the dust being optically thick, we comment on this in Sect. 4.4.4. In general however, the dust mass does seem to indeed trace the total disk mass relatively well for Herbig disks, in contrast with the T Tauri disks.

Lastly, we note that the size of the models are generally smaller than those of the observations for the same flux in a particular CO isotopologue. When selecting all models within a factor of two of the observations, we find that the ¹²CO, ¹³CO, and C¹⁸O radii of the models are on average a factor 1.3, 1.4, and 1.9 times smaller than the observed radii. This is likely due to the models being smooth, with no substructures present. The observations do show some gas structures which increase the emitting area of the disk, and also present weak emission extended structures which increase their 90% radius without much affecting the overall flux. Specifically gas cavities, the main gas structure seen in our data, are generally found to be smaller than continuum cavities (van der Marel et al. 2016; Leemker et al. 2022), and thus are likely not affecting the gas as much as it would to the continuum. Comparing our found gas masses to existing disk specific modeling efforts taking into account the structure of the disk show only differences of a factor of a few. For example the MAPS Herbig disks HD 163296 and HD 31648 were found to have gas masses of 0.14 M_⊙ and 0.16 M_⊙ respectively (Zhang et al. 2021), within a factor of a few from our derived masses (other examples include Tilling et al. 2012; Flaherty et al. 2015, 2020). Scaling the observed radii by the aforementioned factors does indeed not affect the inferred range in disk masses for a given disk. However, when using Fig. 8 one should keep this in mind when selecting which panel to use. Implementing gas and dust structures in the models exceeds the scope of this work, but future work could also consider the effect of these structures on the obtained dust and gas masses.

Fig. 7 Gas masses of the Herbig disks detected in the CO isotopologues as selected from the models presented in Fig. 5. The disks are ordered from left to right by increasing 13CO luminosity. The gray lines indicate the range in possible disk mass based on the disk parameters listed in Table 1, while the black lines indicate the standard deviation. The mean mass of these models is given as the dark blue diamond. After removing models which are potentially gravitationally unstable, the mean mass changes into the light blue diamond. The orange diamonds are taken from the dust masses of Stapper et al. (2022) multiplied by the canonical gas-to-dust ratio of 100. The HD upper limits from Kama et al. (2020) are shown as green downward pointing triangles. The C18O lower limits are directly obtained from the integrated flux. The middle panel shows the resulting gas-to-dust ratio based on the mass range and mean values, and the top panel shows the factor by which the mass could be underestimated when using C18O assuming optically thin emission.

Fig. 8 Range of C18O luminosities for different mass disks, selected for 12CO R90% radii from the DALI models presented in Fig. 5. The shaded area indicates the minimum and maximum values, and the solid line is the mean value.

4.4.2 Mass lower limits from C¹⁸O

In many works (see e.g., Hughes et al. 2008; Loomis et al. 2018; Miley et al. 2018; Booth & Ilee 2020), the disk mass has been estimated by using a flux scaling relation, assuming the emission is optically thin. This formula has been used to obtain a lower limit on the gas mass from ¹²CO, ¹³CO, and C¹⁸O flux. As we find that most Herbig disks are optically thick in C¹⁸O, it is useful to obtain a measure of how much the total gas mass is underestimated when using C¹⁸O as a gas mass tracer. The total number of C¹⁸O molecules in the disk from the C¹⁸O flux can be calculated, assuming optically thin emission (see e.g., Loomis et al. 2018), with $$n_{{\text{c}}^{18}\text{o}}=\frac{4\pi }{hc}\frac{F_v\text{Δ}V{d}^2}{A_{ul}{x}_u},$$(4)

where h, c, and A_u1 are the Planck constant, speed of light, and the Einstein A coefficient for spontaneous emission respectively. x_u is the fractional population of the upper level, d the distance to the source, and F_v∆V is the velocity integrated flux over the disk. Using a ratio between CO and H₂ of 2.7 × 10⁻⁴ (Lacy et al. 1994), and ratios of 77 and 560 of ¹²C/¹³C and ¹⁶O/¹⁸O in the local interstellar medium (Wilson & Rood 1994), in combination with a factor of 2.4 for the mean molecular weight, the total disk mass can be calculated. An excitation temperature is necessary to compute the population levels using the Boltzmann equation as well. We assume T_ex = 40 K which is higher than the C¹⁸O brightness temperatures of some of the Herbig disks (see e.g., Zhang et al. 2021), but is the same as other disks in our sample (e.g., HD 100546, HD 135344B, and HD 169142). We use the line properties (e.g., partition function values, A_u1, E_u, and ց_u) from CDMS (Endres et al. 2016). Using the integrated C¹⁸O fluxes from Table 2, we obtain lower limits on the total gas mass of the disk, see Fig. 7. The formula is also applied to our models, a comparison can be found in Fig. 9 and the upper panel in Fig. 7.

Compared to the masses we find from the models, the masses obtained with Eq. (4) are a factor of 10–100 times lower, see Fig. 7. In the top panel of Fig. 9, the blue shaded region indicates the retrieved masses using Eq. (4) with a temperature of 40 K compared to the true mass of the model. This large difference between the two is primarily due to two reasons. First, for the most massive disks, the C¹⁸O emission is optically thick which reduces the obtained disk gas mass. The disks for which the retrieved masses are closest to the true disk mass are also the largest disks, indicating that these are (marginally) optically thin. For the highest mass disks, even the largest disks are not optically thin as the maximum retrieved mass bend toward relatively lower values. For the lowest mass disks on the other hand (e.g., HD 9672 and HD 141569) photodissociation becomes the dominant process reducing the CO abundance. Consequently, the mass of the disk is always underestimated when connecting C¹⁸O directly to the mass of the disk. As can be seen in both Figs. 7 and 9, the lowest mass disks have very low luminosities and retrieved masses due to photodissociation. The gray shaded area in the top panel of Fig. 9 are the retrieved masses after using the gas temperature in the emitting region of C¹⁸O where 90% of the emission is coming from as weighted by the mass in each cell, as indicated in the bottom panel. Both lower mass disks and smaller disks are warmer in the C¹⁸O emitting region, hence using a temperature of 40 K is not adequate. Especially for the lower mass disks this can underestimate the total disk mass with a factor of a few in addition to the lower CO abundance.

Fig. 9 Gas mass of the models versus the retrieved gas mass from the C18O flux using Eq. (4), and the mass-weighted average of the temperature in the C18O emitting region. For the blue shaded region a temperature of 40 K is used, the gray shaded region uses the temperature of each model shown in the bottom panel. The shaded areas are the minimum and maximum values, with the mean value indicated by the solid line. The dashed line indicates the one-to-one correspondence between the two masses. The typical range in temperatures is indicated by the gray shaded region in the bottom panel.

Fig. 10 Same as Fig. 7, but now for the disks with only upper limits on the 13CO and C18O fluxes.

4.4.3 Upper limits

For the disks with upper limits on both isotopologues no selection can be made based on a region of luminosities in Fig. 5 nor the size of the disk can be used. Using the dust radii of Stapper et al. (2022) in combination with the typical ratio of 2.7 between the gas and dust radii (see Sect. 4.2) does not add any constraints either as for all disks only upper limits on the dust disk size are known. Hence, a selection of all models within the quadrant confined by the upper limits is made. Figure 10 presents the resulting upper limits on the gas masses.

Based on the parameters in Table 1, a few constraints on the possible disk mass can be made. This is especially true for the ¹³CO and C¹⁸O upper limits on the left side of Fig. 10, as these upper limits include models of all masses. While for the undetected or unresolved disks the inclination and radius do not give constraints, the luminosity gives enough constraints to lower the upper bound on some of the disks by multiple orders of magnitude, see Fig. 10. We find that HD 58647 and HR 5999 can have at most a mass of 10⁻⁴ M_⊙ to explain the non-detection in both ¹³CO and C¹⁸O. Interestingly, for both of these disks a higher dust mass is found, suggesting an increased abundance of dust compared to the ISM. For HD 176386 and TY CrA an upper limit of 10^−3.5 M_⊙ is found. Z CMa has a higher upper limit of 10^−2.5 M_⊙. The standard deviation, as given by the black lines in Fig. 10, show tighter constraints. For the other upper limits no constraints on the gas mass could be made, even after selecting for the source-specific parameters. However, we do note that the 10⁻¹ M_⊙ models are by far outnumbered by lower mass models as most, but not all, high mass models can be excluded. This is reflected in the mean gas mass and standard deviation, as these are much lower than the maximum mass possible, see the dark blue diamonds and black lines in Fig. 10.

Based on the C¹⁸O upper limits, a maximum radius can be estimated for these disks. As Fig. 8 shows, an increase in size increases the luminosity of the disk. Hence, upper limits on the ¹²CO radius can be determined for the disks with no detections. We find upper limits of 550 au (HD 58647, HD 176386, HR 5999, TY CrA and Z CMa) and 800 au (BH Cep, HD 200775, KK Oph, MWC 297 and SV Cep). Besides non-detections, we also have four disks which have upper limits on their radius, for which we obtain an additional minimum radius to explain the found luminosities, as a decrease in size decreases the luninosity of the disk (see Fig. 8). BO Cep has at least a size of 125 au, and HD 139614 and V718 Sco have a size of at least 40 au. For HD 104237 we obtain the most stringent lower limit of 15 au.

Similarly, for HD 245185 and VV Ser, which only have ¹²CO observations, the mass and radius can also be constrained. ¹²CO is optically thick for lower disk masses, resulting in the same luminosities for multiple orders of magnitude in mass even for the largest disks. These luminosities can be used to obtain a lower limit on the mass and size of HD 245185. To explain the observed ¹²CO emission, the disk needs to be at least 125 au in size and have a mass lower limit of 2 × 10⁻⁴ M_⊙. Due to the non-detection of ¹²CO for VV Ser, the radius can be relatively well constrained to be less than 220 au in size using the ¹²CO equivalent to Fig. 8.

Lastly, Grant et al. (2023) have shown that the relationship between the accretion rate M and dust disk mass is largely flat at ~10⁻⁷ M_⊙ yr⁻¹ over three orders of magnitude in dust mass. Hence, some disks have very short inferred disk lifetimes. Our gas mass estimates do not resolve this problem, as the inferred disk masses from our models are either low (see Fig. 10), or no observations are available.

4.4.4 Cumulative distributions

Figure 11 shows cumulative distributions of the dust masses from Stapper et al. (2022), together with the gas masses as found by the models and computed with Eq. (4). Following Stapper et al. (2022), we obtained the cumulative distributions using the Python package Lifelines (Davidson-Pilon et al. 2021). The shaded area indicates the 1σ confidence intervals. The dark blue distribution is made from the mean values in Figs. 7 and 10. The observed disks for which the range in model disk masses extends down to the lowest mass in our model grid are considered as an upper limit. The dust distribution in orange is obtained by using the dust masses from Stapper et al. (2022) in addition to the five NOEMA targets presented in this work (see Appendix B). The green cumulative distribution is obtained from computing the gas using Eq. (4). The dust and gas distributions show a similar slope for the highest mass disks, indicating a relatively constant overestimate of the gas mass, or a constant underestimate of the dust mass. For lower mass disks, the distribution is set by the non-detections, resulting in a leveling-off at a p ~ 0.7 of the distribution reflecting the non-detection rate of 31% (11/35).

Using a bootstrapping method (see for details Stapper et al. 2022), a lognormal distribution is fit to the cumulative distributions to obtain a probability distribution, see the middle panel in Fig. 11. The fitting results are presented in Table 4. The best fit distributions are shown as the solid line. The fainter lines are included to demonstrate the range in possible fits. The dust distribution is slightly shifted toward lower masses when adding the five NOEMA targets, but the confidence intervals still overlap (in Stapper et al. 2022 μ = −2.18 ± 0.05 and σ = 0.53 ± 0.07). We find that the mean of the 100× dust mass distribution lies ~0.7 dex lower than the mean gas mass distribution, indicating that we do find more gas present than what would be assumed given 100× dust mass in most Herbig disks. This was also the observed trend in Figs. 7 and 10. Removing the gravitational unstable disks only moves the gas masses slightly down and makes the distribution slightly wider, see Table 4. The distribution of the minimum gas masses from our models lies ~0.2 dex below the mean of the dust mass distribution. Hence, the gas mass is likely higher than this value, consistent with a gas-to-dust ratio of close to or above 100.

To test this, we can sample the distributions obtained by fitting to the cumulative distributions and obtain a gas-to-dust ratio distribution by dividing the resulting gas mass values by the dust masses. After taking 10⁷ samples of each distribution, the gas-to-dust ratio distributions in the right panel of Fig. 11 are obtained. The canonical value of 100 for the gas-to-dust ratio falls within one standard deviation from the mean values of the different gas distributions obtained from our models. The mean of the gas-to-dust ratio made with the minimum gas mass values from the models differs an order of magnitude with the distribution using the mean gas masses. The gas-to-dust ratio distribution obtained from the gas masses based on the C¹⁸O flux is even lower at almost two orders of magnitude lower than those found by our models. This emphasizes the necessity of models to determine the disk mass, otherwise fundamental properties such as the gas-to-dust ratio can be severely underestimated.

Lastly, the higher than 100 gas-to-dust ratio may be originating from either optically thick dust or an increased size of dust grains in these disks. As more mass has grown into larger sized grains, the total mass visible at millimeter wavelengths decreases. The findings by Liu et al. (2022) and Kaeufer et al. (2023) show that an order of magnitude in mass can be hidden in the most massive disks. Taking multiple continuum wavelength dust observations into account our gas mass estimates are indeed close to 100 times the dust mass (Sierra et al. 2021). Thus, this order of magnitude difference between 100× the dust mass and the total disk mass is expected.

Fig. 11 Cumulative distributions of the gas masses from Figs. 7 and 10 obtained from our models (dark blue) and obtained from the C18O flux (green), and dust masses ×100 (orange) from Stapper et al. (2022). The corresponding probability distributions are shown in the middle panel, obtained by fitting a lognormal distribution to the cumulative distributions. The faint lines indicate the possible range in fits. Sampling from the gas and dust distributions a distribution of possible gas-to-dust mass ratios can be made, as shown in the right most panel.

5 Discussion

5.1 Comparison between different CO mass tracers

5.1.1 Rare CO isotopologues

Because of the disks in Fig. 5 being optically thick in both ¹³CO and C¹⁸O, a different way to measure the different masses of the disks is necessary. To be able to properly do this, one needs to have a handle on the size of the disk and. have an optically thin tracer. Here the rare(r) isotopologues C¹⁷O, ¹³C¹⁷O and ¹³C¹⁸O come into play.

Figure 12 shows the C¹⁷O,¹³C¹⁸O, and ¹³C¹⁷O J = 3–2 luminosities plotted against ¹²CO J = 2–1 of the DALI models, as the observations we compare these to have these transitions available. The optically thick tracer on the horizontal axis gives an indication of the size of the disk, while the optically thinner tracer on the vertical axis gives an indication of the mass of the disk. The same trends (e.g., models of the same mass curve downward for an increase in radius) can be seen as found in Sect. 4.4. For ¹²CO the smaller disks are now clearly separated into different areas of this parameter space. A clear stratification can be seen in the luminosities of the optically thinner isotopologues, where each mass has its own range in possible luminosities. However, even for these rare isotopologues, the most massive disks (>10^−1.5 M_⊙) have very similar luminosities. This may be due to the dust opacity playing an important role in reducing the luminosity of the rare isotopologues, which mostly emit from a layer close to the midplane.

We can use observations of rare isotopologues to see if these can assist us in obtaining a better measure of the disk mass. We use the following lines: C¹⁷O J = 3–2 for HD 100546 (Booth, A., priv. comm.), ¹³C¹⁸O J = 2–1 (Zhang et al. 2020a) and ¹³C¹⁷O J = 3–2 (Booth et al. 2019) for HD 163296, ¹³C¹⁸O J = 3–2 (Loomis et al. 2018) for HD 31648, and ¹³C¹⁸O J = 3–2 (Temmink et al. 2023) for HD 142527. We select the closest 100 models in luminosity around the observations as shown in Fig. 12, which reproduce the observed luminosities of the rare isotopologues to within 15%, except for the ¹³C¹⁷O J = 3–2 line of HD 163296 for which our models predict a factor of three lower luminosity. With these models we obtain ranges of possible disk masses very similar to the lower limits found with ¹³CO and C¹⁸O. For the highest mass models, the luminosities of the very rare isotopologues are still very similar for different mass disks.

However, for compacter lower mass disks (≲10⁻²), for which no rare isotopologues have been observed yet, the luminosity would result in a well constrained gas mass. Figure 13 shows the luminosity of the ¹³C¹⁷O J = 2–1 transition for different mass disks, selected again for different sizes in ¹²CO. This figure shows that the gas mass of a disk, if resolved in ¹²CO, can be constrained to within an order of magnitude given its ¹³C¹⁷O flux. An integration time of 1 h with ALMA typically gives a sensitivity of 1 mJy km s⁻¹ which at 150 pc gives a luminosity of 8.5 × 10² Jy km s⁻¹ pc². In Fig. 13 a 3σ detection is indicated with the dashed gray line. This shows that for the vast majority of disks, typically larger than ~300 au (Fig. 2) a detection is expected if the disk is more massive than 10⁻³ M_⊙ within an hour of observing time.

Fig. 12 Similar to Fig. 5, but now with C17O, 13C18O and 13C17O on the vertical axes and 12CO on the horizontal axis. The black dots mark HD 100546 observations of C17O J = 3–2 (Booth, A., priv. comm.), HD 163296 observations of 13C17O J = 3–2 from Booth et al. (2019), HD 31648 observations of 13C18O J = 3–2 from Loomis et al. (2020), and HD 142527 observations of 13C18O J = 3–2 by Temmink et al. (2023). Similar to Fig. 5 probability density curves are shown along the vertical and horizontal axes to indicate where most of the models of a particular mass reside.

5.1.2 Peeling the disk ‘onion’

As the previous section showed, an uncertainty of an order of magnitude is still present when determining the mass of the most massive disks. For some of those massive disks ¹³C¹⁷O is only marginally optically thin. Moreover, extrapolating the mass of the disk from the ¹³C¹⁷O emission which is mostly detected in the inner parts of the disk may not be accurate for the disk as a whole. Hence, this section explores how the disk gas surface density and mass can be constructed by considering, from the inner disk to the outer disk, a series of CO isotopologues with increasing abundance, from ¹³C¹⁷O or ¹³C¹⁸O to ¹²CO. Until now we neglected processes such as radial drift which can enhance the CO/H ratio inside the CO snowline (Zhang et al. 2021), and depleting it outside. This technique is able to take this depletion into account.

The left panel of Fig. 14 shows a disk model cut into sections of different CO isotopologues. The outer parts of the disk consist of the most abundant isotopologue, ¹²CO. The closer to the star, the rarer the isotopologue to ensure that the tracer stays as optically thin as possible. The regions are determined by the R_90% radii of the isotopologues. The right panels show the abundance of the six CO isotopologues looked at in this work, together with their τ = 1 lines and areas within which 50% and 95% of the emission is coming from.

Using Eq. (4) we compute the mass of each shell and add those together to obtain a measure of the mass of the disk. We obtain a total mass of 0.05 M_⊙ when combining the four regions shown in Fig. 14, which is a quarter of the mass of the model. This recovered mass comprises for 97% by mass out of ¹³C¹⁸O, with the last 3% made up of ¹³CO due to the larger emitting area. ¹²CO does not contribute any significant amount due to either high optical thickness or weak emission at large radii. Changing the rarest isotopologue into ¹³C¹⁷O does not change the estimated mass, for it is still dominated by the rarest isotopologue.

Doing the same analysis on the HD 163296 and HD 142527 disks, using an excitation temperature of 40 K (based on the bottom panel of Fig. 9) and the rare isotopologue observations by Booth et al. (2019) and Temmink et al. (2023), a mass of respectively 0.05 M_⊙ and 0.18 M_⊙ were found. This is a factor of three lower than what has been found in MAPS for HD 163296 (Zhang et al. 2021), and based on the spiral present one would expect a factor of 1.5 higher disk mass for HD 142527 (Yu et al. 2019). However, these masses are consistent with those found using ¹³CO and C¹⁸O, see Fig. 7.

In conclusion, rare isotopologues do trace the overall disk mass better than C¹⁸O. However, more modeling of these rare isotopologues needs to be done to properly use them.

Fig. 13 Same as Fig. 8, but for 13C17O J = 2–1. The gray dashed line indicates a typical ALMA 3a detection of 1 mJy km s−1 at a distance of 150 pc, which is achievable within an hour of integration time. The J = 3–2 transition can be found in Fig. C.2.

Fig. 14 DALI model with a mass of 0.2 M⊙. The disk is divided into rings based on the R90% of the isotopologue. The rarest isotopologue (13C18O) is in the center, while 12CO is in the outer ring. The six panels on the right show the abundance of the six ray-traced CO isotopologues. The vertical colored lines correspond to the colors of the circles on the left. The white lines indicate the τ = 1 surface, and the black lines indicate where 50% and 95% of the emission is coming from.

5.2 Comparison to other disk populations

5.2.1 T Tauri disks

Herbig disks are expected to be warmer than T Tauri disks, therefore causing less CO depletion due to freeze-out. Using the dust mass estimates from Ansdell et al. (2016) for Lupus and Pascucci et al. (2016) for Chameaeleon I together with the gas mass estimates by Miotello et al. (2017) and Long et al. (2017) (who use the models from Miotello et al. 2016), gas-to-dust ratios of T-Tauri disks are obtained. Figure 15 shows these gas-to-dust ratios plotted against the inferred dust masses together with the gas-to-dust ratios obtained for the Herbig disks. This comparison is also done in Miotello et al. (2023), but we increase the number of Herbigs by more than a factor of two.

As was mentioned before, the mean gas-to-dust ratio of the Herbig disks lies above the canonical ISM value of 100, but many of the disks are still consistent with a gas-to-dust ratio of 100. The T Tauri disks on the other hand lie at least an order of magnitude below the gas-to-dust ratio of the ISM. Moreover, the data shown in Fig. 15 only includes detected disks. The many non-detections among the T Tauri disks suggests a lack of CO gas in these disks. The fact that orders of magnitude differences in the gas-to-dust ratio are found, over multiple orders of magnitude in dust mass, confirms the expectation that the warmer Herbig disks lack the CO-conversion processes that occur in the cold T Tauri disks.

Fig. 15 Gas-to-dust ratio of the Herbig disks compared to T Tauri disks in Lupus and Chamaeleon I (Ansdell et al. 2016; Pascucci et al. 2016; Miotello et al. 2017; Long et al. 2017). The corresponding horizontal dashed lines are the logarithmic mean values. The canonical value of 100 is indicated as the horizontal black dotted line.

5.2.2 Group | versus group ||

Herbig disks can be separated into two different groups based on their spectral energy distribution (SED): group I have a rising SED in the far-infrared, while group II disks do not (Meeus et al. 2001). Stapper et al. (2022) showed a difference in disk dust mass between group I and group II disks, where the group II disks have a dust distribution very similar to T Tauri disks. Moreover, group I disks are generally large disks with large cavities (transition disks), while group II disks are more compact (Stapper et al. 2022; Garufi et al. 2017).

Similar to the distributions in Fig. 11, we can obtain cumulative distributions for the separate groups as well, see Fig. 16. Both the gas distributions in red (group I) and blue (group II), and the dust distributions in dark gray (group I) and light gray (group II) are shown. The dust mass distributions are relatively well constrained, while the gas distributions much less so. The gas distributions are shown in the left most panel in Fig. 16 for Lupus (orange; Miotello et al. 2017) and Chamaeleon I (green; Long et al. 2017).

Comparing the group I and group II disks, the group II disks are less massive than the group I disks, which is consistent with their dust masses. The probability distributions show an offset in the mean of the distribution of an order of magnitude for both dust and gas (see Table 5). Hence, regardless of the differences in dust and gas masses between the two groups, the relative gas-to-dust ratio remains the same. As the right most panel of Fig. 16 shows, the mean of the gas-to-dust ratio distributions are very similar (see Table 5). While Stapper et al. (2022) found that the group II disks have a very similar dust mass distribution as T Tauri disks, compared to group I disks they have an order of magnitude lower dust mass. We find here that the gas distribution of the group II disks is shifted toward lower gas masses compared to the group I disks by the same factor (see Table 5). This further supports that the gas-to-dust ratio is disk mass independent and rather a result from differences in temperature compared to T Tauri disks.

Some group I disks have been found to have little to no CO depletion (e.g., HD 100546 Booth et al. 2023a; HD 169142 Carney et al. 2018; Booth et al. 2023b), while group II disks do in the outer disk (Zhang et al. 2021). As discussed in Stapper et al. (2022), group I disks show in general large cavities in millimeter emission, while group II disks are more compact. These differences in CO depletion between the two groups could indicate a much higher impact due to radial drift where most CO stayed in the outer disk for group I disks and the CO drifted inward for group II disks. This is further supported by observations of the metallicities of the hosting Herbig stars (Kama et al. 2015; Guzmán-Díaz et al. 2023), for which low metallicities generally coincide with group I disks. Hence, group I disks might be the formation sites of giant exoplanets stopping radial drift, stopping the enrichment of the host star, trapping CO in the outer disk, and creating the quintessential large cavity structure often associated with these disks. The fact that we find higher gas masses of group I disks compared to group II disks only further supports this hypothesis.

Regarding the sizes of the group I and group II disks, Brittain et al. (2023) report a comparison of the ¹³CO and dust radii for 17 Herbig disks. They find that the group II disks are the smallest disks in both gas and dust and that in general the ratio between the dust and gas radii are consistent with the Lupus disks. Using ¹²CO radii we do not report a similar difference between the two groups in gas observations (see Fig. 4) as most of the smallest disks lack any detection of ¹²CO. Hence, this remains inconclusive. A uniform survey of Herbig disks would help in characterizing these differences.

Fig. 16 Cumulative distributions of the dust and gas masses of the Herbigs separated into group I and group II. The dust distributions are shown in the left and middle panel in dark gray for group I and light gray for group II. In addition, the gas distributions from Lupus (orange) and Chamaeleon I (green) is shown (Miotello et al. 2017; Long et al. 2017). The fitted probability distributions are shown in the middle panel. The resulting gas-to-dust ratios are presented in the right panel.

5.2.3 Debris disks

Debris disks are the final stage of planet-forming disks and are sustained by collisional processes producing secondary dust (Hughes et al. 2018). In contrast to later spectral type stars, debris disks around A-type stars (i.e., the evolutionary successors of Herbig disks) are more common to be detected in CO gas compared to later spectral type stars (e.g., Moór et al. 2017, 2020). Whether this is primordial or secondary gas is still heavily debated (see for an overview Hughes et al. 2018). What has been found is that the amount of CO mass in the (debris) disks around A-type stars rapidly decreases between Herbig disks and debris disks around A-type stars (e.g., Moór et al. 2020). In this section we compare the CO masses found by Moór et al. (2020) and Cataldi et al. (2023) to the disk masses we obtained for Herbig disks.

Using the previously adopted CO isotopologue and CO/H₂ ratios, we compute the total disk mass from the CO mass estimates by Moór et al. (2017), see Fig. 17. Moór et al. (2020) showed that most debris disks with an A spectral type are two orders of magnitude less massive than Herbig disks in CO. Here we show that this difference may be even more dramatic, with a difference between the logarithmic means of the two population of four orders of magnitude. In Fig. 17 we additionally show the LTE models from Cataldi et al. (2023) who modeled the CO and CI emission of 14 debris disks to obtain a CO mass in the disk. After converting these into total disk mass, large differences between the Herbig disks and the debris disks can be seen, as is the case for the masses from Moór et al. (2017). The fact that few disks are in between the Herbig disks and debris disks suggests that the disk needs to dissipate quickly over just a few million years.

Lastly we note that if the disk consists of second generation gas, the assumption of CO/H₂ = 10⁻⁴ might not be correct and should rather be closer to one, but the assumption does give an upper limit on the gas mass in the disk. In the case of CO/H₂ = 1 the debris disk masses shown in Fig. 17 are even lower.

Fig. 17 Comparison of the gas masses and ages of debris disks from Moór et al. (2017) and Cataldi et al. (2023) to our Herbig disk gas masses. The horizontal lines indicate the logarithmic means.

Fig. 18 Upper limits on the disk wind mass. The mass of the disk wind of HD 163296 is shown as the dotted horizontal line (Booth et al. 2021). Disks which have an accretion rate higher or lower than that of HD 163296 are shown as blue or red respectively (Guzmán-Díaz et al. 2021). HD 34282 has a lower limit on its accretion rate. Disks with no accretion rate measured are shown in gray.

5.3 Viscous or wind driven evolution

Finally, we discuss the viscous versus wind driven evolution of disks. Wind signatures have been found in embedded disks but for later stage disks, in particular Herbig disks, direct evidence for a disk wind has only been found in HD 163296 in ¹²CO and ¹³CO (Klaassen et al. 2013; Booth et al. 2021). As our work has compiled all Herbig disks with CO observations available, we present in this section the results of a search for wind signatures in other Herbig disks.

Following previous works (Klaassen et al. 2013; Booth et al. 2021), the search for wind signatures is done with visibility spectra. To ensure that the visibility spectrum is dominated by the large-scale structure of the wind, only the short baselines are selected. The 20% shortest baselines are used for the analysis, with two exceptions: AK Sco and HD 142666. Due to their relatively high resolution, a lower percentage is chosen such that all baselines smaller than respectively 80 and 100 meter are selected, corresponding to spatial scales of ~2.0″. To obtain the visibility spectra, the time and baselines are averaged to obtain a better S/N. The velocity is sampled from around -50 km s⁻¹ to 50 km s⁻¹ relative to the system velocities, so that all the possible disk wind channels can be covered (Pascucci et al. 2023).

No wind signatures have been found, either due to low-quality data or no presence of a disk wind. To assess this we obtain a measure of the noise from the visibility spectra of each disk scaled by the square-root of the velocity resolution relative to that of the HD 163296 observation. Using the estimate of the total mass of the disk wind from Booth et al. (2021) of 10⁻³ M_⊙ and the peak flux in their spectrum of 0.14 Jy in ¹³CO, we can scale the upper limits of 5× the noise with the distance to the object relative to HD 163296 (see Table 1), assuming that CO/H₂ is not different in each disk wind. This results in an upper limit on the disk wind in Earth masses for each disk. Figure 18 presents the results.

The upper limits on the majority of the disks fall below the mass of the disk wind of HD 163296. Of these disks, a third have accretion rates lower than HD 163296 (see the red scatter points, from Guzmán-Díaz et al. 2021). This could result in the non-detections as a stronger disk wind induces a higher accretion rate. For these disks we can rule out the sensitivity of the data being the main reason for not detecting a disk wind similar to that of HD 163296. On the other hand, nine disks have upper limits above the mass of HD 163296, which does indicate a lack of sensitivity and a disk wind similar to that of HD 163296 could possibly still be present.

Comparing these results to the gas masses found in Sect. 4.4, it is clear that only a few disks are likely more massive than HD 163296, most of which have a better upper limit on the disk wind mass than lower mass disks. For example HD 290764 and HD 245185 may be good candidates for deeper follow-up surveys of disk winds in Herbig disks due to their high gas mass and accretion rates, but relatively high upper limits. Another option would be HD 34282, which has a high disk mass but a lower accretion rate compared to HD 163296. Notably, Pegues et al. (2023) identify a tentative extended structure in ¹²CO in HD 34282, while not seen in our data this further substantiates a follow-up. For some other disks such as HD 142666 longer and lower resolution observations to obtain better short baseline coverage and higher sensitivities would also be useful in order to confirm a lack of disk wind signatures.

Using the gas radii in combination with the ages from Guzmán-Díaz et al. (2021), we determine if the radii of Herbig disks evolve over their lifetime, and if so, whether a specific evolutionary scenario is favorable. Viscously evolving disks are expected to increase in size over time, while wind driven evolution results in a reduction of the size of the disk (e.g., Manara et al. 2023). Trapman et al. (2022) modeled the evolution of the CO radius in the wind driven case, showing that the size of the disk is indeed expected to decrease in this specific tracer. Though effects of external photoevaporation are important to keep in mind (Trapman et al. 2023). Figure 19 shows the gas radius plotted against the age of the system. The Herbig disks younger than 12 Myr do seem to correspond well with the wind evolution shown by Trapman et al. (2022), the largest and highest mass disks decrease in size by a factor of ~4 between 2 and 10 Myr. For the three disks on the right hand side of Fig. 19, two (HD 34282 and HD 169142) only have upper limits on their age. However, other works do put these at younger ages of ~10 Myr (e.g., Vioque et al. 2018). Samples of older ages are missing, which is needed to constrain the specific evolutionary scenario. More sensitive observations of a larger sample are needed to make progress.

Fig. 19 Gas radii (R90%) of the Herbig disks plotted against the age of the system (from Guzmán-Díaz et al. 2021). The downward facing triangles indicate an upper limit on the gas radius, an upper limit on the age is indicated as a left facing triangle.

6 Conclusion

In this work we analyze the ¹²CO ,¹³CO and C¹⁸O J = 2–1 or J = 3–2 emission in 35 Herbig disks, 30 with ALMA archival data and five new datasets observed with NOEMA. We compare the integrated line luminosities and R_90% radii of ¹²CO and ¹³CO with a large grid of DALI models (Bruderer et al. 2012; Bruderer 2013) to obtain a measure of the gas mass. We can conclude the following:

• 1.

  We detect ¹²CO emission in 20 out of the 27 disks which have the line covered, for ¹³CO in 22 disks out of 33 disks, and for C^l8O in 21 disks out of 33 disks. In total, 15 disks are resolved in ¹²CO, 16 in ¹³CO, and 15 in C¹⁸O. For all resolved disks, the ¹²CO emission extends to larger radii compared to ¹³CO, which in turn is larger than the C¹⁸O disk.

• 2.

  The main model parameters affecting the luminosities of the ¹³CO and C¹⁸O lines in the models are the mass and size of the disk. In addition, the power-law index of the surface density affects the line luminosities due to changing the distribution of the mass in the disk. For the most massive disks we find that the effect of the stellar luminosity, the vertical distribution of the disk mass, and the inclination of the disk have negligible impact on the line luminosities. For low mass disks however, increasing the vertical distribution and stellar luminosity decreases the luminosity of the isotopologue lines by one to two orders of magnitude.

• 3.

  The two deciding processes influencing the line luminosity of the disk is are the line optical depth as seen from the observer and the self-shielding capacity of the CO isotopo-logues. Hence, one can make disks with similar luminosities with very different masses if the disk is optically thick in both ¹³CO and C¹⁸O. When enlarging the disk, first C¹⁸O becomes optically thin and starts to photodissociate, reducing the C¹⁸O line luminosity. At even larger sizes, ¹³CO becomes optically thin as well and subsequently its luminosity reduces.

• 4.

  We find that most of the detected Herbig disks are optically thick in both ¹³CO and C¹⁸O. For almost all disks we can only find a lower disk mass which can reproduce the observed line luminosities. The R_90% size of the disk is an essential observable to constrain the gas mass, as only the most massive disks can be large in ¹²CO.

• 5.

  Comparing the gas masses to those obtained from the number of C¹⁸O molecules with a simple CO/H₂ conversion shows that the disk mass is generally underestimated by at least an order of magnitude. This shows that to obtain the mass of a disk, modeling the disk is vital.

• 6.

  Combining the gas masses of the disks with the dust masses from Stapper et al. (2022), we find that Herbig disks are consistent with the canonical gas-to-dust ratio of 100. In general the ratio is even higher, which may be caused by a combination of dust optical depth and grain growth.

• 7.

  Comparing the gas radii with the dust radii from Stapper et al. (2022) we find a ratio of 2.7, higher compared to the disks in Lupus (a factor of 2.0, Ansdell et al. 2018). Still, the majority of the disks fall well below the factor of four, indicating that this difference may only be due to line optical depth effects rather than radial drift.

• 8.

  To distinguish different disk masses for optically thick disks, a combination of ¹²CO tracing the size of the disk and ¹³C¹⁷O tracing the mass of the disk would make this possible for a large range in masses of disks. However, for the most massive disks dust opacity may inhibit tracing the disk mass with such a rare isotopologue.

• 9.

  Comparing the Herbig gas-to-dust ratios with those in T Tauri disks (Ansdell et al. 2016; Pascucci et al. 2016; Miotello et al. 2017; Long et al. 2017), we find that Herbig disks have a gas-to-dust ratio of almost two orders of magnitude higher over a range of multiple orders of magnitude in dust mass. Hence, this disparity could be caused by a fundamental difference in chemistry due to Herbig disks being much warmer, as proposed in previous works.

• 10.

  The gas and dust masses of the two different Meeus et al. (2001) groups are found both differ by an order of magnitude. This results in the same gas-to-dust ratios for both groups, even though group II disks have similar dust masses as T Tauri disks (Stapper et al. 2022). This further supports the idea that the gas-to-dust ratio is disk mass independent and that the lack of CO emission in T Tauri disks is not due to their lower mass disks, but rather a temperature difference compared to Herbig disks.

• 11.

  The total masses of debris disks are found to be at least four orders of magnitude lower than those of Herbig disks, indicating a rapid dissipation of the material in the disk within a few Myr. Full chemical modeling of debris disks is necessary to explore this further.

• 12.

  A search for disk wind signatures such as those found for HD 163296 (Booth et al. 2021) in the ¹²CO data has resulted in no additional detections. Most disks have sufficiently sensitive data in which a disk wind analogous to that in HD 163296 would have been detected. This lack of disk wind may be in some of the disks related to a difference in accretion rate.

• 13.

  Comparing the gas radii to the age of the system we find that the data seem to support a disk wind driven evolution, but data of older age systems are lacking.

In conclusion, the warmer Herbig disks have significantly larger gas-to-dust mass ratios compared to the colder T Tauri stars, close to, or exceeding, the canonical value of 100.

Acknowledgements

The research of LMS is supported by the Netherlands Research School for Astronomy (NOVA). This paper makes use of the following ALMA data: 2012.1.00158.S, 2012.1.00303.S, 2012.1.00698.S, 2012.1. 00870.S, 2013.1.00498.S, 2015.1.00192.S, 2015.1.00222.S, 2015.1.00986.S, 2015. 1.01058.S, 2015.1.01353.S, 2015.1.01600.S, 2016.1.00110.S, 2016.1.00204.S, 2016.1.00344.S, 2016.1.00484.L, 2016.1.00724.S, 2017.1.00466.S, 2017.1. 00940.S, 2017.1.01404.S, 2017.1.01419.S, 2017.1.01607.S, 2018.1.00814.S, 2018. 1.01055.L, 2018.1.01222.S, 2019.1.00218.S, 2019.1.00579.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. This work is based on observations carried out under project number S21AS with the IRAM NOEMA Interferometer. IRAM is supported by INSU/CNRS (France), MPG (Germany) and IGN (Spain). This work makes use of the following software: The Common Astronomy Software Applications (CASA) package (McMullin et al. 2007), Dust And LInes (DALI; Bruderer et al. 2012; Bruderer 2013), Python version 3.9, astropy (Astropy Collaboration 2013, 2018), cmasher (van der Velden 2020), Lifelines (Davidson-Pilon et al. 2021), matplotlib (Hunter 2007), numpy (Harris et al. 2020), pandas (Pandas development team, T. 2020), scipy (Virtanen et al. 2020) and seaborn (Waskom 2021). Lastly, we thank the referee for their insightful comments which have improved this paper.

Appendix A Datasets used

In Table A.1 the project codes of the used datasets are listed together with their spatial and velocity resolution and the line-free rms noise.

Appendix B NOEMA continuum data

Fig. B.1 Continuum images of the five northern Herbig disks observed with NOEMA.

Figure B.1 presents the continuum images of the five northern Herbig disks observed with NOEMA during the summer semester of 2021 (PI: Cridland, Booth; Project code: S21AS). XY Per was observed on the 17th of November 2021, for 3.4 hours. Both the bandpass and gain calibrator used was 3c84, while the flux calibrator was MWC349. BH Cep, BO Cep, HD 200775, and SV Cep were observed on the 18th of October, for 1.6 hours. The bandpass calibrator was 3C454.3, the gain calibrator was 2010+723, and the flux calibrator was MWC349. The observations were done in the C configuration, with nine antennas.

The data imaging was done using the Common Astronomy Software Applications (CASA) version 5.8.0 (McMullin et al. 2007). Forthe data one round ofphase-only self-calibration was done. After applying the resulting calibration table to all spectral windows, the data were cleaned using the hogbom algorithm and imaged using the multifrequency synthesis spectral definition mode and a Briggs robust weighting of 0.5. The resulting beam and rms of each continuum observation can be found in Table B.1. Four out of five Herbig disks have been detected. XY Per is a binary with the A component being a Herbig star.

The integrated fluxes and their corresponding dust masses obtained by following the same procedure as Stapper et al. (2022) can be found in Table B.1. Because the disks are unresolved, the median upper limits on their size are 245 au and 373 au for the 68% and 90% radii respectively. Comparing the dust masses to the distribution found by Stapper et al. (2022), we find that only BO Cep is more massive than the mean dust mass, which is likely related to it being the only disk in our sample observed with NOEMA for which the ¹³CO and C¹⁸O isotopo-logues are detected. XY Per is slightly less massive than the mean dust mass from Stapper et al. (2022), but still well above the dust mass for which we still detect the CO isotopologues, which we do not detect in this disk possibly caused by the binary nature of this object. SV Cep and HD 200775 both have relatively low disk masses. Interestingly, in the mid-infrared HD 200775 was found to have diffuse emission going out to ~700 au (2″, corrected for the most recent distance estimate) in both north and south direction, and a large tail extending to the north-east (~10″, Okamoto et al. 2009). We do not resolve the continuum emission with a beam of 1.1″ × 0.8″, and thus also do not see any of these structures. Lastly, for BH Cep we determined an upper limit on the dust mass of 1.5 M_⊕ which is higher than some of the detection made by ALMA (Stapper et al. 2022).

Appendix C Figures for the J = 3−2 transition

Figure C.1 presents the data and models for the J = 3−2 transition, similar to Fig. 5. The two disks plotted are HD 135344B and HD 290764. Figures C.2 and C.3 are the same as Figs. 8 and 13 respectively but for the J = 3−2 transition.

Fig. C.1 Same as Figure 5, but for the J = 3−2 transition.

Fig. C.2 Same as Figure 8, but for the J = 3–2 transition.

Fig. C.3 Same as Figure 13, but for the J = 3−2 transition.

Appendix D Parameter overview plots

This appendix presents plots showing in which direction the parameters from Table 3 change the ¹³CO and C¹⁸O luminosities, similar to Fig. 6. Figures D.1 to D.4 are the same as Fig. 6, but for the other R_c values used. Similarly, Figures D.5 and D.6 are the same as Fig. 6, but for the other γ values. Lastly, Figures D.7 to D.13 show the same figures but for the J = 3−2 transition.

Fig. D.1 Same as Fig. 6, but with Rc = 5 au

Fig. D.2 Same as Fig. 6, but with Rc = 10 au

Fig. D.3 Same as Fig. 6, but with Rc = 30 au

Fig. D.4 Same as Fig. 6, but with Rc = 60 au

Fig. D.5 Same as Fig. 6, but with γ = 0.4

Fig. D.6 Same as Fig. 6, but with γ = 1.5

Fig. D.7 Same as Fig. 6, but with the J = 3 − 2 transition.

Fig. D.8 Same as Fig. D.7, but with Rc = 5 au

Fig. D.9 Same as Fig. D.7, but with Rc = 10 au

Fig. D.10 Same as Fig. D.7, but with Rc = 30 au

Fig. D.11 Same as Fig. D.7, but with Rc = 60 au

Fig. D.12 Same as Fig. D.7, but with γ = 0.4

Fig. D.13 Same as Fig. D.7, but with γ = 1.5

Appendix E Disk masses

Table E.1 presents the mean disk masses and range in possible disk masses based on our models for all Herbig disks presented in this work.

References

All Tables

All Figures

Fig. 1 Keplerian masked velocity integrated-intensity maps of the 23 Herbig disks with at least one detection of 12CO, 13CO, or C18O. For each disk, 12CO, 13CO, and C18O are given in the first, second, and third panel from left to right. If nothing is shown, the molecule is not covered. The bar in every first panel from the left indicates a size of 100 au with the corresponding angular size below. The size of the beam is indicated in the bottom left of each panel. A sinh stretch is used to make the fainter parts of the maps better visible. The minimum value is set to zero, and the maximum value in mJy beam−1 km s−1 is indicted on each panel in the top left corner. The disks with no detections in all three isotopologues are not shown here.

In the text

	Fig. 2 Summary of the luminosities and 90% radii for the 12CO, 13CO, and C18O observations as presented in Table 2.
In the text

Fig. 3 Azimuthally averaged radial profiles of the three CO isotopologues in all 23 Herbig disks in which at least one isotopologue is detected. The 1σ uncertainty interval is indicated by the shaded region in the same color as the profile. The vertical solid lines indicate the derived 90% radii for each line. An upper limit on the radius is indicated with a left facing arrow. The beam size is illustrated by the horizontal lines in the top right of each panel.

In the text

Fig. 4 Gas radii versus dust radii of the Herbig disks with detected 12CO emission, both R90%. An upper limit on the gas radius or dust radius, or both, is marked with a triangle pointing in the corresponding direction(s). The colors indicate the Meeus et al. (2001) group of the disk: red is group I, blue group II. Lines for different Rdust to Rgas ratios are shown as well. The black line shows for the resolved disks the relation of the R90% radii. Using R68% does not alter the relation significantly. The yellow region indicates the region where the difference between the dust and gas radii cannot solely be explained by optical depth effects, and radial drift is necessary (Trapman et al. 2019). The red line corresponds to the relation found for the disks in Lupus (Ansdell et al. 2018).

In the text

Fig. 5 C18O luminosity versus 13CO luminosity for the DALI models (colored circles) and the observations (black circles). The colors indicate the gas mass of the disk model. The size of the markers indicate the 90% radius of the 13CO emission. Probability density curves of the models for each gas mass are shown along the x-axis in the top panel and along the y-axis in the right panel, i.e., most models of a particular gas mass reside around the peak of each curve.

In the text

Fig. 6 Overview of the different parameters explored in the DALI models and their effect on the luminosity of 13CO and C18O. The colors indicate the mass of the model. The arrows indicate in which direction the parameter increases in value (see Table 3). The black outlined circles are the same models for each mass in all panels. The black outlined models have the following parameters: γ = 0.8, ψ = 0.2, hc = 0.2 rad, Rc = 200 au, i = 10° and L★ = 10 L⊙. In Appendix D the same plot is shown, but with γ = 0.4 and γ = 1.5, and with Rc = 5 au, Rc = 10 au, Rc = 30 au, and Rc = 60 au, and for the J = 3−2 transition.

In the text

Fig. 7 Gas masses of the Herbig disks detected in the CO isotopologues as selected from the models presented in Fig. 5. The disks are ordered from left to right by increasing 13CO luminosity. The gray lines indicate the range in possible disk mass based on the disk parameters listed in Table 1, while the black lines indicate the standard deviation. The mean mass of these models is given as the dark blue diamond. After removing models which are potentially gravitationally unstable, the mean mass changes into the light blue diamond. The orange diamonds are taken from the dust masses of Stapper et al. (2022) multiplied by the canonical gas-to-dust ratio of 100. The HD upper limits from Kama et al. (2020) are shown as green downward pointing triangles. The C18O lower limits are directly obtained from the integrated flux. The middle panel shows the resulting gas-to-dust ratio based on the mass range and mean values, and the top panel shows the factor by which the mass could be underestimated when using C18O assuming optically thin emission.

In the text

	Fig. 8 Range of C18O luminosities for different mass disks, selected for 12CO R90% radii from the DALI models presented in Fig. 5. The shaded area indicates the minimum and maximum values, and the solid line is the mean value.
In the text

Fig. 9 Gas mass of the models versus the retrieved gas mass from the C18O flux using Eq. (4), and the mass-weighted average of the temperature in the C18O emitting region. For the blue shaded region a temperature of 40 K is used, the gray shaded region uses the temperature of each model shown in the bottom panel. The shaded areas are the minimum and maximum values, with the mean value indicated by the solid line. The dashed line indicates the one-to-one correspondence between the two masses. The typical range in temperatures is indicated by the gray shaded region in the bottom panel.

In the text

	Fig. 10 Same as Fig. 7, but now for the disks with only upper limits on the 13CO and C18O fluxes.
In the text

Fig. 11 Cumulative distributions of the gas masses from Figs. 7 and 10 obtained from our models (dark blue) and obtained from the C18O flux (green), and dust masses ×100 (orange) from Stapper et al. (2022). The corresponding probability distributions are shown in the middle panel, obtained by fitting a lognormal distribution to the cumulative distributions. The faint lines indicate the possible range in fits. Sampling from the gas and dust distributions a distribution of possible gas-to-dust mass ratios can be made, as shown in the right most panel.

In the text

Fig. 12 Similar to Fig. 5, but now with C17O, 13C18O and 13C17O on the vertical axes and 12CO on the horizontal axis. The black dots mark HD 100546 observations of C17O J = 3–2 (Booth, A., priv. comm.), HD 163296 observations of 13C17O J = 3–2 from Booth et al. (2019), HD 31648 observations of 13C18O J = 3–2 from Loomis et al. (2020), and HD 142527 observations of 13C18O J = 3–2 by Temmink et al. (2023). Similar to Fig. 5 probability density curves are shown along the vertical and horizontal axes to indicate where most of the models of a particular mass reside.

In the text

	Fig. 13 Same as Fig. 8, but for 13C17O J = 2–1. The gray dashed line indicates a typical ALMA 3a detection of 1 mJy km s−1 at a distance of 150 pc, which is achievable within an hour of integration time. The J = 3–2 transition can be found in Fig. C.2.
In the text

Fig. 14 DALI model with a mass of 0.2 M⊙. The disk is divided into rings based on the R90% of the isotopologue. The rarest isotopologue (13C18O) is in the center, while 12CO is in the outer ring. The six panels on the right show the abundance of the six ray-traced CO isotopologues. The vertical colored lines correspond to the colors of the circles on the left. The white lines indicate the τ = 1 surface, and the black lines indicate where 50% and 95% of the emission is coming from.

In the text

	Fig. 15 Gas-to-dust ratio of the Herbig disks compared to T Tauri disks in Lupus and Chamaeleon I (Ansdell et al. 2016; Pascucci et al. 2016; Miotello et al. 2017; Long et al. 2017). The corresponding horizontal dashed lines are the logarithmic mean values. The canonical value of 100 is indicated as the horizontal black dotted line.
In the text

Fig. 16 Cumulative distributions of the dust and gas masses of the Herbigs separated into group I and group II. The dust distributions are shown in the left and middle panel in dark gray for group I and light gray for group II. In addition, the gas distributions from Lupus (orange) and Chamaeleon I (green) is shown (Miotello et al. 2017; Long et al. 2017). The fitted probability distributions are shown in the middle panel. The resulting gas-to-dust ratios are presented in the right panel.

In the text

	Fig. 17 Comparison of the gas masses and ages of debris disks from Moór et al. (2017) and Cataldi et al. (2023) to our Herbig disk gas masses. The horizontal lines indicate the logarithmic means.
In the text

	Fig. 18 Upper limits on the disk wind mass. The mass of the disk wind of HD 163296 is shown as the dotted horizontal line (Booth et al. 2021). Disks which have an accretion rate higher or lower than that of HD 163296 are shown as blue or red respectively (Guzmán-Díaz et al. 2021). HD 34282 has a lower limit on its accretion rate. Disks with no accretion rate measured are shown in gray.
In the text

	Fig. 19 Gas radii (R90%) of the Herbig disks plotted against the age of the system (from Guzmán-Díaz et al. 2021). The downward facing triangles indicate an upper limit on the gas radius, an upper limit on the age is indicated as a left facing triangle.
In the text

	Fig. B.1 Continuum images of the five northern Herbig disks observed with NOEMA.
In the text

	Fig. C.1 Same as Figure 5, but for the J = 3−2 transition.
In the text

	Fig. C.2 Same as Figure 8, but for the J = 3–2 transition.
In the text

	Fig. C.3 Same as Figure 13, but for the J = 3−2 transition.
In the text

	Fig. D.1 Same as Fig. 6, but with Rc = 5 au
In the text

	Fig. D.2 Same as Fig. 6, but with Rc = 10 au
In the text

	Fig. D.3 Same as Fig. 6, but with Rc = 30 au
In the text

	Fig. D.4 Same as Fig. 6, but with Rc = 60 au
In the text

	Fig. D.5 Same as Fig. 6, but with γ = 0.4
In the text

	Fig. D.6 Same as Fig. 6, but with γ = 1.5
In the text

	Fig. D.7 Same as Fig. 6, but with the J = 3 − 2 transition.
In the text

	Fig. D.8 Same as Fig. D.7, but with Rc = 5 au
In the text

	Fig. D.9 Same as Fig. D.7, but with Rc = 10 au
In the text

	Fig. D.10 Same as Fig. D.7, but with Rc = 30 au
In the text

	Fig. D.11 Same as Fig. D.7, but with Rc = 60 au
In the text

	Fig. D.12 Same as Fig. D.7, but with γ = 0.4
In the text

	Fig. D.13 Same as Fig. D.7, but with γ = 1.5
In the text