© The Authors 2026

1 Introduction

Previous observations have shown that the circumstellar envelopes (CSEs) around asymptotic giant branch (AGB) stars remain largely spherical on large scales, despite showing asymmetries, and typically expand isotropically at low velocities (10-20 km s⁻¹, see Sahai et al. 2007; Castro-Carrizo et al. 2010; Decin et al. 2020; Danilovich et al. 2025). Their sucessors, the pre-planetary and planetary nebulae (pPNe and PNe), present different morphologies: pPNe tend to present strongly aspherical and (often) axisymmetrical shapes resulting from the interaction of axial outflows (30-100 km s⁻¹) with the CSE formed in the AGB stage; PNe often exhibit shapes that are nearly spherical, axisymmetric, or irregular (see e.g. Ueta et al. 2000; Sahai et al. 2007; Stanghellini et al. 2016). This transition lasts only a few thousand years (Balick & Frank 2002).

An oft-acknowledged explanation for this evolution is that a companion, either a main-sequence star or a degenerate compact object, accumulates material from a rotating disk. Subsequently, extremely fast jets are propelled, akin to the process observed in protostars (Soker 2002; Frank & Blackman 2004; Blackman & Lucchini 2014). Our understanding of the impact of binary systems on post-AGB stars remains limited. Nevertheless, the axisymmetric outflows observed in pPNe and PNe are statistically consistent with the prevalence of binaries in the low- to intermediate-mass star population (De Marco et al. 2013; De Marco & Izzard 2017).

Particularly remarkable, certain post-AGB stars, specifically those in binary systems (binary post-AGB stars), consistently exhibit signs of having disks that orbit their central stars. All of these sources show a notable excess in near-infrared (NIR) and narrow CO line profiles indicative of rotating disks (see Van Winckel 2003; Bujarrabal et al. 2013, and references therein). Analysis of their spectral energy distributions (SEDs) indicates the existence of warm dust in proximity to the stellar system, and interferometric infrared data confirm its disk-like structure (Hillen et al. 2017; Corporaal et al. 2021; Kluska et al. 2022). These disks are considered stable structures, as is indicated by observations showing significant dust processing, including a high degree of crystallization and the presence of large grains (see e.g. Gielen et al. 2011a,b). Additionally, single-dish COline observations affirm the widespread occurrence of rotation (Bujarrabal et al. 2013), as the characteristics of the line profiles resemble those found in young stars encircled by a rotating disk (Guilloteau et al. 2013). Both post-AGB disks and disks around young stars may look similar, but they present significant differences in practically all physical-chemical parameters, such as origin, lifetime, chemical abundances, temperature, density, mass, or even size (Gallardo Cava 2023).

Although the group is very homogeneous, a close inspection of their luminosities has revealed the presence of dusty postRGB primaries in the sample, as they lie below the canonical AGB lower luminosity limit of about 2500 L_⊙ (Kamath et al. 2016). However, it must be noted that these luminosities, derived through SED fitting (Kluska et al. 2022), sometimes present substantial uncertainties due to errors in the distance and reddening corrections (see Van Winckel 2025), so some sources could in fact lie above or below this threshold. In view of this and the structural similarities, we retain the term post-AGB binaries for all dusty evolved binaries discussed in this paper.

The detailed study of these circumbinary post-AGB disks requires interferometric observations (at both radio and opti-cal/IR wavelengths, see e.g. Bujarrabal et al. 2018; Kluska et al. 2019). According to millimeter (mm)-wave interferometric maps, we know that there is an additional extended structure surrounding the rotating disk, which presents low expansion velocities. These outflows may be composed of gas and dust extracted from the disk (see discussion in Bujarrabal et al. 2016). In this work, we refer to the nebula as all the circumstellar material, including both the well-defined rotating structure that forms the disk and the expanding component. Recent studies classify these post-AGB binaries into two subclasses: disk-dominated and outflow-dominated, depending on which structure holds most of the total nebular material. In the disk-dominated subclass, the disk holds the majority of the nebular mass, ranging from 85% to 95% of the total nebular material. In the case of outflow-dominated sources, it is the extended and expanding component that holds the majority of the detected nebular material: 65% to 75% (see Gallardo Cava et al. 2021, 2023).

A remarkable increase in the number of nebulae around post-AGB binaries studied in CO through mm-CO single-dish observations is reported in this work. This is important, as these objects have not been widely studied using single-dish observations, except for the previous CO survey (Bujarrabal et al. 2013), or the more recent survey searching for other species by Gallardo Cava et al. (2022b), which focused on molecular content.

Motivated by these objectives, this paper is laid out as follows. We describe the observational techniques used in Sect. 2. We introduce the observed sample in Sect. 3. In this section, we also provide the main parameters derived from the final CO spectra. In Sect. 4, we calculate the nebular masses of these diskcontaining post-AGB nebulae. Next, in Sect. 5, we compare the nebular mass values obtained in this work with those from earlier studies. In Sect. 6, we discuss the presence of molecules other than CO. Lastly, we provide a concise overview of our results in Sect. 7.

2 Observations and data reduction

In this work, we observed a selected sample of disk-containing post-AGB nebulae through single-dish observations. The CO lines were studied using the IRAM-30 m telescope, while the SiO maser line was investigated using the Yebes-40 m telescope.

2.1 IRAM-30m data

Our data for the ¹²CO and ¹³CO J = 2-1 and J = 1-0 lines have been taken using the IRAM-30 m telescope using the Eight MIxer Receiver (EMIR), which can simultaneously observes in the 1.3 and 2.6 mm bands in dual polarization mode (EMIR, Carter et al. 2012). These observations were performed in two different runs under project 144-17 for 125 hours in total. Data from the first run were obtained between 31 January and 5 February 2018, while data from the second one were obtained between 10 May and 14 May 2018.

We connected the Fast Fourier Transform Spectrometer (FFTS) units to the EMIR receiver with a resolution of 200 kHz per channel, which corresponds to 0.25 and 0.5 km s⁻¹ in the 1.3 mm and 2.6 mm bands, respectively. The spatial resolution, or half power beam width (HPBW), of the observations is ~ 22″ and ∼11″ at the frequencies of CO J = 1-0 and J = 2-1 lines, respectively. Our sources were observed using the wobblerswitching mode, which provides flat baselines. The subreflector was alternatively placed between ON and OFF positions every 2 seconds, with the OFF reference being symmetrically located at ±90″ in the azimuth direction. We simultaneously observed using the vertically and horizontally polarized receivers, with no significant variations in their relative calibrations.

