Open Access
Issue
A&A
Volume 665, September 2022
Article Number A45
Number of page(s) 11
Section Planets and planetary systems
DOI https://doi.org/10.1051/0004-6361/202244104
Published online 08 September 2022

© L. M. Pirovano et al. 2022

Licence Creative CommonsOpen Access article, published by EDP Sciences, under the terms of the Creative Commons Attribution License (https://creativecommons.org/licenses/by/4.0), which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.

This article is published in open access under the Subscribe-to-Open model. Subscribe to A&A to support open access publication.

1 Introduction

Water is a key element in the physical and chemical evolution of protoplanetary discs and in the formation of planets (e.g. van Dishoeck et al. 2021). Ro-vibrational and rotational transitions of H2O and OH in discs have been found by several authors using Spitzer (e.g., Carr & Najita 2008; Salyk et al. 2008; Pontoppidan et al. 2010), Keck/NIRSPEC (e.g. Salyk et al. 2008; Mandell et al. 2008; Doppmann et al. 2011), VLT/CRIRES (e.g., Fedele et al. 2011; Banzatti et al. 2017), VLT/VISIR (Pontoppidan et al. 2010), and Berschel (e.g., Hogerheijde et al. 2011; Riviere-Marichalar et al. 2011; Meeus et al. 2012; Fedele et al. 2012, 2013a; Podio et al. 2013). The detection of H2O and OH triggered intensive modelling campaigns designed to improve our understanding of the origin and abundance of water in discs. Various chemical and physical-chemical models have been developed, including different processes: Glassgold et al. (2009) and Najita et al. (2011) focus on the formation of H2O in the inner disc and take into account the effect of the energetic radiation (Ultraviolet and X-rays); Bethell & Bergin (2009) and Ádámkovics et al. (2014) included molecular self-shielding to explain the high abundance of H2O in the inner region of T Tauri discs and the physical-chemical models of Woitke et al. (2009); Bruderer et al. (2012); Du & Bergin (2014) include gas-phase H2O formation in the entire disc and in the disc atmosphere including desorption from icy grains.

Observations of multiple transitions of H2O in protoplanetary discs have the potential to unveil the abundance and distribution of water vapour – both warm and cold – in discs at the time of planet formation, as shown by Zhang et al. (2013) for the disc around TW Hya. These authors find a high concentration of water vapour at r ~ 4 au, corresponding to the water condensation front (or snow line). This latter experiment was possible thanks to the combination of low- and high-J H2O transitions from Berschel and Spitzer.

This paper presents a H2O ‘line mapping’ of the two Herbig Ae systems where rotational H2O emission has been detected with Berschel, namely HD 100546 and HD 163296. Both discs show substructures (cavities and dust rings) in the dust continuum with ALMA, which may or may not trap icy pebbles in the outer disc that could otherwise drift inward. Our analysis is based on a mixture of H2O lines targeted with the Photodetector Array Camera and Spectrometer (PACS; Poglitsch et al. 2010) and the Heterodyne Instrument for the Far Infrared (HIFI; de Graauw et al. 2010) on Berschel. The aim of this paper is to determine the abundance distribution of H2O vapour in HD 100546 and HD 163296 by comparing the observed H2O spectrum to a set of physical–chemical models of discs.

2 Water reservoirs in discs

According to physical-chemical models (e.g. Woitke et al. 2009; Bruderer et al. 2012; Du & Bergin 2014), there are three main reservoirs of gaseous H2O in discs:

  • s1

    is the inner reservoir (inside the snow line) where H2O molecules are formed through gas-phase reactions. The high kinetic temperature in this region allows the energy barrier of the radical-H2 reactions to be overcome (e.g. Glassgold et al. 2009). At the same time, this reservoir is fully shielded from the ultraviolet (UV) radiation coming from the star, meaning that H2O molecules are protected from photodissociation. The H2O abundance in this region can be further enhanced if icy planetesimals migrate inwards through the snow line and H2O molecules are released via thermal desorption.

  • s2

    is the (cold) photodesorption layer in the outer disc where the H2O vapour is released from the icy grains. The gas temperature in this region is too low (<100 K) and radical H2 reactions are inhibited, meaning that the H2O vapour reservoir is only supplied by photodesorption.

  • s3

    is the warm reservoir which is located at intermediate layers above the disc midplane and that can extend to several astronomical units (au) from the star. The production of H2O molecules is driven by gas-phase reactions. UV photons are not fully absorbed and photodissociation of H2O molecules is relevant in this region: because of this, the amount of H2O vapour in this layer is also regulated by the self-shielding (e.g. Bethell & Bergin 2009).

The three H2O-rich regions are described in Sect. 4. Each region contributes to the H2O line flux, but while the high- J rotational lines and the ro-vibrational lines arise primarily from reservoirs s1 and s3, the low- J rotational transitions and in particular those connecting to the ground state, are more sensitive to the cold photodesorption layer s2.

3 Observations

3.1 PACS spectra

HD 163296 and HD 100546 were observed with Herschel/PACS as part of the DIGIT (Green et al. 2013) and GASPS (Dent et al. 2013) key programs (Table 1). The PACS spectra used here are from the DIGIT program and the observations are presented in Fedele et al. (2012, 2013a). The spectral window of Herschel/PACS (50-200 µm) includes several rotational H2O transitions ranging in upper level energy from Eup = 114 K (212–101 at 179.52 µm) to Eup = 1552 K (918−909 at 62.93 µm). The fluxes and upper µmits of the H2O transitions are reported in Fedele et al. (2012, 2013a). A number of transitions are detected towards HD 163296 (Fedele et al. 2012; Meeus et al. 2012), while no H2O emission is detected in the DIGIT and GASPS programs towards HD 100546 (Fedele et al. 2013a). The PACS spectrum of HD 100546 was taken early in the mission (Sturm et al. 2010) and has a lower signal-to-noise ratio (S/N) compared to that of HD 163296. The upper µmits reported by Fedele et al. (2013a) are almost an order of magnitude higher than the line fluxes measured for HD 163296.

A deep PACS observation of HD 100546 was executed as part of a calibration program targeting the H2O transition 321 − 212 at 75.38 µm. Figure 1 shows the pipeline-reduced spectrum (taken in ‘linespec’ mode). The CO J = 35−34 transition at 74.9 µm is detected along with two other emission lines at 75.01 µm and 75.41 µm. These two emission lines could be due to ro-vibrational transitions of NH3 or SO2 (further analysis is deferred to a forthcoming paper). Notably, the H2O 321 – 212 is not detected. The 3σ upper µmit (Table 2) is computed adopting a line FWHM of 0.019 µm as measured on the CO J = 35–34 transition and accounting for the flux calibration error provided by the Herschel pipeline.

thumbnail Fig. 1

Herschel/PACS deep integration of HD 100546 at 75 µm. The CO J = 35–34 line is detected along with two other emission lines likely due to NH3. The blue dashed line indicates the position of the H2O 321 – 212 transition.

Table 1

Observations log.

3.2 HIFI spectra

The H2O ground state transitions were observed with Herschel/ HIFI during an open call proposal (Table 1). The line fluxes and upper µmits were published by Du et al. (2017). Here we report a new, more stringent upper µmit to the 110 – 101 transition towards HD 163296 based on deeper integration (Table 1). We downloaded the level 2 data product from the data archive. The WBS spectrum is converted from antenna temperature to mean-beam temperature after correcting for the beam and forward efficiency (Roelfsema et al. 2012). Finally, the spectra of the horizontal and vertical polarization are averaged together. The reduced and flux-calibrated spectrum is presented in Fig. 2. The ground-state H2O line remains undetected even in this case. We compute the upper limit and converted from K km s−1 to W m−2 using the formula Fvdv=2.73×1029v3θ2Tvdv,$\int {{F_v}{\rm{d}}v = 2.73 \times {{10}^{ - 29}}{v^3}{\theta ^2}\int {{T_v}{\rm{d}}v,} } $(1)

where ν is the line frequency (GHz) and θ the half power beam width (arcsec). Following Du et al. (2017), we adopted a line full width at half maximum of 9 km s−1 (e.g. Thi et al. 2004), which is similar to the low- J CO lines, and the line upper µmits include a 20% flux calibration uncertainty.

Table 2

H2O line fluxes.

thumbnail Fig. 2

Herschel/HIFI spectrum of HD 163296110 − 101 (at 556 GHz) at rest-frame velocity.

4 Analysis

The analysis presented in this paper is based on simulations with the DALI thermo-chemical models of discs (Bruderer et al. 2012; Bruderer 2013).

4.1 Disc physical setup

DALI takes as input the stellar spectrum (in this case corresponding to the blackbody emission at T = Teff) and a power-law surface density structure with an exponential tail: gas(R)c(RRc)γexp[ (RRc)2γ ],${\sum _{{\rm{gas}}}}\left( R \right){\sum _{\rm{c}}}{\left( {{R \over {{R_{\rm{c}}}}}} \right)^{ - \gamma }}\exp \left[ { - {{\left( {{R \over {{R_{\rm{c}}}}}} \right)}^{2 - \gamma }}} \right],$(2)

where R is the radial distance from the star, Rc the critical radius, and Σc is the gas surface density at R = Rc. The dust surface density is Σgasgd, where Δgd is the gas-to-dust-mass ratio. The density in the vertical direction is assumed to follow a Gaussian distribution with scale height h = H/R: h=hc(rRc)ψ,$h = {h_{\rm{c}}}{\left( {{r \over {{R_{\rm{c}}}}}} \right)^\psi },$(3)

where hc is the scale height at R = Rc and ψ the flaring angle. The dust grain size distribution is described using two grain size populations as in D'Alessio et al. (2006), and dust mass absorption cross sections from Andrews et al. (2011). The grain size distribution is approximated by a power-law distribution with exponent q (equal for both populations). The settling of the large grains is controlled by the parameters χ (scale height of the large grain population compared to the gas one) and flarge (small-to-large grain mass ratio).

DALI solves the 2D dust continuum radiative transfer and determines the dust temperature and radiation field strength at each disc position. In the standard procedure, DALI iteratively solves the gas thermal balance and chemistry. The H2O spectrum is obtained by computing the non-local thermodynamic equilibrium (NLTE) level populations, including infrared pumping and by fixing the ortho-to-para ratio to 3. The collisional rate coefficients are from the LAMDA database (Schöier et al. 2005, and references therein).

4.1.1 The spectral energy distribution

The goal of this first step is to constrain the disc dust temperature and density gradients using the spectral energy distribution (SED) as an indicator. A reference physical model is obtained from the comparison of the observed SEDs with a grid of synthetic SEDs based on a grid of DALI models.

The initial setup for the parameterisation of HD 163296 is taken from Kama et al. (2020). From the reported disc gas mass of 6.7 × 10−2 M we derive a value of = 6.8 g cm−2 at Rc = 125 au, which was obtained by solving the mass integral. The stellar parameters have been changed to L = 17 L and M = 2 M, adapted from Vioque et al. (2018). With this setup, a first grid of simulations changing Σc, Δgd, hc, and ψ was run to determine the disc reference model. The SED is shown in Fig. 3 (lower left panel).

In the case of HD 100546, the initial setup is taken from Kama et al. (2016). The initial values are Σc = 18.2 g cm−2 at Rc = 50 au and Δgd = 40 which yield a total gas mass of 3.2 × 10−2 M and total dust mass of 8.1 × 10−4 M. The Kama et al. (2016) setup was derived using a distance of 97 pc; the stellar and disc parameters were updated to account for the new distance estimate of d = 108 pc (Gaia Collaboration 2021). The setup includes an inner dust gap between 4 au and 13 au. Recently, Fedele et al. (2021) reported a new analysis of ALMA 870 dust continuum observations that reveals the presence of an outer gap in the disc extending from ~40 to ~150 au (hinted at by Walsh et al. 2014) that divides the disc inner bright ring, centred at nearly 28 au, and the disc outer faint ring, centred at ~200 au. To reproduce the observed double-ring structure, the density structure was modified using different scaling factors (δ) for the gas and the small and the large dust grains to reduce the densities in the inner dust cavity and in the outer gap: ni={ ni×δi,cvif Rcav,in<R<Rcav,outniif Rcav,outRRgap,in ni×δi,gapif Rgap,in< R<Rgap,out niif rRgap,out, ${n_i} = \left\{ {\matrix{ {{n_{\rm{i}}} \times {\delta _{{\rm{i,cv}}}}} \hfill &amp; {{\rm{if}}\,{R_{{\rm{cav,in}}}} &lt; R &lt; {R_{{\rm{cav,out}}}}} \hfill \cr {{n_{\rm{i}}}} \hfill &amp; {{\rm{if}}\,{R_{{\rm{cav,out}}}} \le R \le {R_{{\rm{gap,in}}}}{\rm{ }}} \hfill \cr {{n_{\rm{i}}} \times {\delta _{{\rm{i,gap}}}}} \hfill &amp; {{\rm{if}}\,{R_{{\rm{gap,in}}}} &lt; \,R &lt; {R_{{\rm{gap,out}}}}{\rm{ }}} \hfill \cr {{n_{\rm{i}}}} \hfill &amp; {{\rm{if }}r \ge {R_{{\rm{gap,out}}}}} \hfill \cr } ,} \right.$(4)

where i = gas, small grains, large grains. The values of all the δ factors are reported in Table 3, while the obtained gas and dust densities are shown in the top left panels of Fig. 4. The small dust grains are assumed to be dynamically coupled to the gas, and so we impose δgas = δsmall at every radius, while δlarge,cave = δlarge,gap = 0 to clear the inner cavity and the outer gap from the large dust grains. ALMA observations of CO isotopo-logues suggest a gas density drop of one order of magnitude in the outer gap (Booth private communication) and so we reduced the gas density by the same amount. With this density setup, a grid of simulations was created varying Σc, Δgd, hc, and ψ to match the observed SED.

We note that the disc around HD 163296 is also characterised by dust gaps and rings (e.g. Isella et al. 2016; Andrews et al. 2018), but in contrast to HD 100546, the gaps are relatively narrow and their presence does not seem to affect the H2O distribution (see discussion further below).

thumbnail Fig. 3

DALI disc structure for HD 163296. Top panels: gas and dust density structure, the gas temperature, and the H2O abundance structure. The three main H2O reservoirs are indicated by the labels s1, s2, and s3. Bottom panels: model predictions of the SED (left), CO ladder (middle), and H2O fluxes (right). The observed fluxes and 3σ upper µmits are shown as (grey) circles and open triangles, respectively. The representative model shown as a red line and red squares) reproduces the SED (left) and the CO rotational ladder (middle) well but is not able to reproduce the H2O line fluxes and upper µmits. We note that the H2O abundances shown here are those of the representative disc model based on the full gas-grain chemistry calculation.

thumbnail Fig. 4

Same as Fig. 3 but for HD 100546.

Table 3

DALI model parameters.

4.1.2 CO rotational ladder

In a second step, we modelled the CO rotational ladder to constrain the 2D gas temperature structure: the fluxes of the optically thick CO rotational transitions are indeed a powerful tool to trace the gas temperature throughout the radial and vertical extent of the disc (Fedele et al. 2013b, 2016). Observational data of the CO rotational emission in HD 163296 are taken from Fedele et al. (2016), while data of HD 100546 are taken from Meeus et al. (2013); van der Wiel et al. (2014); Fedele et al. (2016). The CO collisional rate coefficients are from Yang et al. (2010).

To model the CO rotational ladders, a grid of DALI models was created starting from the parameters that best fit the SED and by changing parameters that mostly affect the gas structure, namely: Σc, Δgd, and ψ. These parameters mostly change the shape of the CO rotational ladder; Σc and Δgd are changed simultaneously to lower the total gas mass, maintaining the disc dust mass at a constant level, thus fixing the shape of the SED. The parameters of the reference models are given in Table 3 and the synthetic CO ladders are presented in Figs. 3 and 4 (middle panels) for the two discs, respectively.

4.2 The H2O line fluxes

Figures 3 and 4 (lower right panel) show the flux of the H2O rotational transitions obtained with the reference thermo-chemical models compared to the observational data. The DALI chemical network used here is based on the UMIST06 database (Woodall et al. 2007) and contains ten elements (H, He, C, N, O, Mg, Si, S, Fe, and PAH), 109 species, and 1463 reactions. The chemical network is initialised with atomic abundances, with a carbon-to-oxygen abundance ratio of 0.468.

In the case of HD 163296, the reference thermo-chemical model (red squares) almost matches the low flux of the ground-state transitions but it underestimates the flux of the excited transitions by up to two orders of magnitude. On the other hand, the HD 100546 reference model predicts a an overly strong emission in the 321 –212 transition, while it almost matches the flux of the ground-state transitions.

The mismatch between the observed and synthetic H2O line fluxes cannot be due to the global physical properties, such as the density and temperature structures, as these are constrained by the SED and by the CO ladder. Similarly, the discrepancy cannot be ascribed to the global chemical composition or to the initial chemical abundances (such as the carbon/oxygen mass ratio) as in this case all lines would change accordingly. The two opposite trends of the H2O rotational fluxes are most likely the outcome of a difference in H2O abundance structure. To verify this hypothesis, we modelled the H2O line fluxes adopting a parametric abundance distribution rather than a full gas-grain chemistry.

thumbnail Fig. 5

Top: results of the H2O parametric models: (left) HD 163296: the models that best reproduce the observations are those with a low H2O content in the cold reservoir (x3 10−9) and with no H2O in the warm layer beyond the snow line; (right) HD 100546: models with x2~24×109$x_2^\prime \~2 - 4 \times {10^{ - 9}}$ and low values of x1 (≤10−9) and x3 (≤10−10) match the observed flux of the ground state transitions as well as the upper µmit of the high- J lines. Bottom: H2O abundance structure based on the final parametric model showing the three H2O reservoirs.

4.2.1 Water abundance parametric distribution

Starting from the dust and gas temperatures derived above, we switched off the chemical solver and the thermal balance modules. This allows us to parametrise the H2O abundance. The three H2O reservoirs can be parametrised as follows: x(H2O)={ x1for Tdust>150K,ngas>109cm3x2for Tdust<100K,Av>3.0magx3for Tgas>300K,ngas>106cm3, $x\left( {{{\rm{H}}_{\rm{2}}}{\rm{O}}} \right) = \left\{ {\matrix{ {{x_1}} \hfill &amp; {{\rm{for}}\,{T_{{\rm{dust}}}} > 150{\rm{K}},{n_{{\rm{gas}}}} > {{10}^9}{\rm{c}}{{\rm{m}}^{ - 3}}} \hfill \cr {{x_2}} \hfill &amp; {{\rm{for}}\,{T_{{\rm{dust}}}} &lt; 100{\rm{K}},Av > 3.0mag} \hfill \cr {{x_3}} \hfill &amp; {{\rm{for}}\,{T_{{\rm{gas}}}} > 300{\rm{K}},{n_{gas}} > {{10}^6}{\rm{c}}{{\rm{m}}^{ - 3}}} \hfill \cr } ,} \right.$(5)

where x1, x2, and x3 correspond to the H2O abundance in s1, s2, and s3, respectively. Such a parametrisation is tuned to match the prediction of thermo-chemical models (Figs. 3 and 4). Once the three abundance regions had been defined, we performed a grid of DALI simulations varying the values x1, x2, and x3.

In parallel with the DALI model grid, we investigated the H2O line flux ratios by performing a set of NLTE simulations with RADEX (van der Tak et al. 2007). This is useful for checking the trend of different molecular transitions at different densities and temperatures and can be used to properly set the DALI model grid.

4.2.2 H2O abundance in HD 163296

In the case of HD 163296, the non-detection of the ground-state H2O transitions points to a low H2O abundance in the cold photodesorption (s2) layer. On the other hand, the detection of high- J H2O lines suggests a medium-to-high abundance in the warm and hot layers. We run an initial grid of models varying x1 and x3 (between 10−4–10−8) and x2 (10−8–10−12). The results are reported in Appendix A. In all cases, the models do not match the observed trend. Models with high H2O abundance (≳10−5) in s1 and s3 almost match the flux of the most energetic transitions (though they still do not match the observations) but these models overestimate the upper µmits of the low- J lines independently of the value of x2 (Fig. A.1). This implies that the warm layer s3 contributes substantially to the ground-state transitions. Using RADEX we investigated how the physical conditions (gas temperature and density) affect the line flux ratio between the detected mid- J (e.g. 4i4 – 303) and undetected low- J (110 – 101) lines: the line flux ratio increases with increasing gas temperature and decreasing gas density. In particular, we find reasonable agreement with the observed flux ratio lower µmit (≳10) at ngas ≳ 108 cm−3 and Tgas ≳ 250 K; these conditions are met only in s1 and in the innermost part of the warm layer s3, at r ≲ 10 au within the snow line, right above the hot inner reservoir s1. For HD 163296, we therefore modified the parametrisation of x3 as follows: x(H2O)=x3forTgas>300K, ngas>108cm3,$x\left( {{{\rm{H}}_{\rm{2}}}{\rm{O}}} \right) = x_3^\prime {\rm{for}}{T_{{\rm{gas}}}} > 300{\rm{K, }}{n_{{\rm{gas}}}} > {10^8}{\rm{c}}{{\rm{m}}^{ - 3}},$(6)

and we performed a new grid of models varying x1, x2, and x3$x_3^\prime $. Fig. 5 (left panels) shows the grid results: models with x3~3×$x_3^\prime \~3 \times $ 10−5 and x2 < 10−10 are in good agreement with the observed trend. We note that the line fluxes are now insensitive to the value of x1, and so the actual abundance in s1 is unconstrained.

In conclusion, we find that HD 163296 hosts a H2O-rich inner disc inside the H2O snow line (r ~ 5 AU), while the disc outer region is H2O-poor. The high- J H2O lines detected by PACS require a warm inner reservoir with a H2O abundance x3$x_3^\prime $ of a few 105. On the contrary, from the non-detection of the low- J H2O lines, we constrain the H2O abundance in the photodesorption layer – defined for 3 < AV < 5 outside the disc snow line − to be x2 < 10−9 and the H2O abundance in the disc warm layer to be negligible. Based on the revised parametrisation, we computed the synthetic spectrum of H2O with DALI. The results are reported in Appendix B.

4.2.3 H2O abundance in HD 100546

In the case of HD 100546, the detection of the ground-state H2O lines implies a contribution from the cold photodesoption layer (s2). Contrary to HD 163296, the high- J transitions in the PACS range remain undetected with flux upper µmits of the order of a few 10−17 Wm−2 (Meeus et al. 2012; Fedele et al. 2013a). The deep PACS observation of the 321 – 212 transition at 75.38 µm is an order of magnitude more sensitive than the previous observations and yet the line is not detected. These non-detections suggest a low H2O abundance in the warm layer of HD 100546. We ran an initial grid of parametric models, varying the H2O abundance in the three reservoirs in the ranges 10−8−10−12 (in s1 and s3) and 10−8–10−10 (in s2). The results of the model grid are shown in Fig. A.2. As expected, the fluxes of the lowest transitions (up to the 321 − 212) are mostly dominated by the abundance in the cold reservoir, while both s1 and s3 contribute to the flux of the high- J lines.

Of key importance is the upper µmit of the 321 − 212 transition as all the models in the initial grid overestimate the upper µmit (Fig. A.2). Using RADEX, we inspected how the line flux ratios (e.g. 321 − 212/110 − 101 ≲ 8) vary with gas temperature and density, and H2O column density. The only way to match the observed flux ratio is to lower the gas temperature (≲80 K) in the line-emitting region. This imposes a very low abundance in s3 (otherwise the upper level is easily populated) as its upper layers are known to be warm from the analysis of the CO ladder. Interestingly, the independent analysis of the line velocity profile of the HIFI H2O spectra allows us to pose a robust constraint on the line-emitting region in the cold reservoir: as reported by van Dishoeck et al. (2021), the ground-state H2O emission comes from ~40 au outward. Based on this observational finding, we slightly modified the parametrisation of s2 to be µmited beyond 40 au (equivalent to Tdust < 80 K): x(H2O)=x2forTdust<80K, AV>3.0 mag.${\rm{x}}\left( {{{\rm{H}}_{\rm{2}}}{\rm{O}}} \right) = x_2^\prime {\rm{for}}{T_{{\rm{dust}}}} &lt; 80{\rm{K, }}{A_{\rm{V}}} > 3.0\,{\rm{mag}}.$(7)

With this assumption, the models with x2~3×109$x_2^\prime \~3 \times {10^{ - 9}}$ and low values of x1 (≤10−9) and x3 (≤10−10) are in good agreement with the observational trend (Fig. 5, right). Figure 6 shows the observed and synthetic line velocity profile of the ground-state transitions. The DALI synthetic spectra nicely reproduce the observed line profile. Also, in this case we computed the H2O spectrum based on the parametrised distribution and the result is shown in Appendix B.

thumbnail Fig. 6

Velocity profile of the ground-state H2O lines in HD 100546 and comparison with the DALI model spectrum. The spectra are continuum-subtracted and normalised. The DALI model nicely reproduces the observed profile of both transitions.

thumbnail Fig. 7

Top: schematic of the H2O distribution in the two discs. Middle: surface brightness profile (azimuthally averaged) of the milllmetre dust continuum at 1.2 mm for HD 163296 (from DSHARP survey, Andrews et al. 2018; Huang et al. 2018) and at 0.87 mm for HD 100546 (Pineda et al. 2019; Fedele et al. 2021) in the inner 100 au. Bottom: water-emitting region (cyan) overlaid on the ALMA dust continuum image. Notably in the case of HD 163296, the size of the H2O-emitting region corresponds to the narrow dust gap at 10 au. This may be due to dust growth at the border of the snow line.

5 Discussion and conclusion

The dichotomy of the gas-phase H2O abundance structure in HD 163296 and HD 100546 can be explained through the dynamical history of the two discs, which influences their chemistry.

Banzatti et al. (2020) proposed a link between the inner disc chemistry and the migration of solid pebbles for H2O while Booth et al. (2019) and Zhang et al. (2021) speculate that the enhanced CO column density in the disc inner 100 au of HD 163296 is due to the inward drift of icy dust particles that release the CO to the gas phase when they cross the CO frost line at r ≈ 100 au. A similar mechanism might explain the high H2O abundance within the disc snow line. Notably, a narrow dust gap is visible in the ALMA dust continuum at r = 10 au (Fig. 7). This spatial coincidence might be due to dust growth at the border of the snow line as suggested by Zhang et al. (2015) in the case of HL Tau. Further analysis is needed to verify such a hypothesis. In this regard, Notsu et al. (2019) find that the we should observe a change in the milµmetre dust opacity at the snow line position.

In HD 100546, the H2O emission region overlaps with a disc outer gap between r ~ 40 au and r ~ 150 au detected with ALMA at 870 µm (Fedele et al. 2021). Notably, this is spatially coincident with the 3 µm H2O ice absorption band detected by Honda et al. (2016) that extends from r ~ 40 au to r ~ 120 au. The high H2O abundance in the photodesorption layer can be explained through the spatial coincidence of the dust gap and the H2O-emitting region: the presence of the outer disc gap lets the UV stellar photons penetrate to the colder layers to make H2O ice subµmate through photodesorption. This is shown in the schematics presented in Fig. 7. This scenario can also explain the origin of the crystalline forsterite emission at 69 µm: the blackbody fitting suggests that the emission arise from dust grains at ~70 K at radial distances r > 50 au (Sturm et al. 2010). As for the H2O ice, the crystalline forsterite is directly illuminated by the stellar UV photons that enter in the dust gap.

Interestingly, Booth et al. (2021) report the detection of gas-phase CH3OH towards HD 100546: contrary to H2O, the CH3OH emission peaks in the inner 60 au of the disc and extends out to nearly 300 au. According to Booth et al. (2021), the CH3OH emission is due to desorption of the primordial ice inherited from the cold dark cloud: in particular, the bulk of the emission from the inner disc is due to thermal desorption mostly from the inner cavity at 13 au, while the spatially extended emission is due to photodesorption, similarly to the scenario proposed here for H2O.

It is worth noting that the only other disc with a detection of cold gas-phase H2O, TW Hya (Hogerheijde et al. 2011), also shows CH3OH emission (Walsh et al. 2016).

The low H2O abundance in the outer warm molecular layer in both HD 163296 and HD 100546 implies an O-poor molecular layer compared to the standard thermochemical predictions. This is in agreement with the findings of low gas-phase CO abundances and large carbon/oxygen mass ratios observed in several discs (e.g. Bergin et al. 2010; Favre et al. 2013; McClure et al. 2016; Kama et al. 2016; Miotello et al. 2017; Semenov et al. 2018).

In both cases, oxygen in the outer disc is incorporated into dust particles settled to the midplane (e.g. Krijt et al. 2020). In the case of HD 163296, the radial migration of icy pebbles to the inner disc decreases the H2O content in the outer disc and enhances the gas-phase reservoir inside the H2O snow line. On the other hand, in HD 100546, the presence of the outer wide gap (see Fig. 7 bottom panels) prevents or slows down the radial drift of the ice dust grains towards the high-temperature region, resulting in a (relatively) high gas-phase H2O abundance in the outer disc. Compared to the wide and deep gap of HD 100546, the three outer gaps observed in HD 163296 (Fig. 7) are not wide and deep enough to stop the inward drift of icy grains. Our results show how the chemistry of discs is strictly connected to the dynamical history of the system and to the presence of giant protoplanets.

Acknowledgements

D.F. acknowledges the support of the Italian National Institute of Astrophysics (INAF) through the INAF Mainstream projects ARIEL and the “Astrochemical Link between Circumstellar Disks and Planets”, “Protoplanetary Disks Seen through the Eyes of New- generation Instruments” and by the PRIN-INAF 2019 Planetary Systems At Early Ages (PLATEA). This project has received funding from the European Union’s Horizon 2020 research and innovation programme under the Marie Sklodowska-Curie grant agreement No 823823 (DUSTBUSTERS). This work has made use of data from the European Space Agency (ESA) mission Gaia (https://www.cosmos.esa.int/gaia), processed by the Gaia Data Processing and Analysis Consortium (DPAC, https://www.cosmos.esa.int/web/gaia/dpac/consortium). Funding for the DPAC has been provided by national institutions, in particular the institutions participating in the Gaia Multilateral Agreement. HIFI was designed and built by a consortium of institutes and university departments from across Europe, Canada and the US under the leadership of SRON Netherlands Institute for Space Research, Groningen, The Netherlands with major contributions from Germany, France and the US. Consortium members are: Canada: CSA, U. Waterloo; France: CESR, LAB, LERMA, IRAM; Germany: KOSMA, MPIfR, MPS; Ireland, NUI Maynooth; Italy: ASI, IFSI-INAF, Arcetri-INAF; Netherlands: SRON, TUD; Poland: CAMK, CBK; Spain: Observatorio Astronömico Nacional (IGN), Centro de Astrobiologia (CSICINTA); Sweden: Chalmers University of Technology - MC2, RSS & GARD, Onsala Space Observatory, Swedish National Space Board, Stock-holm University - Stockholm Observatory; Switzerland: ETH Zürich, FHNW; USA: Caltech, JPL, NHSC. PACS has been developed by a consortium of institutes led by MPE (Germany) and including UVIE (Austria); KUL, CSL, IMEC (Belgium); CEA, OAMP (France); MPIA (Germany); IFSI, OAP/OAT, OAA/CAISMI, LENS, SISSA (Italy); IAC (Spain). This development has been supported by the funding agencies BMVIT (Austria), ESA-PRODEX (Bel-gium), CEA/CNES (France), DLR (Germany), ASI (Italy), and CICYT/MCYT(Spain).

Appendix A H2O model grid

The results of the parametric model grid for HD 163296 and HD 100546 are shown in Figures A.1 and A.2. The models shown here refer to the original definition of the H2O reservoir (Equation 5). Only a portion of the grid is shown here.

thumbnail Fig. A.1

Results of the parametric model grid for HD 163296.

In the case of HD 163296, none of the models are able to reproduce the observed trend: on the one hand, the high-J line requires an high abundance (≳ 10−5, Figure A.1) in the warm layer (s3) but these models overestimate the upper µmit of the two ground-state transitions. Lowering the value of x2 does not help.

Figure A.2 shows a subset of the model grid for two different values of x2, 10−9 and 4 × 10−9. In the first case, the models underestimate the flux of the ground-state transitions. Increasing the H2O abundance in s2 is better suited to fit the flux of the cold transitions but note that these models fail to reproduce the upper limit to the 321 − 212 transition.

thumbnail Fig. A.2

Same as Figure A.1 but for HD 100546.

Appendix B Synthetic H2O spectrum

Figure B.1 shows the synthetic H2O spectra for both discs based on the DALI parametric distribution reported in Figure 5. The spectroscopic data for the ray tracing are from the LAMDA database (Schöier et al. 2005, and reference therein). Submilµmetre observations with e.g. ALMA are potentially able to detect emission of H2O isotopologues, such as H218O${\rm{H}}_{\rm{2}}^{{\rm{18}}}{\rm{O}}$. According to our predictions, the brightest H2O submilµmetre transition is the 414 − 321 (lower energy level of 212 K); while the main isotopologue cannot be observed from the ground because of atmospheric absorption, the corresponding H218O${\rm{H}}_{\rm{2}}^{{\rm{18}}}{\rm{O}}$ transition at 390,607 GHz falls in the ALMA band 8. Assuming a H2O / H218O${\rm{H}}_{\rm{2}}^{{\rm{18}}}{\rm{O}}$ abundance ratio of 560, the predicted flux of the H218O${\rm{H}}_{\rm{2}}^{{\rm{18}}}{\rm{O}}$ 414 − 321 transition is 7.5 × 10−24 Wm−2 and 1 × 10−24 Wm−2 for HD 163296 and HD 100546, respectively. These lines are too faint to be detected with current instrumentation in a reasonable amount of time.

thumbnail Fig. B.1

Synthetic H2O spectrum of HD 163296 and HD 100546 based on the H2O parametric distribution reported in Figure 5.

References

  1. Ádámkovics, M., Glassgold, A. E., & Najita, J. R. 2014, ApJ, 786, 135 [CrossRef] [Google Scholar]
  2. Andrews, S. M., Wilner, D. J., Espaillat, C., et al. 2011, ApJ, 732, 42 [Google Scholar]
  3. Andrews, S. M., Huang, J., Pérez, L. M., et al. 2018, ApJ, 869, L41 [NASA ADS] [CrossRef] [Google Scholar]
  4. Banzatti, A., Pontoppidan, K. M., Salyk, C., et al. 2017, ApJ, 834, 152 [NASA ADS] [CrossRef] [Google Scholar]
  5. Banzatti, A., Pascucci, I., Bosman, A. D., et al. 2020, ApJ, 903, 124 [NASA ADS] [CrossRef] [Google Scholar]
  6. Bergin, E. A., Hogerheijde, M. R., Brinch, C., et al. 2010, A&A, 521, L33 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  7. Bethell, T., & Bergin, E. 2009, Science, 326, 1675 [NASA ADS] [CrossRef] [Google Scholar]
  8. Booth, A. S., Walsh, C., Ilee, J. D., et al. 2019, ApJ, 882, L31 [Google Scholar]
  9. Booth, A. S., Walsh, C., Terwisscha van Scheltinga, J., et al. 2021, Nat. Astron., 5, 684 [NASA ADS] [CrossRef] [Google Scholar]
  10. Bruderer, S. 2013, A&A, 559, A46 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  11. Bruderer, S., van Dishoeck, E. F., Doty, S. D., & Herczeg, G. J. 2012, A&A, 541, A91 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  12. Carr, J. S., & Najita, J. R. 2008, Science, 319, 1504 [Google Scholar]
  13. D’Alessio, P., Calvet, N., Hartmann, L., Franco-Hernández, R., & Servín, H. 2006, ApJ, 638, 314 [Google Scholar]
  14. de Graauw, T., Helmich, F. P., Phillips, T. G., et al. 2010, A&A, 518, L6 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  15. Dent, W. R. F., Thi, W. F., Kamp, I., et al. 2013, PASP, 125, 477 [NASA ADS] [CrossRef] [Google Scholar]
  16. Doppmann, G. W., Najita, J. R., Carr, J. S., & Graham, J. R. 2011, ApJ, 738, 112 [NASA ADS] [CrossRef] [Google Scholar]
  17. Du, F., & Bergin, E. A. 2014, ApJ, 792, 2 [NASA ADS] [CrossRef] [Google Scholar]
  18. Du, F., Bergin, E. A., Hogerheijde, M., et al. 2017, ApJ, 842, 98 [Google Scholar]
  19. Favre, C., Cleeves, L. I., Bergin, E. A., Qi, C., & Blake, G. A. 2013, ApJ, 776, L38 [Google Scholar]
  20. Fedele, D., Pascucci, I., Brittain, S., et al. 2011, ApJ, 732, 106 [NASA ADS] [CrossRef] [Google Scholar]
  21. Fedele, D., Bruderer, S., van Dishoeck, E. F., et al. 2012, A&A, 544, L9 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  22. Fedele, D., Bruderer, S., van Dishoeck, E. F., et al. 2013a, A&A, 559, A77 [CrossRef] [EDP Sciences] [Google Scholar]
  23. Fedele, D., Bruderer, S., van Dishoeck, E. F., et al. 2013b, ApJ, 776, L3 [NASA ADS] [CrossRef] [Google Scholar]
  24. Fedele, D., van Dishoeck, E. F., Kama, M., Bruderer, S., & Hogerheijde, M. R. 2016, A&A, 591, A95 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  25. Fedele, D., Toci, C., Maud, L., & Lodato, G. 2021, A&A, 651, A90 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  26. Gaia Collaboration (Brown, A. G. A., et al.) 2021, A&A, 649, A1 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  27. Glassgold, A. E., Meijerink, R., & Najita, J. R. 2009, ApJ, 701, 142 [NASA ADS] [CrossRef] [Google Scholar]
  28. Green, J. D., Evans, N. J., II, Jørgensen, J. K., et al. 2013, ApJ, 770, 123 [Google Scholar]
  29. Hogerheijde, M. R., Bergin, E. A., Brinch, C., et al. 2011, Science, 334, 338 [Google Scholar]
  30. Honda, M., Kudo, T., Takatsuki, S., et al. 2016, ApJ, 821, 2 [NASA ADS] [CrossRef] [Google Scholar]
  31. Huang, J., Andrews, S. M., Dullemond, C. P., et al. 2018, ApJ, 869, L42 [NASA ADS] [CrossRef] [Google Scholar]
  32. Isella, A., Guidi, G., Testi, L., et al. 2016, Phys. Rev. Lett., 117, 251101 [NASA ADS] [CrossRef] [Google Scholar]
  33. Kama, M., Bruderer, S., van Dishoeck, E. F., et al. 2016, A&A, 592, A83 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  34. Kama, M., Trapman, L., Fedele, D., et al. 2020, A&A, 634, A88 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  35. Krijt, S., Bosman, A. D., Zhang, K., et al. 2020, ApJ, 899, 134 [NASA ADS] [CrossRef] [Google Scholar]
  36. Mandell, A. M., Mumma, M. J., Blake, G. A., et al. 2008, ApJ, 681, L25 [Google Scholar]
  37. McClure, M. K., Bergin, E. A., Cleeves, L. I., et al. 2016, ApJ, 831, 167 [Google Scholar]
  38. Meeus, G., Montesinos, B., Mendigutía, I., et al. 2012, A&A, 544, A78 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  39. Meeus, G., Salyk, C., Bruderer, S., et al. 2013, A&A, 559, A84 [CrossRef] [EDP Sciences] [Google Scholar]
  40. Miotello, A., van Dishoeck, E. F., Williams, J. P., et al. 2017, A&A, 599, A113 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  41. Najita, J. R., Ádámkovics, M., & Glassgold, A. E. 2011, ApJ, 743, 147 [NASA ADS] [CrossRef] [Google Scholar]
  42. Notsu, S., Akiyama, E., Booth, A., et al. 2019, ApJ, 875, 96 [Google Scholar]
  43. Pineda, J. E., Szulágyi, J., Quanz, S. P., et al. 2019, ApJ, 871, 48 [Google Scholar]
  44. Podio, L., Kamp, I., Codella, C., et al. 2013, ApJ, 766, L5 [NASA ADS] [CrossRef] [Google Scholar]
  45. Poglitsch, A., Waelkens, C., Geis, N., et al. 2010, A&A, 518, L2 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  46. Pontoppidan, K. M., Salyk, C., Blake, G. A., & Käufl, H. U. 2010, ApJ, 722, L173 [NASA ADS] [CrossRef] [Google Scholar]
  47. Riviere-Marichalar, P., Ménard, F., Thi, W. F., et al. 2011, A&A, 538, L3 [Google Scholar]
  48. Roelfsema, P. R., Helmich, F. P., Teyssier, D., et al. 2012, A&A, 537, A17 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  49. Salyk, C., Pontoppidan, K. M., Blake, G. A., et al. 2008, ApJ, 676, L49 [NASA ADS] [CrossRef] [Google Scholar]
  50. Schöier, F. L., van der Tak, F. F. S., van Dishoeck, E. F., & Black, J. H. 2005, A&A, 432, 369 [Google Scholar]
  51. Semenov, D., Favre, C., Fedele, D., et al. 2018, A&A, 617, A28 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  52. Sturm, B., Bouwman, J., Henning, T., et al. 2010, A&A, 518, L129 [CrossRef] [EDP Sciences] [Google Scholar]
  53. Thi, W.-F., van Zadelhoff, G.-J., & van Dishoeck, E. F. 2004, A&A, 425, 955 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  54. van der Tak, F. F. S., Black, J. H., Schöier, F. L., Jansen, D. J., & van Dishoeck, E. F. 2007, A&A, 468, 627 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  55. van der Wiel, M. H. D., Naylor, D. A., Kamp, I., et al. 2014, MNRAS, 444, 3911 [NASA ADS] [CrossRef] [Google Scholar]
  56. van Dishoeck, E. F., Kristensen, L. E., Mottram, J. C., et al. 2021, A&A, 648, A24 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  57. Vioque, M., Oudmaijer, R. D., Baines, D., Mendigutía, I., & Pérez-Martínez, R. 2018, A&A, 620, A128 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  58. Walsh, C., Juhász, A., Pinilla, P., et al. 2014, ApJ, 791, L6 [Google Scholar]
  59. Walsh, C., Juhász, A., Meeus, G., et al. 2016, ApJ, 831, 200 [Google Scholar]
  60. Woitke, P., Kamp, I., & Thi, W. 2009, A&A, 501, 383 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  61. Woodall, J., Agúndez, M., Markwick-Kemper, A. J., & Millar, T. J. 2007, A&A, 466, 1197 [CrossRef] [EDP Sciences] [Google Scholar]
  62. Yang, B., Stancil, P. C., Balakrishnan, N., & Forrey, R. C. 2010, ApJ, 718, 1062 [Google Scholar]
  63. Zhang, K., Pontoppidan, K. M., Salyk, C., & Blake, G. A. 2013, ApJ, 766, 82 [NASA ADS] [CrossRef] [Google Scholar]
  64. Zhang, K., Blake, G. A., & Bergin, E. A. 2015, ApJ, 806, L7 [Google Scholar]
  65. Zhang, K., Booth, A. S., Law, C. J., et al. 2021, ApJS, 257, 5 [NASA ADS] [CrossRef] [Google Scholar]

All Tables

Table 1

Observations log.

Table 2

H2O line fluxes.

Table 3

DALI model parameters.

All Figures

thumbnail Fig. 1

Herschel/PACS deep integration of HD 100546 at 75 µm. The CO J = 35–34 line is detected along with two other emission lines likely due to NH3. The blue dashed line indicates the position of the H2O 321 – 212 transition.

In the text
thumbnail Fig. 2

Herschel/HIFI spectrum of HD 163296110 − 101 (at 556 GHz) at rest-frame velocity.

In the text
thumbnail Fig. 3

DALI disc structure for HD 163296. Top panels: gas and dust density structure, the gas temperature, and the H2O abundance structure. The three main H2O reservoirs are indicated by the labels s1, s2, and s3. Bottom panels: model predictions of the SED (left), CO ladder (middle), and H2O fluxes (right). The observed fluxes and 3σ upper µmits are shown as (grey) circles and open triangles, respectively. The representative model shown as a red line and red squares) reproduces the SED (left) and the CO rotational ladder (middle) well but is not able to reproduce the H2O line fluxes and upper µmits. We note that the H2O abundances shown here are those of the representative disc model based on the full gas-grain chemistry calculation.

In the text
thumbnail Fig. 4

Same as Fig. 3 but for HD 100546.

In the text
thumbnail Fig. 5

Top: results of the H2O parametric models: (left) HD 163296: the models that best reproduce the observations are those with a low H2O content in the cold reservoir (x3 10−9) and with no H2O in the warm layer beyond the snow line; (right) HD 100546: models with x2~24×109$x_2^\prime \~2 - 4 \times {10^{ - 9}}$ and low values of x1 (≤10−9) and x3 (≤10−10) match the observed flux of the ground state transitions as well as the upper µmit of the high- J lines. Bottom: H2O abundance structure based on the final parametric model showing the three H2O reservoirs.

In the text
thumbnail Fig. 6

Velocity profile of the ground-state H2O lines in HD 100546 and comparison with the DALI model spectrum. The spectra are continuum-subtracted and normalised. The DALI model nicely reproduces the observed profile of both transitions.

In the text
thumbnail Fig. 7

Top: schematic of the H2O distribution in the two discs. Middle: surface brightness profile (azimuthally averaged) of the milllmetre dust continuum at 1.2 mm for HD 163296 (from DSHARP survey, Andrews et al. 2018; Huang et al. 2018) and at 0.87 mm for HD 100546 (Pineda et al. 2019; Fedele et al. 2021) in the inner 100 au. Bottom: water-emitting region (cyan) overlaid on the ALMA dust continuum image. Notably in the case of HD 163296, the size of the H2O-emitting region corresponds to the narrow dust gap at 10 au. This may be due to dust growth at the border of the snow line.

In the text
thumbnail Fig. A.1

Results of the parametric model grid for HD 163296.

In the text
thumbnail Fig. A.2

Same as Figure A.1 but for HD 100546.

In the text
thumbnail Fig. B.1

Synthetic H2O spectrum of HD 163296 and HD 100546 based on the H2O parametric distribution reported in Figure 5.

In the text

Current usage metrics show cumulative count of Article Views (full-text article views including HTML views, PDF and ePub downloads, according to the available data) and Abstracts Views on Vision4Press platform.

Data correspond to usage on the plateform after 2015. The current usage metrics is available 48-96 hours after online publication and is updated daily on week days.

Initial download of the metrics may take a while.