Pointing and focus were checked every four hours, with particular attention during the temperature variations associated with dusk and dawn. All corrections were negligible compared to the beam size. The flux calibration was derived using the chopper-wheel method, repeated at ∼15-minute intervals. Supplementary corrections were obtained from observations of the standard calibrators NGC 7027 and CW Leo (IRC+10216).

2.2 Yebes-40 m data

We conducted observations in the 7 mm band (Q band) using the Yebes-40 m telescope. Five sources from our sample were observed between 20 November 2024 and 13 January 2025, with a total on-source time of 50 hours under project code 24B011: HD213985, RVTau, UMon, HR4049, and IRAS 11472–0800. Q-band data for HD 52961 were obtained from project 20B006 and are already published. Detailed information on these sources is found in the next section.

The received signal was processed using the FFTS backend units. The Q band provided an instantaneous bandwidth of 18 GHz, a spectral resolution of 38 kHz, and an HPBW ranging from 37″ to 49″. Observations were carried out using the wobbler-switching technique, with a typical offset of around 300″ in the azimuthal direction, resulting in flat baselines. Spectra were simultaneously obtained in both vertical and horizontal linear polarizations, with no significant differences found in their flux density calibrations. For additional technical details, see Tercero et al. (2021).

2.3 Data reduction

The data reduction process was carried out with CLASS within the GILDAS¹ software package. For each source, we implemented the standard data reduction procedure, involving the exclusion of bad scans, averaging of valid ones (for both polarizations), and the subtraction of a first-order polynomial baseline. Also, the derived calibration corrections imply that the accuracy of the absolute scale is within a margin of <20%.

The observational data obtained from both the IRAM-30 m and Yebes-40 m telescopes are expressed in units of T^*_A, which represents the antenna temperature corrected for atmospheric attenuation, rearward spillover, and ohmic losses. This calibration provides a consistent measure of the received signal, but further conversion is required for physical interpretation. Thus, the data were converted to flux density (Jy) or the main-beam temperature (T_mb), depending on the analysis needs, using the telescope-specific efficiency parameters provided by IRAM² and Yebes³.

3 Source sample and CO data

Ten sources have been observed in this second single-dish CO survey project: HD 213985, TW Camelopardalis, RV Tauri, IRAS05208–2035, HD52961, UMonocerotis, HR4049, IRAS 11472–0800, V Vulpeculae, and AU Pegasi (see Table 1). This work provides a significant expansion in the number of disk-containing post-AGB sources observed in mm-waves through CO single-dish observations (see Bujarrabal et al. 2013, for details of the first single-dish CO survey).

These ten new sources are identified as binary post-AGB stars. Their evolved nature is evidenced by their high luminosities, which range from ∼70 to ∼12600 L_⊙ with a mean value of ∼2700 L_⊙ (private communication), and their intermediate spectral types (Oomen et al. 2019; Kluska et al. 2019). Additionally, all of them show significant NIR excess in their SED, together with far-infrared (FIR) excess that is indicative of material ejected by the star.

In Bujarrabal et al. (2013), 24 such circumbinary-disk-containing post-AGB nebulae were observed in CO lines; 3 of them are also observed in detail in this work. The 7 remaining new sources represent a 30% increase in the CO sample. Together, these 31 objects studied in CO lines represent ∼37% of the total sample of this type of disk sources, which currently numbers 85.

Nevertheless, we recall that the luminosities derived from SED fitting can sometimes suffer from substantial uncertainties, mainly due to errors in distance estimates and reddening corrections (see Kluska et al. 2022; Van Winckel 2025). This implies that some of the objects might actually fall below the canonical lower luminosity limit for AGB stars (∼2500L_⊙; Kamath et al. 2016), in which case they could be classified as post-RGB rather than post-AGB binaries (Moltzer et al. 2025). We therefore adopt the generic term post-AGB binaries throughout this paper for consistency.

We delve into the key findings obtained from the CO survey for each of these new sources:

• HD 213985: This source exhibits a narrow and faint ¹²CO J = 2-1 line profile. This tentative detection is centered at –29km s⁻¹ (see Fig. 1). Although this velocity differs from that reported by Oomen et al. (2018) (see Appendix B for a detailed discussion⁴), the molecular line profile suggests that the emission could still originate from circumstellar material physically linked to this source.

• TW Cam: we find no CO emission for this source. Strong ISM contamination is detected in the range –35 to –40 km s⁻¹.

• RV Tau: this source was previously observed with no detections. Our new and improved observations benefit from a longer integration time, allowing us to clearly detect a narrow ¹²CO J = 2-1 line (see Fig. 2).

• IRAS 05208–2035: this source does not show any observable CO emission.

• HD 52961: This object shows intense ¹²CO and ¹³CO lines (see Fig. 3). It exhibits a narrow line profile, characteristic of rotating disks, along with relatively wide line wings.

  Fig. 3 Newly detected CO line profiles in HD 52961. Axes as in Fig. 1.

  The J = 1-0 lines (both in ¹²CO and ¹³CO) show the distinctive double peak commonly seen in rotating disks. From both observational and theoretical perspectives, it is widely accepted that these profiles effectively indicate the presence of rotating disks (Bujarrabal et al. 2013; Gallardo Cava et al. 2021). The additional intensity in the line wings appears difficult to account for solely based on disk emission, so the presence of a relatively massive outflow could explain this CO line profile (as in the case of R Sct or AI CMi, see Bujarrabal et al. 2013).

• U Mon: the first CO survey did not yield any detections for this source. The superior quality of these recent observations has facilitated the clear detection of the ¹²CO J = 2-1 line. We are able to tentatively discern the double peak in this narrow CO line, which does not appear to show line wings (see Fig. 4), in a similar way to that found in AC Her (see Bujarrabal et al. 2013; Gallardo Cava et al. 2021).

• HR 4049: this source was already observed in the first survey, yielding weak CO line detections. Our new observations show significant improvement over the previous study: the signal-to-noise ratios for the ¹²CO J = 2-1, ¹²CO J = 1-0, and ¹³CO J = 1-0 transitions are now 8, 2.5, and 3.4, representing improvement factors of 1.5, 1.3, and 1.3, respectively. The resulting higher-quality data, particularly for the ¹²CO J = 2-1 and ¹³CO J = 1-0 lines, more useful for our study, lead to a much more accurate outcome (see Fig. 5).

• IRAS 11472–0800: this source also shows the narrow CO line profile of rotating disks, where the double peak can be tentatively seen (see Fig. 6). Note that the x-axis shows a velocity width of 200 km s⁻¹.

• V Vul: there is no detected CO emission from this source.

• AU Peg: no CO emission was detected from this source.

As a summary of what has been discussed, we have detected CO emission in HD213985, RVTau, UMon, HR4049, HD52961, and IRAS 11472–0800. Only HD52961 presents intense ¹²CO and ¹³CO lines compared with the rest of the surveyed sources and in line with other intense objects of our first CO survey, such as R Sct or the Red Rectangle. All of these CO-emitting sources show a narrow component in their line profiles, suggesting the presence of a rotating structure in their nebulae (see Figs. 1, 2, 3, 4, 5, and 6). For the sources without CO emission, we have estimated upper limits.

The main observational parameters are given in Table 2. The table lists the peak flux, I [Jy], its associated noise, σ (Jy), the integrated line intensity, ∖ I dV (Jy km s⁻¹), its associated uncertainty, σ(∖I dV) (Jy km s⁻¹), the spectral resolution, δV (km s⁻¹), the velocity centroid of the CO line measured with respect to the local standard of rest, V_LSR(km s⁻¹), and some relevant comments.

Uncertainty in the integrated line intensity for CO detections was estimated using the following expression: $$\sigma (\int I\text{d}V)=\sigma \delta V\sqrt{n_{\text{ch}}},$$(1)

where n_ch is the number of channels covered by the detected CO line. For the undetected lines, we used an upper limit of 3σ × √n_ch × δV, where n_ch corresponds to the number of channels typically spanned by the ¹²CO J = 2-1 line. When a CO line is detected, n_ch is taken as the number of channels covered by that line in the source. In contrast, when no line is detected, n_ch is taken as the average number of channels measured from all the detected lines in this work. The number of channels for this line is 52 for HD 52961, 7 for HD 213985, 11 for RV Tau, 18 for U Mon, 19 for HR 4049, and 8 for IRAS 11472–0800, with an average of 20 channels across all sources.

Fig. 1 Newly detected 12CO J = 2-1 line profile in HD213985. The x-axis indicates velocity with respect to VLSR and the y-axis represents the detected flux measured in janskys.

Fig. 2 Newly detected 12CO J = 2-1 line profile in RV Tau. Axes as in Fig. 1.

Fig. 4 Newly detected 12CO J = 2-1 line profile in U Mon. Axes as in Fig. 1.

Fig. 5 Newly detected CO line profiles in HR 4049. Axes as in Fig. 1.

Fig. 6 Newly detected 12CO J = 2-1 line profile in IRAS 11472-0800. Axes as in Fig. 1.

4 Total mass of circumbinary-disk-containing post-AGB nebulae

This section describes the mathematical formalism adopted for calculating the masses of CO-emitting regions using observational data. Subsequently, we shall discuss and interpret this data to extract physical information from the observed objects.

4.1 Calculation methodology

The methodology used in this work to estimate masses of diskcontaining nebulae is the same to that used in Bujarrabal et al. (2013); Gallardo Cava et al. (2022b). First, the main-beam temperature (T_mb) is defined as: $$T_{\text{mb}}(V_{\text{LSR}})=S_{\nu }\tau \varphi (V_{\text{LSR}}){\mathrm{\Omega }}_{\text{S}}{\mathrm{\Omega }}_{\text{mb}}^{-1},$$(2)

where each term is defined as follows:

• Ω_S is the solid angle subtended by the source, expressed as $${\mathrm{\Omega }}_{\text{S}}=2\pi (1-\cos \alpha )\approx \pi {\alpha }^2\approx \pi \frac{r^2}{d^2}\approx \frac{s}{d^2},$$(3)

  where α is the angle formed by the distance, d, and the radius, r, of the observed source. Considering that d will always take large values due to long distances, α will always have small values. Thus, we can approximate $\tan \alpha \approx \frac{r}{d}$ and $\cos \alpha \approx 1-\frac{{\alpha }^2}{2}$. The quantity s is the surface area of the observed source derived from its radius and its distance.

• Ω_mb is the solid angle subtended by the telescope main beam: $${\mathrm{\Omega }}_{\text{mb}}=\frac{\pi }{4\ln 2}{\theta }_{\text{HPBW}}^2\approx 1.133{\theta }_{\text{HPBW}}^2,$$(4)

  in which θ_HPBW is the HPBW.

• Φ (V_LSR) represents the normalized line profile in terms of Doppler-shifted velocity and it takes a value of 1 when it is integrated in terms of velocity.

• S_ν is the source function and it is given by $$S_{\nu }=\frac{h\nu }{k_{\text{B}}}(\frac{1}{\exp \frac{h\nu }{k_{\text{B}}{T}_{\text{ex}}}-1}-\frac{1}{\exp \frac{h\nu }{k_{\text{B}}{T}_{\text{BG}}}-1}),$$(5)

  where k_B is the Boltzmann constant, T_ex represents the excitation temperature of the molecular transition, and T_BG is the background temperature, with a constant value of 2.73 K.

• The characteristic optical depth, τ, is $$\tau =\frac{c^3}{8\pi {\nu }^3}{A}_{\text{ul}}{g}_{\text{u}}(x_{\text{l}}-x_{\text{u}})n_{\text{T}}XL,$$(6)

  where X represents the abundance of the molecule relative to the total number of particles (¹²CO and ¹³CO in this work), n_T is the total number density of particles, and L represents the typical length of the studied source along the line of sight. Thus, this opacity expression can be rewritten in terms of the mass, M, taking into account that M = n_T Ls = n M_mol L s, with M_mol being the total mass of the molecule-rich emitting gas and n the number of particles. A_ul is the Einstein coefficient of the molecular transition, where u and l are the upper and lower levels of the transition, respectively. The statistical weight of the upper level (g_u) is $$g_{\text{u}}\equiv g_J=2J+1.$$(7)

  Upper and lower level populations can be expressed as $$x_J\sim \frac{\exp (-E_J{T}_{\text{rot}}^{-1})}{F(T_{\text{rot}})},$$(8)

  in which the energy of the level is calculated assuming a simple linear rotor and taking into account the rotational constant of the molecule (B_rot is around 2.8 and 2.6 K for ¹²CO and ¹³CO, respectively): $$E_J=B_{\text{rot}}J(J+1).$$(9)

  The rotational temperature (T_rot) is the typical temperature value that represents the average of the excitation temperatures of the relevant rotational transitions used, i.e., low-J transitions in our case. This temperature is used to estimate the partition function: $$F(T_{\text{rot}})=\sum _J{g}_J\exp (-\frac{E_J}{T_{\text{rot}}})=\frac{T_{\text{rot}}}{B_{\text{rot}}}.$$(10)

  Finally, integrating Eq. (2) in terms of velocity, we get $$\begin{array}{rl}\int T_{\text{mb}}& \left(V_{\text{LSR}}\right)\text{d}V\\ & =S_{\nu }\frac{c^3}{8\pi {\nu }^3}{A}_{\text{ul}}{g}_{\text{u}}(x_{\text{l}}-x_{\text{u}})X\frac{4\ln 2}{\pi }{\theta }_{\text{HPBW}}^2{M}_{\text{mol}}{d}^{-2}.\end{array}$$(11)

Therefore, by following this procedure, we obtain an analytical expression that directly relates the area (velocity integrated main beam temperature) of the molecular rotational line to the nebular mass of the CO-emitting region, additionally being inversely proportional to the square of its distance: $\text{Area}\equiv \int T_{\text{mb}}(V_{\text{LSR}})\text{d}V\propto M_{\text{mol}}{d}^{-2}$.

The application of this procedure relies on the assumption of optically thin emission, a valid consideration in our case due to the notably weak emission of the molecular lines. Other factors can influence the opacity, such as the inclination angle of the disk. In systems where the disk is seen edge-on, the CO lines tend to be more optically thick, particularly near the systemic velocity. As in the previous CO single-dish work, we assume T_rot = 50 K, which is expected for kinetic temperatures between 70 and 100 K and total densities between approximately 10⁴ and 10⁵ cm⁻³ (see details in Bujarrabal et al. 2013). The gas mass values were derived by assuming an average particle mass of 3 × 10⁻²⁴ g that takes into account molecular hydrogen (H₂), atomic hydrogen, and helium. We also assumed a relative abundance of X(¹²CO) = 2 × 10⁻⁴, which may be a low value when compared to AGB CSEs (2-4 × 10⁻⁴ in O-rich stars and 6-8 × 10⁻⁴ in C-rich stars, see Teyssier et al. 2006; Ramstedt et al. 2008), but is in line with previous works on this type of objects and with the molecular study performed on these sources (see e.g. Bujarrabal et al. 2013; Gallardo Cava et al. 2021). We also considered a ¹²CO/¹³CO abundance ratio of 10, which is also the typical ratio found in these objects (see e.g. Bujarrabal et al. 2016, 2017, 2018; Gallardo Cava et al. 2021, 2023). When possible, we used the ¹³CO J = 1-0 line, which is the best for our purposes, because this transition is expected to be optically thinner than the J = 2-1 and those of ¹²CO. When only ¹²CO J = 2-1 was detected, the nebular mass was estimated based on this molecular line. Considering that this line is optically thicker than ¹³CO J = 1-0, the calculated mass values are underestimated, but can be corrected, as is discussed in Sect. 4.2. The assumption that the ¹²CO line is optically thicker than that of ¹³CO is valid because, for all sources without ¹³CO detection, the $I_{{}^{12}\text{CO}}^{\text{peak}}/3{\sigma }_{{}^{13}\text{CO}}$ ratio is always much lower than the assumed CO isotopic ratio of 10 (i.e., ¹²CO/¹³CO).

For the sources where both lines are detected, we find a ¹²CO/¹³ CO integrated line ratio with a mean value of about 2 for the J = 1-0 transition, and 6.6 for the ratio between the ¹²CO J = 2-1 and the ¹³CO J = 1-0 lines. These values are consistent with those found in other binary post-AGB disks (Gallardo Cava et al. 2022b), but they are lower than the ratios typically measured in other post-AGB stars. This suggests that our objects may present a higher ¹³CO abundance, with the binary nature of the sources likely being the origin of such unusually low ¹²CO/¹³CO ratios.

The derived H₂, atomic hydrogen, and helium masses correspond exclusively to the gas regions where CO emission is detected, as CO serves as our tracer for the total gas content. Contributions from other molecules are assumed to be negligible, since SiO, HCN, SiS, and CS have very low abundances compared to molecular hydrogen, on the order of 10⁻⁸−10⁻⁹, 10⁻⁹−10⁻¹⁰, 10⁻¹¹, and 10⁻¹⁰, respectively (see Gallardo Cava et al. 2022b, for complete description).

4.2 The derived nebular masses

Based on our CO data (see Table 2), we can estimate the total mass of the pPNe around our binary post-AGB stars using Eq. (11). We have compiled in Table 3 the values of the nebular masses derived using each of the rotational lines observed in this work: $M_{{}^{12}\text{CO}J=2-1}$, $M_{{}^{12}\text{CO}J=1-0}$, $M_{{}^{13}\text{CO}J=2-1}$, and $M_{{}^{13}\text{CO}J=1-0}$, according to ¹²CO J = 2-1, ¹²CO J = 1-0, ¹³CO J = 2-1, and ¹³ CO J = 1-0, respectively. These calculations have been made for both disk sources with molecular CO detection and those without, for which we present upper limits.

We begin by presenting the results obtained for the sources HD 52961 and HR 4049. These were derived directly from the ¹³CO J = 1-0 transition, ensuring that the derived masses are more accurate and representative of the total content in these circumbinary disk-containting post-AGB nebulae.

HD 52961 presents a nebular mass of 2.3 × 10⁻² M_⊙. This source show line profiles with a narrow component and relatively wide wings (see Fig. 3). These wings are relatively intense in ¹²CO J = 2-1 and cover a velocity-range of 20 km s⁻¹. We would need high-resolution maps and accurate models to determine the nebular structure and the contribution of the disk to the nebula, allowing us to classify it as either disk-dominated or outflow-dominated.

HR4049 presents a total nebular mass of 6.7 × 10⁻³ M_⊙. This source also exhibits the typical narrow line profiles of rotating gas, further displaying the characteristic double peak of Keplerian disks in its ¹²CO J = 2-1 line profile (see Fig. 5).

In the case of HD213985, RVTau, UMon, and IRAS 11472–0800, we have only detected the ¹²CO J = 2-1 line. This transition tends to be optically thick (specially in disk-containing nebulae around binary post-AGB stars), so the derived nebular mass would be underestimated. Overall, those sources of this work and the previous one (Bujarrabal et al. 2013) with well-detected ¹²CO and ¹³CO lines adhere to the following empirical law: $$\beth =\frac{M_{{}^{13}\text{CO}J=1\text{--}0}}{M_{{}^{12}\text{CO}J=2\text{--}1}}\approx \mathrm{13.6.}$$(12)

Thus, we can provide a more accurate nebular mass values for these four sources based on their assumed ¹³CO J = 1-0 emission, the optically thinnest line (in a similar way of Bujarrabal et al. 2013). In this way, the masses derived from the assumed ¹³CO J = 1-0 emission, $M_{{}^{13}\text{CO}J=1\text{--}0}^{\text{EL}}$, must be $\beth$ times those derived from the ¹² CO J = 2-1 line: $$M_{{}^{13}\text{CO}J=1\text{--}0}^{\text{EL}}=\beth \times M_{{}^{12}\text{CO}J=2\text{--}1}.$$(13)

As a result, we present the nebular masses of these objects: HD213985 has a mass of 2.5 × 10⁻⁴ M_⊙, RVTau has 2.5 × 10⁻⁴ M_⊙, UMon has 1.8 × 10⁻⁴ M_⊙, and IRAS 11472–0800 has 2.3 × 10⁻² M_⊙. These results are good estimates, since $M_{{}^{13}\text{CO}J=1-0}^{\text{EL}}\approx M_{{}^{13}\text{CO}J=1\text{--}0}$ (see Table 3). These four sources show narrow ¹²CO J = 2-1 line profiles of just ∼10kms⁻¹ in range that strongly suggest that the emission from their rotating disks dominate their respective nebulae (see Figs. 1, 2, 4, and 6). High-resolution interferometric observations are necessary to resolve the disk and disk structure. Using the single-dish data and analyses, we estimate only approximate disk masses.

For the non-CO-detected sources, we can estimate upper limits to the nebular mass by using the respective upper limits on the integrated line intensities (see Table 2). In a very similar way to those sources detected only in ¹²CO J = 2-1, we can estimate their mass upper limit using this same transition, while applying the $\beth$ correction factor for these objects. This approach yields a more stringent constraint on their nebular mass compared to directly using the upper limit derived from the ¹³CO J = 1-0 line. We think that rescaling the ¹²CO J = 2-1 result is a more restrictive approach and likely provides a more stringent estimate of the actual mass. In this way, we find that TW Cam has an upper limit to its nebular mass of 4.9 × 10⁻⁴ M_⊙, IRAS05208–2035 of 3.4 × 10⁻⁴ M_⊙, VVul of 2.2 × 10⁻⁴ M_⊙, and AU Peg of 1.4 × 10⁻⁴ M_⊙. These nebular mass upper limits are given in Table 3. The root mean square in the spectra is sufficiently low thanks to the depth of these observations. Despite the non-detection of CO, the nebular mass limit values are considered reliable and consistent with results derived from positive detections.

According to what has been discussed, the most accurate nebular mass values for these CO-emitting circumbinary-disk-containing sources can be found in Table 4. In addition, leveraging the high sensitivity of the observations, upper mass limits for disk objects without molecular emission are listed, and these values enhance the derived nebular mass results.

5 Diversity and trends in nebular masses

In this section, we present a comprehensive analysis of all available nebular mass determinations, including both our new measurements and those from previous studies. Our mass determinations, based on both single-dish and interferometric data (Bujarrabal et al. 2013, 2016, 2017, 2018; Gallardo Cava et al. 2021, 2022a, 2023, among others), incorporate two key improvements: (1) updated distance measurements from Bailer-Jones et al. (2021), and (2) new mass estimates for sources in Bujarrabal et al. (2013) that previously lacked CO detections and upper limits, calculated using the same methodology we apply to non-detections in this study (Sect. 4.1). Complete details of all mass determinations from previous works are provided in Appendix A.

5.1 Direct measures

The full sample of 31 sources observed at millimeter wavelengths to estimate their nebular masses is shown in Fig. 7. These histograms illustrate the wide variety of results, including both secure detections and upper limits derived from sources without observed molecular emission.

It is remarkable that this group of disk-containing post-AGB binaries, seemingly similar in nature, displays such a pronounced spread in nebular mass, or in some cases a complete lack of detectable molecular content. This contrast is evident even within subsamples: both disk- and outflow-dominated sources show mass differences exceeding a factor of 20. When the entire sample is considered, including objects with unknown disk-to-outflow mass ratios, the disparity becomes even more dramatic, surpassing a factor of 600.

Beyond this broad dynamic range, the mass distribution of CO-emitting sources is strongly right-skewed. The mean nebular mass is μ̄_det = 2.3 × 10⁻² M_⊙, which is significantly higher than the median value of μ̂_det = 1.3 × 10⁻² M_⊙, while the mode, μ^*_det = 8.7 × 10⁻³ M_⊙, closely reflects typical low-mass objects. This pronounced asymmetry, together with a large standard deviation of σ_det = 2.9 × 10⁻² M_⊙, reveals a heavy tail toward higher masses driven by a few extreme outliers (see Table 5).

In addition, we derived the geometric mean mass for the CO-emitting sources, obtaining μ_{g, det} = 8.0 × 10⁻³ M_⊙. This metric offers a more robust characterization of the sample by reducing the impact of extreme values. While the arithmetic mean is strongly affected by high-mass outliers, the geometric mean provides a more representative central tendency for the bulk of the nebular mass distribution, confirming that most systems cluster at lower masses than is indicated by the arithmetic mean. Based on the disk-mass percentage classification by Gallardo Cava et al. (2023), we find that disk-dominated sources present a geometric mean mass of 1.1 × 10⁻² M_⊙, while outflow-dominated sources exhibit a geometric mean mass of 2.9 × 10⁻² M_⊙, respectively. Taken together, these direct statistics highlight the strong diversity in nebular mass among post-AGB binaries, pointing to a wide range of physical conditions, or possibly rapid evolutionary tracks, within this class of systems.

Fig. 7 Top panel : histogram of the logarithm of disk masses for detected sources, overlaid with the probability density function (PDF; red line) and cumulative distribution function (CDF; dashed blue line) resulting from a censored normal likelihood fit that accounts for non-detections (upper limits). Shaded gray regions mark the ±1σlog and ±2σlog intervals around the fitted mean (vertical black line). The vertical dash-dotted line shows the geometric mean of the detections (green), while vertical dashed lines indicate the positions of the upper limits (in violet). Bottom panel : representation of individual sources, showing the log-mass of each object. Names appear to the left of each bar. Sources with upper limits are marked with arrows pointing left from their log-mass value. The shaded regions, vertical lines, and color coding are consistent with the top panel.

5.2 Normal distribution of nebular mass

To extend our mass distribution analysis, we incorporated upper limits by assuming a normal distribution for the logarithm of the masses (log-norm from now on), which is empirically justified by the relatively symmetric spread of log₁₀-transformed detected masses (see Fig. 7 top) and the consistent trend of upper limits with the observed distribution.

A log-normal distribution is well suited for this dataset, as the masses span multiple orders of magnitude and are strictly positive, effectively capturing the multiplicative nature of mass formation and the skewness of the data (Limpert et al. 2001). Our maximum likelihood approach incorporates the full dataset: for detections we include asymmetric mass uncertainties through distribution convolution, while for non-detections we properly account for upper limits using censored data analysis. This comprehensive treatment prevents bias and provides an accurate estimate of the underlying mass distribution for the entire population of CO-detected and non-CO-detected sources. In Fig. 7 top, we present the resulting log-normal fit.

This approach yields a log-normal distribution with parameters μ_log = −2.7 and σ_log = 1.0 (see Table 5). This fit reveals fundamental characteristics of the population that are not apparent from direct detections alone. The corresponding median mass, μ_log = 2.0 × 10⁻³ M_⊙, agrees well with the median of detected values, suggesting that it robustly represents the typical object in the overall disk-containing population. However, the log-normal model provides richer information on the shape of the distribution. The ±1σ_log interval, [1.8 × 10⁻⁴, 2.2 × 10⁻²] M_⊙, shows that 68% of sources fall within nearly an order of magnitude, while the broader ±2σ_log range, [1.6 × 10⁻⁵, 2.4 × 10⁻¹] M_⊙, highlights the asymmetry of the distribution, with a long tail toward higher masses. This explains why the arithmetic mean of detected sources, μ_det = 2.3 × 10⁻² M_⊙, is skewed upward relative to the median, as it is sensitive to massive outliers (see Table 5).

When comparing these results with the statistics from direct detections only, several important differences emerge. The detection-only arithmetic mean is notably higher than the median of the log-normal fit, reflecting how the inclusion of censored data affects population estimates. The geometric mean of detections (μ_{g, det} = 8.0 × 10⁻³ M_⊙), being less sensitive to extreme values, lies between the log-normal median and the detection mean, further highlighting the positive skewness of the distribution. The standard deviation of direct detections (σ_det = 2.9 × 10⁻² M_⊙), while large, underestimates the full spread implied by the log-normal standard deviation apparent in the sigma parameter of the log-normal fit and its associated confidence intervals. Most importantly, the detection-only analysis misses the true dynamic range of the population, particularly the substantial number of low-mass objects implied by the lower bounds of the log-normal fit.

These differences have important implications for interpreting the post-AGB disk population. The log-normal fit reveals that most objects have relatively low masses clustered around the median, along with a smaller subpopulation of significantly more massive disks that skew the arithmetic mean upward. This explains why the mean derived from detections alone exceeds both the median and the geometric mean, reflecting a bias toward more easily detectable, massive sources in the observational sample. By integrating both detections and upper limits into a statistically consistent framework, the log-normal model offers a more complete and accurate characterization of the underlying mass distribution. It shows that the majority of objects are less massive than direct detections suggest, and that relying solely on simple summary statistics underestimates both the spread and central tendency of the data. The consistency in median values across methods supports the robustness of this parameter, while the divergence in other metrics underscores the need for rigorous statistical treatment in studies of circumbinary disk-bearing pPNe. Ultimately, the results provide compelling evidence that the mass distribution in these systems is highly asymmetric, shaped by a substantial population of low-mass, often undetected, disks.

5.3 Uncertainties

We are well aware that there is a nebular mass problem. What we derive can be considered to be a lower limit to the total mass lost by the primary, as the derived values are very low when compared to the amount of mass the star must have lost to become a post-AGB star (however, this may not be the case for post-RGB primaries). This missing mass could reside in non-molecular components. In Gallardo Cava et al. (2023), the possible contributions from atomic material, the collimated low-density high-velocity jets launched by the companion during the late AGB and early post-AGB stages, and the frequently debated presence of very extended halos originating during the AGB phase, were analyzed. It was concluded that these non-molecular structures are unlikely to significantly contribute to the total nebular mass.

In our targets, one possibility can be that a significant fraction was accreted by the companion star during or after the massloss phase of the primary, potentially accelerating its evolution. This mass-loss problem is not exclusive to binary post-AGB stars. It is also present in many evolved nebulae (Santander-García et al. 2022), where we are not able to determine where this missing nebular mass resides.

The uncertainties associated with the calculated masses were determined by considering that the main sources of error are the absolute calibration uncertainty of the observations and the distance to each source. For the absolute calibration error, we found values of approximately 10% and 15% for the J = 1-0 and J = 2-1 transitions, respectively. Regarding distances, this parameter significantly contributes to the overall uncertainties due to the inverse-square dependence of the mass on distance. The asymmetry in their uncertainties propagates into similarly asymmetric errors in our results (see Table 4).

6 Searching for molecules other than CO

6.1 Observational results

The CO-emitting sources analyzed in this work have been extensively observed in the Q band using the Yebes-40 m telescope, as part of a dedicated effort to investigate their molecular content beyond the standard CO lines. Complementary observations were also carried out with the IRAM-30 m telescope during the CO campaigns.

The primary objective of the Yebes-40 m Q-band observations was the search for SiO maser emission, along with the detection of the SiO thermal line. Despite the depth and sensitivity of our observational setup, no SiO emission, either maser or thermal, was detected in any of the targeted sources. This result, however, is consistent with expectations: previous comprehensive surveys already indicated the molecular scarcity that characterizes these post-AGB nebulae hosting rotating disks (see Gallardo Cava et al. 2022b, for details).

Even in the absence of clear detections, the non-detections themselves are scientifically valuable. The stringent upper limits derived from these observations provide meaningful constraints on the presence of SiO. The derived upper limits for each source are summarized in Table C.1.

6.2 Lack of molecular content

Here, we discuss several possible reasons that might explain the low or apparent absence of molecular content in some of these circumbinary-disk-containing objects.

One possible reason for the apparent lack of CO content in these sources could be the presence of very low-mass rotating disks. In this case, future, deeper, and more sensitive observations might detect weak molecular emission, revealing these faint rotating disks. Additionally, the absence of CO emission in these sources could be attributed to photodissociation effects. Although our post-AGB stars are not hot enough efficiently photodissociate CO, the binary system can still destroy molecules in the inner regions of the nebula due to the UV radiation from the accretion disk around the companion or the interaction with the ejected ionized wind. Additionally, CO can be subject to freeze-out and interstellar-UV photodissociation in the outer and colder parts of the rotating disk. Another possibility for explaining this poor molecular content is that the rotating disk is not detected simply because it is too small. However, this scenario is unlikely, since for the disk to be both massive and undetected it would need to be optically thick, which would actually make it more detectable due to the higher densities.

The absence of molecular content can also be observed in protoplanetary disks. The aforementioned phenomena explaining the lack of CO have been widely confirmed in these disks around young stars (see Piétu et al. 2014; Barenfeld et al. 2017; Miotello et al. 2021, and references therein).

Another more plausible possibility is that CO escapes from the disk through outflows, while the dust, particularly large grains, remains gravitationally bound (Gallardo Cava et al., in prep). In this scenario, the CO non-detections may correspond to older disks that have experienced prolonged molecular gas depletion due to sustained outflow activity.

In summary, it would be possible to differentiate between very low-mass rotating disks and very high-density compact disks. For both cases, and even for those showing CO emission, the additional effect of photodissociation occurring in both the inner and outer regions of the disk must be considered.

Fig. 8 Distribution of 85 known Galactic binary post-AGB stars based on CO observations. The gray segment represents stars not observed in CO. Among those observed (encircled in blue), the dark blue sector corresponds to CO emitters, while the light blue sector shows nondetections.

7 Conclusions

To date, 85 Galactic post-AGB binary objects with circumbi-nary disks (studied through optical and infrared observations) are known (Kluska et al. 2019). In this work, we have observed ten sources in mm-waves through single-dish observations using the IRAM-30 m radiotelescope (see Table 1). These new sources significantly expand the total sample of objects observed in mm-wave in the search of CO. We note that three of the objects in this study were previously observed in the earlier survey; however, these new observations are much deeper and sensitive, leading to an increase in the percentage of objects studied at these wavelengths from 28% to 37%.

In this work, we have confirmed that six out of the ten studied objects do emit in CO: HD213985, RVTau, HD52961, UMon, HR4049, and IRAS 11472–0800. These six new detections significantly increase the sample of objects that emit in CO, rising from 18% to 25% of the 85 disk-containing post-AGB nebulae. Regarding the other four new objects that do not emit in molecules in a detectable way, TW Cam, IRAS 05208–2035, V Vul, and AU Peg, despite the absence of CO emission, we have achieved very promising upper limits. The main parameters derived from the observations presented in this work can be consulted in Table 2. Of the 31 objects studied in mm-waves (Bujarrabal et al. 2013, 2016, 2017, 2018; Gallardo Cava et al. 2021, 2023, and this work), 20 show CO emission, representing a 65% success rate in molecular detection (see Fig. 8).

From these new CO data, we calculated the derived nebular mass for each of the observed rotational lines (see Table 3). We observe that the mass increases as the studied line becomes less opaque. Therefore, we used the ¹³CO J = 1-0 data whenever possible. However, we find that for all disk-containing sources emitting in both ¹²CO and ¹³CO, the ratio of derived nebular masses from the ¹³CO J = 1-0 and ¹²CO J = 2-1 lines is 13.6. Thus, for sources with ¹²CO detections but without ¹³CO emission, we have provided a more accurate estimate of the total nebular mass corresponding to its presumed emission in this molecule, which is reasonably consistent with the upper limits obtained. For sources without molecular emission, we have calculated upper limits for the nebular mass for each of the observed rotational lines. Our best-estimated nebular mass values for this sample of ten objects can be found in Table 4.

Moreover, we have compared the nebular mass values from this work with all other existing values from previous studies using the most recent distances. According to Fig. 7, the data show a normal distribution when examined on a logarithmic scale with high variability, as these types of disk-containing sources exhibit a wide range of total nebular mass values (including both the rotating and expanding components). This dispersion is evident even in sources where the contribution of the disk to the total nebular mass is similar. However, overall, we can find mass differences that exceed a factor of 600. The full CO-studied dataset presented in Table 4 constitutes the most comprehensive and up-to-date compilation of nebular masses for post-AGB sources.

Additionally, we have presented a fit of our results to a lognormal distribution. This modeling approach incorporates the upper limits derived for the entire sample, recognizing their statistical and scientific value. However, this model exhibits a somewhat larger dispersion, as it inherently gives more weight to extreme values. Despite this limitation, the log-normal approach offers a significant advantage: it accounts for the rotating disks for which we only have upper limits to the mass, ensuring that these disk-containing sources are not excluded from the overall statistical analysis.

The typical mass of these objects, accounting for both detections and non-detections through censored data analysis, is approximately 2.0 × 10⁻³ M_⊙. This value corresponds to the median of the log-normal fit, making it a robust indicator that is not biased by the presence of high-mass outliers.

The CO-emitting sources analyzed in this work have been deeply observed in the Q band with the aim of detecting molecular species beyond CO, with particular focus on SiO maser and thermal emission. Despite the sensitivity of our observations, no SiO emission was detected in any of the targets. This lack of detection reinforces the idea that these disk-bearing post-AGB objects are severely deficient in molecular content, an intriguing characteristic whose origin remains an open question. A plausible explanation for the CO non-detections and lack of other molecules is the preferential removal of molecular gas through outflows, while the dust, particularly larger grains, remains gravitationally bound, leading to compact disk sources. This scenario implies that some of the undetected systems may be more evolved disks that have experienced prolonged gas depletion over time.

Finally, this work lays a valuable foundation for future research, as it paves the way for detailed investigations into the intrinsic nature of these systems. Although these systems appear remarkably similar at first glance, they exhibit substantial differences across various physical and observational properties. The only aspect we can say with certainty that they share is the presence of an evolved binary central star and a surrounding circumbinary disk.

Acknowledgements

We are grateful to the anonymous referee for the valuable recommendations and comments. This work is part of the coordinated research project NEBULAe WEB, grants PID2O19-1o52O3GB-C21 funded by MICIU/AEI/10.13039/501100011033, and coordinated project MESON, grant PID2023-146056NB-C21, funded by MICIU/AEI/10.13039/501100011033 and ERDF/EU. This work is based on observations carried out under project number 144-17 with the 30 m telescope. IRAM is supported by INSU/CNRS (France), MPG (Germany) and IGN (Spain). IRAM, the Institut de radioastronomie millimétrique, is an international research institute funded the French Centre National de la Recherche Scientifique (CNRS), the German Max-Planck Gesellschaft (MPG) and the Spanish Instituto Geográfico Nacional (IGN). This work is also based on observations carried out with the Yebes-40 m telescope (20B006). The 40 m radiotelescope at Yebes Observatory is operated by the Spanish Geographic Institute (IGN; Ministerio de Transportes y Movilidad Sostenible). This work has made use of GILDAS software for data reduction and analysis (https://www.iram.fr/IRAMFR/GILDAS)

References

Appendix A Additional nebular mass calculations

In this appendix, we present the derived mass values for those seven sources presented in Fig. 7 that were observed in Bujarrabal et al. (2013), and for which only the observational σ was provided without an estimation of the upper limit for the nebular mass. These objects were observed in the ¹²CO J = 2-1 and ¹²CO J = 3-2 lines with the APEX telescope.

Analogous to the method used for the sources in this work, we first calculate an upper limit for the integrated area using Eq. 1. In this case, the number of channels is again the average of the channels that typically covered the observed line, while the remaining values can be found in the original article. The method for calculating the mass is the same as described in Sect. 4.1. In this way, the masses derived from the ¹²CO J = 2-1 line should be corrected by multiplying them by a factor of 13.6 to match the masses derived from the assumed ¹³CO J = 1-0 emission.

The upper limits to the nebular mass for each of these sources, along with the upper limit of the integrated area of the ¹²CO J = 2-1 line used to calculate the mass, can be found in Table A.1.

Appendix B Heliocentric velocity comparison

In this appendix, we present a comparison between the heliocentric velocities derived from our CO spectral line observations and the systemic velocities previously reported (Oomen et al. 2018).

In most cases, the velocities derived from our CO observations are in excellent agreement with those reported by Oomen et al. (2018), showing very low dispersion (see Table B.1 and Fig. B.1). The mean μ_v and scatter σ_v of the residuals further confirm the strong consistency between both datasets. This agreement reinforces the reliability of the velocity measurements in tracing the systemic motion of post-AGB objects.

One notable exception is HD 213985, for which the COderived velocity significantly deviates from the value listed by Oomen et al. (2018). The ¹²CO J = 2-1 line profile observed for this source is very narrow and symmetric, characteristic of emission from a Keplerian rotating disk (see Fig. 1). We are confident that the emission does not originate from ISM features. In typical ISM conditions, the lower-energy ¹²CO J = 1-0 transition is significantly stronger than ¹²CO J = 2-1 (e.g., Falgarone et al. 1998). In our observations, the ¹²CO J = 1-0 line is not detected, while the ¹²CO J = 2-1 line is clearly present, which is incompatible with an ISM origin and supports a circumstellar disk scenario.

Fig. B.1 Comparison between heliocentric velocities derived from CO observations, Vhel,CO, and those reported by Oomen et al. (2018), Vhel. Each point corresponds to an individual disk-containing source, with the color indicating the residual velocity difference (Vhel,CO - Vhel). The dashed line represents the one-to-one relation (Vhel,CO = Vhel). The average offset μV and scatter σV of the residuals are shown in the bottom-right box.

Appendix C Yebes-40m and IRAM-30m results

In this appendix, we present the results for the most relevant lines for those sources of this work in which CO lines have been detected: SiO (both thermal and maser emission), ²⁹SiO,³⁰SiO, CS, SO, C¹⁷O, and C¹⁸O. The results are summarized in the following table (see Table C.1):

All Tables

All Figures

	Fig. 1 Newly detected 12CO J = 2-1 line profile in HD213985. The x-axis indicates velocity with respect to VLSR and the y-axis represents the detected flux measured in janskys.
In the text

	Fig. 2 Newly detected 12CO J = 2-1 line profile in RV Tau. Axes as in Fig. 1.
In the text

	Fig. 3 Newly detected CO line profiles in HD 52961. Axes as in Fig. 1.
In the text

	Fig. 4 Newly detected 12CO J = 2-1 line profile in U Mon. Axes as in Fig. 1.
In the text

	Fig. 5 Newly detected CO line profiles in HR 4049. Axes as in Fig. 1.
In the text

	Fig. 6 Newly detected 12CO J = 2-1 line profile in IRAS 11472-0800. Axes as in Fig. 1.
In the text

Fig. 7 Top panel : histogram of the logarithm of disk masses for detected sources, overlaid with the probability density function (PDF; red line) and cumulative distribution function (CDF; dashed blue line) resulting from a censored normal likelihood fit that accounts for non-detections (upper limits). Shaded gray regions mark the ±1σlog and ±2σlog intervals around the fitted mean (vertical black line). The vertical dash-dotted line shows the geometric mean of the detections (green), while vertical dashed lines indicate the positions of the upper limits (in violet). Bottom panel : representation of individual sources, showing the log-mass of each object. Names appear to the left of each bar. Sources with upper limits are marked with arrows pointing left from their log-mass value. The shaded regions, vertical lines, and color coding are consistent with the top panel.

In the text

	Fig. 8 Distribution of 85 known Galactic binary post-AGB stars based on CO observations. The gray segment represents stars not observed in CO. Among those observed (encircled in blue), the dark blue sector corresponds to CO emitters, while the light blue sector shows nondetections.
In the text

Fig. B.1 Comparison between heliocentric velocities derived from CO observations, Vhel,CO, and those reported by Oomen et al. (2018), Vhel. Each point corresponds to an individual disk-containing source, with the color indicating the residual velocity difference (Vhel,CO - Vhel). The dashed line represents the one-to-one relation (Vhel,CO = Vhel). The average offset μV and scatter σV of the residuals are shown in the bottom-right box.

In the text