Chemically tracing the water snowline in protoplanetary disks with HCO +
M. Leemker, M. L. R. van 't Hoff, L. Trapman, M. L. van Gelder, M. R. Hogerheijde, D. Ruíz-Rodríguez, E. F. van Dishoeck
AAstronomy & Astrophysics manuscript no. Leemker_et_al_accepted © ESO 2020November 26, 2020
Chemically tracing the water snowline in protoplanetary disks withHCO + M. Leemker , M. L. R. van ’t Ho ff , , L. Trapman , M. L. van Gelder , M. R. Hogerheijde , , D. Ruíz-Rodríguez , andE. F. van Dishoeck , Leiden Observatory, Leiden University, P.O. box 9513, 2300 RA Leiden, The Netherlandse-mail: [email protected] University of Michigan, Department of Astronomy, 1085 S. University, Ann Arbor, MI 48109, USA Anton Pannekoek Institute for Astronomy, University of Amsterdam, Science Park 904, 1090 GE Amsterdam, The Netherlands National Radio Astronomy Observatory, 520 Edgemont Road, Charlottesville, VA 22903-2475, USA Max-Planck-Institut für Extraterrestrische Physik, Giessenbachstrasse 1, 85748 Garching, GermanyNovember 26, 2020
ABSTRACT
Context.
The formation of planets is expected to be enhanced around snowlines in protoplanetary disks, in particular around the watersnowline. Moreover, freeze-out of abundant volatile species in disks alters the chemical composition of the planet-forming material.However, the close proximity of the water snowline to the host star combined with the di ffi culty of observing water from Earth makesa direct detection of the water snowline in protoplanetary disks challenging. HCO + is a promising alternative tracer of the watersnowline. The destruction of HCO + is dominated by gas-phase water, leading to an enhancement in the HCO + abundance once wateris frozen out. Aims.
Following earlier observed correlations between water and H CO + emission in a protostellar envelope, the aim of this researchis to investigate the validity of HCO + and the optically thin isotopologue, H CO + , as tracers of the water snowline in protoplanetarydisks and the required sensitivity and resolution to observationally confirm this. Methods.
A typical Herbig Ae disk structure is assumed and its temperature structure is modelled with the thermochemical code
DALI . Two small chemical networks are then used and compared to predict the HCO + abundance in the disk; one without water andone including water. Subsequently, the corresponding emission profiles are modelled for the J = − CO + andHCO + , which provides the best balance between brightness and optical depth e ff ects of the continuum emission, and is less a ff ectedby blending with complex molecules. Models are then compared with archival ALMA data. Results.
The HCO + abundance jumps by two orders of magnitude over a radial range of 2 AU outside the water snowline, whichin our model is located at 4.5 AU. We find that the emission of H CO + and HCO + is ring-shaped due to three e ff ects: destructionof HCO + by gas-phase water, continuum optical depth, and molecular excitation e ff ects. Comparing the radial emission profiles for J = − (cid:48)(cid:48) .
05 beam reveals that the presence of gas-phase water causes an additional drop of only ∼
13% and 24%in the center of the disk, for H CO + and HCO + , respectively. For the much more luminous outbursting source V883 Ori, our modelspredict that the e ff ects of dust and molecular excitation are not limiting HCO + as a snowline tracer if the snowline is located at radiilarger than ∼
40 AU. Our analysis of recent archival ALMA band 6 observations of the J = − + is consistentwith the water snowline located around 100 AU, further out than was previously estimated from an intensity break in the continuumemission. Conclusions.
The HCO + abundance drops steeply around the water snowline, when water desorbs in the inner disk, but continuumoptical depth and molecular excitation e ff ects conceal the drop in HCO + emission due to the water snowline. Therefore locating thewater snowline with HCO + observations in disks around Herbig Ae stars is very di ffi cult, but it is possible for disks around outburstingstars such as V883 Ori, where the snowline has moved outwards. Key words.
Astrochemistry; protoplanetary disks; ISM: molecules; submillimeter: planetary systems
1. Introduction
High resolution observations of protoplanetary disks show thatmany of them have rings, gaps, and other substructures (e.g., An-drews et al. 2010; van der Marel et al. 2015; Fedele et al. 2017;Huang et al. 2018; Long et al. 2018; and Andrews 2020 for a re-view). Di ff erent explanations have been proposed for these sub-structures such as planets (e.g. Bryden et al. 1999; Zhu et al.2014; Dong et al. 2018) and snowlines (e.g. Banzatti et al. 2015;Zhang et al. 2015; Okuzumi et al. 2016). A snowline is the mid-plane radius in a protoplanetary disk where 50% of one of themajor volatiles is in the gas phase and 50% is frozen out ontothe dust grains. These snowlines are related to planet formation, because dust properties change around the H O and CO or CO snowlines (Pinilla et al. 2017). Even though water ice mantlesmay not always aid dust coagulation by collisions (Kimura et al.2020), the sublimation, condensation, and di ff usion of gas-phasewater enhances the surface density around the water snowlineaiding planet formation and potentially triggering the streaminginstability (e.g. Stevenson & Lunine 1988; Dr ˛a˙zkowska & Alib-ert 2017; Schoonenberg & Ormel 2017).Snowlines do not only a ff ect the formation of planets but theyalso a ff ect their chemical composition. The sequential freeze-outof the major volatile species changes the bulk chemical compo-sition, often measured as the C / O ratio, of the planet-forming
Article number, page 1 of 19 a r X i v : . [ a s t r o - ph . E P ] N ov & A proofs: manuscript no. Leemker_et_al_accepted material with the ice becoming more oxygen-rich than the gas(Öberg et al. 2011; Eistrup et al. 2016, 2018). The water snow-line is of particular importance as water is crucial for the originof life.Knowledge of the water snowline location is thus essentialfor the understanding of planet formation and composition. Yetobserving water snowlines is challenging because their expectedlocation is only a few AU from the host star for T Tauri disks(Harsono et al. 2015) and ∼
10 AU for disks around more lu-minous Herbig Ae stars. Therefore, observations with very highspatial resolution are necessary. On top of that, observing wa-ter snowlines directly from the ground by observing gas-phasewater in protoplanetary disks is challenging because of the ab-sorption of water in the Earth’s atmosphere. The only detectionsof water in protoplanetary disks probe the inner most hot parts( (cid:46)
O, reduces the atmo-spheric absorption, but detecting these molecules at a signifi-cant signal-to-noise ratio greatly increases the required observ-ing time. Space-based telescopes circumvent this problem, butso far, lacked the spatial resolution to resolve the water snow-line.Chemical tracers provide an alternative way of locatingsnowlines. Some molecules show ring-shaped emission evenif the density of the disk is smooth. One example of such amolecule is CN, which is enhanced if a strong UV field is present(e.g. van Zadelho ff et al. 2003; Teague et al. 2016; Cazzolettiet al. 2018) and other examples are molecules that trace snow-lines. Chemical imaging has been used to locate the CO snowlinein the disks around TW Hya and HD 163296 using N H + andDCO + (Mathews et al. 2013; Qi et al. 2013, 2015, 2019). Bothtracers show ring-shaped emission but detailed chemical mod-elling is needed to infer the location of the CO snowline fromthese observations (Aikawa et al. 2015; van ’t Ho ff et al. 2017;Carney et al. 2018).Water can be traced chemically using the HCO + ion, whichis destroyed by gas-phase water (Phillips et al. 1992; Bergin et al.1998):HCO + + H O → CO + H O + . (1)Based on this reaction, HCO + is expected to be abundant whenwater is frozen out onto the grains and HCO + is depleted whenwater is in the gas phase. This anti-correlation between waterand HCO + has been verified observationally by the optically thinisotopologue H CO + and an isotopologue of water, H O, in aprotostellar envelope (van ’t Ho ff et al. 2018a). Due to the highluminosity of young stellar objects at this early stage, the snow-line is located further away from the star than in protoplanetarydisks (Visser & Bergin 2012; Jørgensen et al. 2013; Vorobyovet al. 2013; Visser et al. 2015; Cieza et al. 2016; Hsieh et al.2019). Similarly, the higher luminosity of Herbig Ae stars com-pared to T Tauri stars increases the temperature in the proto-planetary disks around them, which locates the water snowlinefurther out in disks around the former type. Therefore, we focuson protoplanetary disks around Herbig Ae stars in this work.The aim of this paper is to investigate HCO + as a tracer of thewater snowline in protoplanetary disks. Compared to protostel-lar envelopes the 2D temperature and density structure of diskscomplicates matters. First, the location of the water snowline in protoplanetary disks around Herbig Ae stars is expected around10 AU, in contrast to (cid:38)
100 AU in protostellar envelopes. Second,the higher column density in disks increases the optical depth ofthe continuum and line emission. This complicates locating thewater snowline as emission from HCO + can be absorbed by dustparticles, mimicking the e ff ect of depletion due to gas-phase wa-ter on the HCO + emission. Finally, the HCO + emission can alsobe ring-shaped when the J -level population of the low levels,that can be observed with ALMA, decreases towards the centerof the star due to the increase in temperature and density in theinner disk. This also mimics the e ff ect of gas-phase water on theHCO + emission.The validity of HCO + as a tracer of the water snowline inprotoplanetary disks is investigated by modelling the physicaland chemical structure of a disk around a typical Herbig Ae starwith a luminosity of 36 L (cid:12) as described in Section 2. The pre-dicted abundance structure of HCO + is presented in Section 3.1,and the corresponding emission profiles of HCO + and H CO + are described in Section 3.2 and 3.3. Archival ALMA observa-tions of H CO + in the outbursting source V883 Ori are dis-cussed in Section 4. Finally we conclude in Section 5 that it isdi ffi cult to use HCO + as a tracer of the water snowline in disksaround Herbig Ae stars, but that it is possible for outburstingsources such as V883 Ori.
2. Protoplanetary disk model
The HCO + and H CO + radial emission profiles are modelled inseveral steps. First, a density structure for a typical disk arounda Herbig Ae star is assumed. Second, the gas temperature is cal-culated using the thermochemical code DALI (Bruderer et al.2009, 2012; Bruderer 2013). The gas density and gas tempera-ture are then used as input for two chemical models, one withoutwater and one with water, that predict the HCO + abundance inthe disk. These abundance profiles together with the density andtemperature structure of the disk are used to model the HCO + and H CO + emission profiles with DALI . The density structure of the disk is modelled with the thermo-chemical code
DALI following the approach of Andrews et al.(2011). This approach is based on the self-similar solution for aviscously evolving disk, where the gas surface density of the diskoutside the dust sublimation radius follows a power law with anexponential taper (Lynden-Bell & Pringle 1974; Hartmann et al.1998): Σ gas ( R ) = Σ c (cid:32) RR c (cid:33) − γ exp − (cid:32) RR c (cid:33) − γ , (2)with Σ c the gas surface density at the characteristic radius R c and γ the power law index. A sublimation radius of R subl = . × cm − . An overview of all parameters used for the densityand temperature structure can be found in Table 1.The gas follows a Gaussian distribution in the vertical direc-tion, where the temperature is calculated explicitly at each loca-tion in the disk. The scale height of the gas is set by the flaringindex ψ and the characteristic scaleheight h c at R c , h = h c (cid:32) RR c (cid:33) ψ . (3) Article number, page 2 of 19. Leemker et al.: Chemically tracing the water snowline in protoplanetary disks with HCO + Table 1:
DALI model parameters for the disk around a typicalHerbig Ae star.Model parameter Value
Physical structureR subl R c
50 AU Σ c − M disk .
01 M (cid:12) γ h c ψ Dust properties χ f ls ∆ gas / dust Stellar spectrum (1)
Type Herbig L (cid:63)
36 L (cid:12) L X − T e ff T X ζ c . r . − Stellar propertiesM (cid:63) . (cid:12) Observational geometryi ◦ d
100 pc
Notes. a ( b ) represents a × b . (1) Spectrum of HD 100546 (Kama et al.2016), which is well approximated by a 10 K black body spectrum.
The resulting gas density of our model is shown in the top panelof Fig. 1.The dust surface density is modelled with
DALI by scalingthe gas surface density with the disk-averaged gas-to-dust massratio, ∆ gas / dust , which is set to the ISM value of 100. The verticalstructure of the dust is modelled using two populations of dustgrains following the approach of D’Alessio et al. (2006). Smallgrains with sizes between 5 nm and 1 µ m are well mixed with thegas and therefore follow the same vertical structure of the gas.Large grains with sizes from 5 nm to 1 mm are settled to the diskmidplane. This is modelled by reducing the scale height of thelarge grains with a factor χ <
1. The fraction of the dust mass inlarge grains is controlled by f ls and the size distribution of bothgrain populations is assumed to be proportional to a − . with a the size of the grains (MRN distribution, Mathis et al. 1977).The gas temperature needs to be computed separately as theabundances of the molecules that act as coolants of the gas in thedisk a ff ect the gas temperature, which in turn a ff ects the chem-istry. Therefore, DALI solves for the gas temperature by iterat-ing over the chemistry and heating by the photoelectric e ff ectand the cooling by molecules in the gas. The dust temperature iscomputed using Monte Carlo continuum radiative transfer. Thegas and dust temperature calculated by DALI are shown in themiddle and bottom panel of Fig. 1. The gas temperature exceedsthe dust temperature in the surface layers only, whereas the gasand dust temperatures are equal in deeper layers of the disk. Fig. 1: Assumed gas density (top) and computed gas and dusttemperature (middle and bottom panels) in the
DALI model for atypical disk around a Herbig Ae star with a luminosity of 36 L (cid:12) .Only the region where n (gas) > cm − is shown. The positionof the star is indicated with the symbol of a black star and the T dust =
150 K line, approximating the water snow surface, isindicated with the grey solid line.
The HCO + abundance with respect to the total number of hydro-gen atoms in the disk is modelled with a small chemical network,which is solved time-dependently up to 1 Myr with the python Article number, page 3 of 19 & A proofs: manuscript no. Leemker_et_al_accepted
Table 2: Reaction coe ffi cients of the reactions used in this paper. The last two columns indicate in which chemical networks thereactions are used.reaction α β γ k (100 K) k (150 K) units ref. NW W ζ c . r . H + c . r . → H + + e − −
17) 0 0 1( −
17) 1( −
17) s − k H + + H → H + + H 2 . −
9) 0 0 2 . −
9) 2 . −
9) cm s − (a) k H + + e − → H + H 2 . − − . . −
8) 3 . −
8) cm s − (b) k CO + H + → HCO + + H . − − . − . . −
9) 1 . −
9) cm s − (c) k e − HCO + + e − → CO + H 2 . − − . . −
7) 3 . −
7) cm s − (d) k H O HCO + + H O → CO + H O + . − − . . −
9) 3 . −
9) cm s − (e) k H O + + e − → H O + H 7 . − − . . −
7) 1 . −
7) cm s − (f) k H O + H + → H O + + H . − − . . −
8) 8 . −
9) cm s − (g, h) k f H O → H O (ice) 1 . − n (H )10(12)cm − . . − n (H )10(12)cm − . − n (H )10(12)cm − s − (i) k d H O (ice) → H O 2 . . . −
13) 5 . −
5) s − (i) Notes. a ( b ) represents a × b . NW is the chemical network without water, whereas W includes water. All reactions are given in the form k = α ( T gas / β exp( − γ/ T gas ), where the values of α , β , and γ for k , k , k , k e − , k H O , k and k are taken from the rate12 UMIST database(McElroy et al. 2013) and k d uses T dust instead of T gas . References: (a) Theard & Huntress (1974), (b) McCall et al. (2004), (c) Klippenstein et al.(2010), (d) Mitchell (1990), (e) Adams et al. (1978), (f) Novotný et al. (2010), (g) Kim et al. (1974), (h) Anicich et al. (1975) and (i) Appendix Aand references therein. H H + H + H OH O (ice) H O + HCO + CO W NW ζ c . r . e − H e − H + e − H + H O Fig. 2: A schematic view of chemical networks NW (no water)and W (water) used to predict the abundance of HCO + . The firstnetwork, network NW, only includes the reactions enclosed inthe dashed box. The second network, network W, includes allreactions present in the figure. The destruction of HCO + by gas-phase water is indicated with the thick red arrow.function odeint . This allows us to more easily study the ef-fect of di ff erent parameters than with a computationally moreexpensive full chemical model. The HCO + abundance in proto-planetary disks is mainly controlled by three reactions. The mainformation route of HCO + involves gas-phase CO:CO + H + → HCO + + H . (4)Here, H + is produced by cosmic ray ionization of molecular hy-drogen at a rate of ζ c . r . . On the other hand, HCO + is destroyed bygas-phase water (Eq. 1) and by dissociative recombination with The odeint function is part of the Scipy package in python and usesthe LSODA routine from the ODEPACK library in FORTRAN. an electron:HCO + + e − → CO + H . (5)Therefore, a jump in the HCO + abundance around the watersnowline is expected if electrons are not the dominant destruc-tion mechanism of HCO + .To investigate the relationship between HCO + and H O, twosmall chemical networks are used. The first network, network nowater (NW), is shown in the dashed box in Fig. 2. This networkincludes three reactions to model the ionization in the disk (re-actions ζ c . r . , k and k in Table 2), and two reactions to model theHCO + abundance in the absence of gas-phase water (reactions k and k e − (Eq. 4 and 5) in Table 2). As reaction k H O (Eq. 1) isnot included in this network, HCO + is not expected to trace thewater snowline in network NW. Therefore, this network servesas a baseline to quantify the e ff ect of the gas-phase water on theHCO + abundance. The second network, network water (W), isshown in the solid box in Fig. 2 and contains all reactions presentin Table 2. Network W includes those in network NW togetherwith reactions to model the e ff ect of water on the HCO + abun-dance. The most important reaction for our study is indicated inred in Fig. 2 and is the destruction of HCO + by gas-phase water(reaction k H O or Eq. 1). The other reactions present in networkW include the freeze-out and desorption of water (reactions k f and k d ) and the formation and destruction of gas-phase waterfrom and to H O + (reactions k and k ). Therefore, HCO + is ex-pected to trace the water snowline in network W. The equationsfor the freeze-out and desorption rates of water can be found inAppendix A.Initially, the abundance of all gas- and ice-phase species areset to 0 except for the abundance of gas-phase CO, which isset to 10 − both in network NW and W, and the abundance ofgas-phase water in network W, which is set to 3 . × − , ap-propriate for a dark cold cloud (McElroy et al. 2013), see alsoTable B.1. The rate coe ffi cient for the freeze-out of water de-pends on multiple parameters, including the number density andsize of the grains. Therefore, we assume in chemical networkNW and W a typical grain number density of 10 − × n (H ) anda grain size of 0.1 µ m, which set the surface area available forchemistry on grains. Assuming a single grains size for the chem-istry is an approximation, but the e ff ect of grain growth in disks,which decreases this area, is cancelled by the dust settling in Article number, page 4 of 19. Leemker et al.: Chemically tracing the water snowline in protoplanetary disks with HCO + disks, which increases this area as there is more dust availablefor chemistry in the midplane (Eistrup et al. 2016). Finally, theH CO + abundance is taken to be a factor 70 smaller than theabundance of HCO + , corresponding to the typical C / C iso-tope ratio (Milam et al. 2005).
Predictions for multiple transitions of HCO + and H CO + weremade using the ray tracer in DALI , where the outputs of thechemical networks discussed in Section 2.2 were used to setthe abundance and both line and continuum optical depth are in-cluded. The line radiative transfer does not assume LTE. Instead,the excitation is calculated explicitly, where we used the col-lisional rate coe ffi cients in the LAMDA database (Botschwinaet al. 1993; Flower 1999; Schöier et al. 2005).Furthermore, the disk is assumed to be located at a distanceof 100 pc and the inclination is taken to be 38 ◦ . Changing theinclination does not significantly change our results. The radialemission profiles are taken to be along the semi-major axis ofthe disk. The focus in the paper is on the J = − + and H CO + respec-tively, because this transition provides the best balance betweenthe brightness of the line and the e ff ects of the continuum op-tical depth. Moreover, the higher transitions could be blendedwith emission from complex organic molecules. Predictions forthe weaker J = − J = − J = − ff ected by the opticaldepth of the continuum emission and emission from complexorganic molecules, can be found in Appendix C. To mimic highresolution ALMA observations, the emission is convolved witha 0 (cid:48)(cid:48) .
05 beam, unless denoted otherwise.
3. Modelling results
The predicted abundances of HCO + by chemical networks NWand W are presented in Fig. 3. The water snow surface in net-work W is computed as the surface where 50% of the total waterabundance is in the gas-phase and 50% is frozen-out onto thegrains and is indicated by the black line. In this model, the watersnowline, the midplane location where water freezes-out, is lo-cated at 4.5 AU. Comparing the two chemical networks outsidethe water snow surface shows a great similarity in the predictedHCO + abundances. This is due to the fact that water is frozen outoutside the water snow surface in network W, hence there is verylittle gas-phase water available for the destruction of HCO + . Innetwork W, the HCO + abundance drops inside the snowline, andis in general at least two orders of magnitude lower than in net-work NW. Up to 2 AU outside the water snowline HCO + is stille ffi ciently destroyed by the small amount of water that is presentin the gas-phase. Similar e ff ects were found for N H + and CO(Aikawa et al. 2015; van ’t Ho ff et al. 2017). As HCO + is e ffi -ciently destroyed by gas-phase water, it thus is a good chemicaltracer of the water snowline.The morphology of the HCO + distribution predicted bychemical network W is similar to the morphology predicted byfull chemical networks (Walsh et al. 2012, 2013; Agúndez et al.2018). Chemical network W agrees with the full chemical net-works quantitatively inside the water snowline where all studiespredict a low HCO + abundance of (cid:46) − . Moreover, the fullnetworks and network W all show a jump of at least one or-der of magnitude in the HCO + abundance over a radial range Fig. 3: HCO + abundance calculated using chemical networkNW (top) and network W (bottom). The water snow surface ismarked with the black line in the bottom panel. The position ofthe star in indicated with the symbol of a black star.of 5 AU outside the water snowline. In addition, these modelspredict a layer starting at z / r ∼ . + abundancereaches a high abundance of 10 − − − (Walsh et al. 2012,2013 and Agúndez et al. 2018). Yet, this layer is not expected tocontribute much to the column density and emission of HCO + ,because the densities are low in this layer. The midplane abun-dances at 10 AU are ∼ − in the full networks and 10 − in network W, though the gradient in the HCO + abundance isvery steep between the water snowline at 4.5 AU and a radiusof 10 AU in network W complicating comparison. An HCO + abundance of 10 − is reached at 7 AU. The HCO + abundanceat 100 AU lies within the range of HCO + abundances predictedby the full chemical models. Our small network is thus suited tostudy HCO + as a tracer of the water snowline. HCO + only acts as a tracer of the water snowline if the destruc-tion of HCO + by electrons is not dominant over the destructionby gas-phase water. The dissociative recombination of HCO + with an electron is included in both chemical networks. The ratesof the reactions in Eq. 1 and Eq. 5, R H O and R e − , respectively, Article number, page 5 of 19 & A proofs: manuscript no. Leemker_et_al_accepted
Fig. 4: Relative reaction rates of HCO + destruction by electrons( R e − ) and water ( R H O ). Gas-phase water is the dominant destruc-tion mechanism of HCO + in the blue regions and electrons arethe dominant destruction mechanism in the red regions. The wa-ter snow surface is indicated by the black line.are given by: R H O = k H O n (HCO + ) n (H O) and (6) R e − = k e − n (HCO + ) n (e − ) . (7)Using Table 2, the ratio of these reaction rates is given by: R e − R H O = (cid:32) T gas (cid:33) − . n (e − ) n (H O) . (8)Therefore, gas-phase water is the dominant destruction mecha-nism of HCO + if the number density of gas-phase water is abouttwo orders of magnitude larger than the number density of elec-trons.Fig. 4 shows this ratio for chemical network W. The re-gion where electrons are the dominant destruction mechanism( R e − > R H O ) is indicated in red and the region where wateris dominant ( R e − < R H O ) is indicated in blue. The first regionexists mostly outside the water snow surface as water is frozenout. The latter region exists mostly inside and above the watersnow surface, because there is plenty of water to destroy HCO + in these regions. Therefore HCO + is a good chemical tracer ofthe water snowline even though electrons are its dominant de-struction mechanism in most of the disk. O abundance, and cosmic rayionization rate on the HCO + abundance In this section we investigate the choice of initial conditions forchemical network W listed in Table B.1. The HCO + abundanceoutside the water snow surface in chemical network W can beapproximated analytically (for details see Appendix B.1): x (HCO + ) = (cid:114) ζ c . r . k e − n (H ) . (9)Based on this equation, it is expected that the abundance ofHCO + scales as the square root of the cosmic ray ionization ratein the disk region outside the water snowline. This analyticalprediction is consistent with the predictions made by the numer-ical solution to chemical network W (top left panel in Fig. B.1), where the HCO + column density decreases by one order of mag-nitude if the cosmic ray ionization rate decreases by two ordersof magnitude.Furthermore, the HCO + abundance is expected to be inde-pendent of the initial CO and H O abundance outside the watersnow surface based on the analytical approximation. This is inline with the predictions by the numerical solution to chemicalnetwork W (see HCO + column density in Fig. B.1). The onlymajor di ff erence between the analytical approximation and thenumerical solution occurs inside the water snowline. If the initialabundance of gas-phase water increases up to ∼ × − as ex-pected for water ice in cold clouds, the column density of HCO + inside the water snowline decreases. This is because the higherabundance of gas-phase water destroys more HCO + . There is nosignificant e ff ect on the HCO + column density outside the wa-ter snowline as all water is frozen out in that region of the disk.In summary, HCO + shows a steep jump in its column densityaround the water snowline for all initial conditions discussed inthis section. Therefore, the results in the following sections donot depend critically on the choice of initial conditions. Furtherdetails on the initial conditions can be found in Appendix B.2. + and H CO + emission Observations do not trace the local abundances, only the emis-sion which is closely related to the column densities if the emis-sion is optically thin. The column densities of HCO + and thepredicted radial emission profiles of HCO + and H CO + in net-work NW and W are shown in Fig. 5. Similar to the abundanceplots of HCO + , the total HCO + column density jumps by a fac-tor ∼
230 over a radial range of 3 AU outside the water snowlinein network W (solid blue line, top panel). Therefore, the HCO + column density shows a clear dependence on the presence of gas-phase water and the high HCO + abundance in the surface layersof the disk does not contribute much to the column density.On the other hand, the shapes of the radial emission profilesof the HCO + J = − ff erent fromthe HCO + column densities where only network W predicts theHCO + column density to peak outside the water snowline. Thedrop in the column density around the water snowline translatesinto a relative maximum di ff erence of only 24% between the ra-dial emission profiles of the J = − + . Thismaximum relative di ff erence is defined as the maximum di ff er-ence between the flux predicted by network NW and W com-pared to the maximum flux predicted by network W and is typi-cally located at the water snowline.HCO + is more abundant than H CO + , hence H CO + emis-sion is expected to be optically thin, whereas the HCO + emissionis optically thick. To test if this a ff ects the ability to trace thewater snowline, we predict emission from the less abundant iso-topologue H CO + . The radial emission profiles of H CO + forboth chemical networks are shown in the bottom panel of Fig. 5.Qualitatively the profiles are similar to the predicted emissionfor HCO + . Though the expected flux is about 39 times lowerfor H CO + than for HCO + . This di ff erence in flux is less thana factor 70, which indicates that the HCO + emission is indeedat least partially optically thick. The relative maximum di ff er-ence between chemical network NW and W is 13% for H CO + ,which is almost half of the corresponding number for HCO + . Sothe e ff ect of the water snowline is also seen in the emission ofthe optically thin H CO + , but even less prominently. The rel-ative di ff erence for H CO + is smaller than for HCO + because Article number, page 6 of 19. Leemker et al.: Chemically tracing the water snowline in protoplanetary disks with HCO + Fig. 5: Top panel: total HCO + column densities (top panel). Mid-dle and bottom panel: HCO + J = − CO + J = − ff ects of continuum opticaldepth (see Section 3.2.1). The water snowline is indicated withthe black dashed line. The radial emission profiles are convolvedwith a 0 (cid:48)(cid:48) .
05 beam. H CO + emits from a lower layer in the disk due to its lower op-tical depth. Therefore more dust is present above the layer whereH CO + emits, hence the e ff ect of dust is larger for H CO + thanfor HCO + .These models show very similar radial emission profiles pre-dicted for network NW and W and only a small di ff erence influx. Both networks predict ring shaped emission, regardless ofthe presence of the water snowline. This decrease in line flux to-wards the center of the disk can have di ff erent origins. The firstone is a decrease in the abundance of the observed molecule,discussed in the previous section, which is the e ff ect we aim toobserve. However, the optical depth of the continuum emissionand molecular excitation e ff ects can change the emission as well,potentially obscuring the e ff ect we would like to trace. The continuum a ff ects the line intensity when the optical depthof the continuum emission is larger than the optical depth of theline or when both are optically thick (Isella et al. 2016). To in-vestigate the e ff ect of the continuum emission, the J = − + and H CO + was ray traced in a disk modelwhere the dust opacities (dust mass absorption coe ffi cients κ ν )are divided by a factor of 10 , to cancel the e ff ect of dust on theHCO + emission compared to the typical Herbig Ae disk model.The results are shown as the dotted lines in Fig. 5. Comparingthese radial emission profiles with the fiducial model shows thatthe e ff ect of dust is small for network W in the inner ∼ ∼
14 AU (seethe orange line in Fig. 6). The reason for this is that the HCO + abundance is low in the inner disk for network W. Therefore, lit-tle emission is expected, hence there is also little emission thatcan be a ff ected by the dust.However, the dust can mimic the e ff ect of gas-phase water,which greatly complicates matters. That the e ff ect of the dust issignificant, can be seen when the radial emission profiles pre-dicted by network NW in the fiducial disk model and in themodel with very low dust opacities (solid versus dotted red lines)are compared. This shows that the dust optical depth is largelyresponsible for the drop in emission in the center in network NWand for the fact that the radial emission profiles for network NWand W look very similar in the fiducial disk model. Therefore, tolocate the water snowline using HCO + , a disk with a high gas-to-dust mass ratio would be most suitable. However, recent workby Kama et al. (2020) has shown that disks with high gas-to-dustmass ratios around Herbig Ae stars are rare. Another option isto target HCO + in even warmer sources than Herbig Ae stars.Fig. 6 shows that in more luminous sources, such as outburstingHerbig sources, the snowline is expected to shift to radii wherethe dust is optically thin at the wavelengths of the J = − J = − + . Note however that the τ dust = Decreasing the dust opacities by a factor of 10 cannot fullyexplain the unexpected ring-shaped emission found by networkNW. A second e ff ect that contributes to the decrease in flux in the Article number, page 7 of 19 & A proofs: manuscript no. Leemker_et_al_accepted
Fig. 6: τ dust = + J = − J = − J = − J = − J -level populationsof HCO + . The column densities and midplane populations ofseveral J -levels of HCO + in chemical network NW are shownin Fig. 7. The column density of the J = + inchemical network NW peaks at a radius of ∼
30 AU, while thetotal HCO + column density, N tot , peaks on-source (red line intop panel of Fig. 5).The radius where the column density of the J = J = − + associated with network NW with lowdust opacities starts to drop. Therefore, both the optical depth ofthe continuum emission as well as the temperature dependenceof the J -level population contribute to the drop in emission seenin the center. Finally there could also be a chemical e ff ect in network NW asthe midplane abundance of HCO + decreases towards the centerof the disk. Yet, the total column density of HCO + keeps in-creasing up to a radius of 0.07 AU, which is too small to explainthe decrease in flux out to ∼
50 AU seen in the radial emissionprofile associated with network NW in the fiducial model or to ∼ + column density. + rings? In summary, all three e ff ects contribute to the decrease in flux inthe inner regions of the disk. Observing di ff erent transitions ofHCO + and H CO + will give a di ff erent balance of these e ff ects.Fig. 7 shows that higher J -levels are more populated in the innerdisk because the inner disk is too warm and dense to highly pop-ulate the J = + . However, the emitting frequencyof HCO + increases as higher transitions are observed, so also thecontinuum optical depth increases. Fig. 7: HCO + column density (top) and midplane populations(bottom) for the J = J = J = J = ffi cult, yet crucial, to disentangle the ef-fects of the water snowline and of the population of HCO + andthe optical depth of the continuum emission. The e ff ect of thecontinuum optical depth can be quantified by observing anothermolecule that emits from the same region as HCO + or H CO + but does not decrease in column density towards the star such asa CO isotopologue. The excitation and column density of HCO + need to be inferred from detailed modelling. In conclusion, theresults in this section show that even though the HCO + abun-dance changes by at least two orders of magnitude around thewater snow surface, it is observationally complicated to verifythis. The discussion above has shown that locating the water snow-line using the HCO + radial emission profile is di ffi cult. The ve-locity resolved line profile provides an alternative method as theremoval of HCO + in the inner disk is expected to remove fluxat the highest Keplerian velocities. Accordingly, the line profilepredicted by chemical network W is expected to be narrowerthan the line profile predicted by network NW. Article number, page 8 of 19. Leemker et al.: Chemically tracing the water snowline in protoplanetary disks with HCO + Fig. 8: Left panels: line profile for the J = − + (top) and H CO + (bottom) for chemical network NW(red) and W (light blue). Right panels: di ff erence between theline profiles in the left column. The black arrows indicate wherethe e ff ect of the water snowline is expected. Note the di ff erencesin the vertical axes.The line profiles predicted by the two chemical networks forthe J = − + and H CO + are shown in leftcolumn of Fig. 8. To show the di ff erence between the modelswe subtracted the spectrum associated with network W from thespectrum associated with network NW. The result is shown inthe right column of Fig. 8. This column shows that even thoughthe line profiles look almost identical in the left column, thereis an additional bump in the di ff erence between them in the linewings (indicated with the black arrows). However, this di ff erenceis only 1.75 mJy for HCO + and 28 µ Jy for H CO + i.e. 0.6%and 0.3% of emission, which cannot be detected with a reason-able signal-to-noise ratio. Moreover, detailed comparison withthe line profile associated with a molecule that does not show adecrease in the column density in the inner disk would be neededas only the flux associated with network W would be observedand not the flux associated with network NW. Hence, it is verydi ffi cult to use the line profiles to locate the water snowline indisks around Herbig Ae stars.
4. Comparison with observations
The best targets to observe the water snowline using HCO + arewarm disks, because these disks have a water snowline at a largeenough radius for ALMA to resolve. In addition, disks withoutdeep gaps in the gas and dust surface density around the expectedlocation of the snowline are preferred, as these gaps complicatethe interpretation of the HCO + emission.One of the most promising targets to locate the water snow-line in a disk is the young outbursting star V883 Ori, with a current luminosity of ∼
218 L (cid:12) (Furlan et al. 2016). During theoutburst, the luminosity is greatly increased, which heats thedisk, shifting the water snowline outwards compared to regularT Tauri and even Herbig Ae disks. This not only allows for ob-servations with a larger beam, it also mitigates the e ff ects of thedust as the optical depth of the continuum emission decreaseswith radius. Previous observations of V883 Ori have inferred thelocation of the water snowline at 42 AU based on a change inthe dust emission (Cieza et al. 2016). An extended hot inner re-gion is consistent with the detection of many complex organicmolecules (Lee et al. 2019). However, an abrupt change in thecontinuum optical depth or the spectral index are not necessarilydue to a snowline, as shown for the cases of CO , CO and N (Huang et al. 2018; Long et al. 2018; van Terwisga et al. 2018).Methanol observations suggest that the water snowline may belocated at a much larger radius of ∼
100 AU if the methanol ob-servations are optically thin and most of its emission originatesfrom inside the water snowline (van ’t Ho ff et al. 2018b). Obser-vations of HCO + have the potential to resolve this discrepancy. CO + observations A promising dataset to locate the water snowline in V883Ori is the band 7 observation by Lee et al. (2019) (projectcode: 2017.1.01066.T, PI: Jeong-Eun Lee), which contains theH CO + J = − CO + observations in a protoplanetary disk, to our knowledge,at su ffi cient spatial and spectral resolution with su ffi cient sensi-tivity to potentially resolve the water snowline. The synthesizedbeam of these observations is 0 (cid:48)(cid:48) . ∼
400 pc; Kounkel et al. 2017), and the spec-tral resolution is 0.25 km s − . The continuum is created using allline-free channels, which are carefully selected to exclude anyline emission. The line data are continuum subtracted using thesecontinuum solutions. Using the CASA 5.1.1 tclean task with aBriggs weighting of 0.5, line images are made with a mask ofabout 2" in diameter centered on the peak of the continuum.Following the approach of Lee et al. (2019), the line obser-vations are velocity-stacked using a stellar mass, M (cid:63) = (cid:12) ,an inclination of 38 ◦ , and a position angle of 32 ◦ (Cieza et al.2016). This technique reduces line blending because it makesuse of the Keplerian velocity of the disk to calculate the Dopplershift of the emission in each pixel (Yen et al. 2016, 2018). ThisDoppler shift is then used to shift the spectrum of each pixel tothe velocity of the star, before adding the spectra of all pixels.The noise level of ∼
15 mJy per spectral bin is determined usingan empty region in the stacked image.The observed spectrum of the inner 0 (cid:48)(cid:48) . CO + J = − , - 17 , , E transition of acetalde-hyde (CH CHO) ( v =
0) at 346.9955 GHz and the two super-imposed 18 , - 17 , , E and 18 , - 17 , , E transitions ofCH CHO ( v T =
2) at 346.99991 GHz and 346.99994 GHz, re-spectively (Jet Propulsion Laboratory (JPL) molecular database;Pickett et al. 1998). In between these acetaldehyde transitions, aclear excess in emission is visible due to the J = − CO + . To quantify this excess emission, two Gaussianprofiles are fitted with curve_fit to model the acetaldehydeemission. To reduce the number of free parameters in the fit,the line frequencies are fixed to their respective rest frequencieslisted above and the widths of the lines are fixed to 2 km s − fol- The curve_fit function is part of the Scipy package in python.Article number, page 9 of 19 & A proofs: manuscript no. Leemker_et_al_accepted
Fig. 9: Top panel: observed and stacked spectrum of V883Ori (black) and a fit of the CH CHO emission (red); bottompanel: spectrum where the CH CHO emission is subtracted. TheH CO + J = − CHO lines (vertical pink lines). The 3 σ noise level is indicated with the horizontal grey line.lowing Lee et al. (2019). In addition, emission within 1 km s − from the rest frequency of the J = − CO + is excluded from the fit. The result is shown as the red line inthe top panel of Fig. 9. Subtracting this fit to the acetaldehydeemission from the observed spectrum clearly reveals emissionof the J = − CO + . Therefore we concludethat H CO + is detected at > σ but is blended with lines that areexpected to peak on-source, preventing a clean image.The line blending of H CO + J = − CO + at su ffi cient spatial and spectral resolution andsensitivity to resolve the water snowline. The J = − CO + is likely blended with CH OCHO as both of them arebright in the Class 0 source B1-c (van Gelder et al. 2020). How-ever, the H CO + J = − ff et al. in prep.) andwould therefore be the best line to target in sources with brightemission from complex organics. + observations Another promising dataset for HCO + is a recent band 6 dataset(project code 2018.1.01131.S, PI: D. Ruíz-Rodríguez; Ruíz-Rodríguez et al. in prep.). The J = − + is detected with a high signal-to-noise ratio. The spatial resolu-tion of ∼ (cid:48)(cid:48) . ffi cient to resolve the water snowline if it islocated within a radius of ∼
90 AU (180 AU diameter).An image of the HCO + emission will be presented in Ruíz-Rodríguez et al. (in prep.). Here, a normalized, deprojected az-imuthal average of the observed emission from the product dataof the J = − + is presented in the top panelof Fig. 10. This azimuthal average shows that the HCO + emis- Fig. 10: Top: normalized azimuthal average of the observedHCO + J = − -17 flux (red; van’t Ho ff et al. 2018b). Bottom: normalized azimuthal average ofthe observed HCO + flux (black; as in top panel) and modelledHCO + flux for a snowline at 47 AU (blue), 76 AU (orange) and119 AU (green). The snowline locations of the models are indi-cated with dashed lines in corresponding colors. The position ofthe star is indicated with a black star and the beam is indicatedby the black bar in the bottom right corner.sion is ring shaped. Comparing the band 6 HCO + emission withthe band 7 emission from complex organic molecules presentedin van ’t Ho ff et al. (2018b) and Lee et al. (2019) , reveals thatthe emission from complex organic molecules such as methanol,13-methanol, acetaldehyde and H C O peak inside the HCO + ring and in some cases is even centrally peaked. This is a strongindication that the lack of HCO + emission in the center is notonly due to the optical depth of the continuum and that HCO + is indeed tracing the water snowline in V883 Ori. Such an anti-correlation between HCO + and methanol has also been seen inprotostellar envelopes (Jørgensen et al. 2013). + images: locating the snowline To quantify the location of the water snowline, we model theHCO + and H CO + emission in a representative model for V883Ori. This model reproduces the previously observed flux of the Article number, page 10 of 19. Leemker et al.: Chemically tracing the water snowline in protoplanetary disks with HCO + Fig. 11: Integrated intensity maps for the J = − + predicted by network NW (top row) and network W (bottomrow). The snowline is indicated with a white ellipse and is located at 47 AU (left column), 76 AU (middle column) and 119 AU(right column). The cartoons above the individual panels provide a sketch of the model, where blue indicates water ice, red indicatesgas-phase water, orange indicates a high abundance of HCO + and gray indicates a low abundance of HCO + . The position of the staris marked with a yellow star in each cartoon and with white star in each panel. Note that the size of the star increases with increasingluminosity of the star. The 0 (cid:48)(cid:48) . J = − O (van ’t Ho ff et al. 2018b), the HCO + and H CO + total flux discussed in the previous Sections, andmm continuum fluxes within a factor of ∼
2. The
DALI modelparameters are presented in Table D.1. The main changes com-pared to the model for the typical Herbig Ae disk include themass and radius of the disk and the mass and luminosity of thestar. The characteristic disk radius is set to 75 AU to match theradial extent of emission from the J = − Opresented by van ’t Ho ff et al. (2018b). The disk mass is esti-mated using the 9.1 mm continuum observations of the VAN-DAM survey (Tobin et al. 2020) and the relation between thecontinuum flux and the disk mass (Hildebrand 1983): M = D F ν κ ν B ν ( T dust ) , (10)with D the distance to V883 Ori, F ν the observed continuumflux, κ ν the opacity and B ν ( T dust ) the Planck function for a dusttemperature T dust . As V883 Ori is an outbursting source, a dusttemperature of 50 K is assumed. Following Tychoniec et al. (2020), a dust opacity of 0.28 cm g − at a wavelength of 9.1 mmis used (Woitke et al. 2016). This gives an estimated disk massof 0.25 M (cid:12) .The snowline has been estimated at 42 AU from dust (Ciezaet al. 2016), but can be as far out as 100 AU based on CH OH(van ’t Ho ff et al. 2018b). Moreover, the outburst likely beganbefore 1888 (Pickering 1890). Approximately 25 years ago theluminosity was measured to be ∼
200 L (cid:12) higher than the currentluminosity (Strom & Strom 1993; Sandell & Weintraub 2001;Furlan et al. 2016) and it could have been much higher in thepast. With a freeze-out time scale of 100-1000 yrs, the snow-line may thus not be at the location expected from the currentluminosity (Jørgensen et al. 2013; Visser et al. 2015; Hsieh et al.2019). We therefore use three luminosities of 2 × , 6 × and 1 . × L (cid:12) in our models, which put the snowline at 47,76, and 119 AU, respectively.The integrated intensity maps of the J = − + for these models are shown in Fig. 11. The correspondingfigures for the weaker HCO + and H CO + J = − Article number, page 11 of 19 & A proofs: manuscript no. Leemker_et_al_accepted are presented in Fig. D.1 and Fig. D.2. These lines are not ex-pected to be contaminated by emission from complex organicmolecules. The results in this section are convolved with a small0 (cid:48)(cid:48) . ff erent luminosities. The total HCO + emission does not depend strongly on the luminosity becausethe population of the J = ff ect as the HCO + emission is marginally optically thick.Similar to the models for a disk around a typical Herbig Ae star,the moment 0 maps of network NW show ring shaped emissiondespite the fact that there is no snowline present in these models.The lack of emission in the center is dominated by absorptionby dust. On top of that, the column density of the J = + decreases in the inner parts of the disk. These two pointsare di ffi cult to disentangle because the continuum optical depthis frequency and hence J -level dependent. However, most im-portantly, the location of the HCO + ring does not change as afunction of luminosity in network NW.This is di ff erent in the bottom row of Fig. 11 where emissionassociated with network W is shown. The water snowline is indi-cated with a white ellipse and shifts outwards as the luminosityof the star increases. This is also reflected in the HCO + emis-sion as the ring of HCO + shifts outwards together with the watersnowline. Therefore, the location of the HCO + ring can be usedto locate the water snowline provided that another molecule isobserved to prove that the decrease in the HCO + flux in the cen-ter is not solely due to molecular excitation e ff ects and absorp-tion by dust. The excitation e ff ect needs to be inferred from diskmodelling. The e ff ect of the optical depth of the continuum emis-sion can be estimated by observing another molecule, that emitsfrom the same disk region as HCO + and whose column den-sity does not drop in the inner disk. The most obvious moleculefor this purpose is C O. However, the HCO + J = − O J = − J = − ff ect of the duston the HCO + emission cannot be fully traced by the J = − J = − O. More rare isotopologues such asC O or C O are more suited for this purpose, as well as othermolecules such as complex organics as discussed in Section 4.2.The model predictions for the observed flux of the J = − + in V883 Ori are compared with the observedflux in the bottom panel of Fig. 10. The model results are con-volved to the same spatial resolution as the observations beforecalculating the deprojected azimuthal average. Even though thebeam is too large to resolve the snowline if it is located inside ∼
90 AU, this figure clearly shows that the observed HCO + fluxdrops o ff steeper in the inner parts than that predicted by themodel with a snowline at 47 AU. In addition, the HCO + ringshifts outwards with increasing snowline location. The observedpeak location and gap depth in the center best match with asnowline between 76 and 119 AU.A central cavity (approx. 1 beam in diameter at 0 (cid:48)(cid:48) . × (cid:48)(cid:48) . ∼
60 AU resolution in radius) wasalso observed for CO J = − + is solely due tothe continuum as the emission from complex organic moleculespeaks inside the HCO + ring (0 (cid:48)(cid:48) . × (cid:48)(cid:48) .
17; Lee et al. 2019). Simi-larly, the high resolution methanol emission (0 (cid:48)(cid:48) . × (cid:48)(cid:48) .
14; van ’tHo ff et al. 2018b) and the continuum emission (0 (cid:48)(cid:48) .
03; Cieza et al.2016) peak well inside the HCO + cavity (Fig. 10). Furthermore, the continuum becomes optically thin well within 100 AU in allthree of our models (Figs. 10 and 11), so the di ff erence betweenthe models is due to di ff erent snowline locations. Therefore, ouranalysis is consistent with a water snowline at ∼ ff ect that needs to be taken into account in outburstingsources like V883 Ori is viscous heating. This e ff ect heats thedisk midplane, hence the water snowline could be at a largerradius than expected based on radiative heating alone. The J = − + is marginally optically thick in ourmodels. Targeting a low J transition of HCO + or targeting theoptically thin H CO + , will allow to directly trace the midplane.Another e ff ect that could occur due to viscous heating is self-absorption of HCO + as the surface layers of the disk could havea lower temperature than the viscously heated midplane. Thiscould be resolved by observing multiple transitions of HCO + ashigher J lines are more optically thick and hence trace a layerhigher in the disk.Taken together, our analysis of the observations suggests thatHCO + is tracing the water snowline in V883 Ori and that thesnowline is located outside the radius where the change in thecontinuum opacity is observed. Therefore, V883 Ori is the firstexample that shows that a change in the dust continuum opac-ity is not necessarily related to the water snowline, similar towhat has been shown for the cases of CO, CO and N snow-lines (Huang et al. 2018; Long et al. 2018; van Terwisga et al.2018).
5. Conclusions
Chemical imaging with HCO + is in principle a promisingmethod to image the water snowline as it has been used suc-cesfully in a protostellar envelope (van ’t Ho ff et al. 2018a). Thiswork examines its application in older protoplanetary disks. TheHCO + abundance is modelled using two chemical networks. Thefirst one, network NW, does not include reactions with water, incontrast to the second one, network W. Predictions for the radialemission profiles of HCO + and H CO + are made to examinethe validity of HCO + as a tracer of the water snowline. More-over, archival observations of V883 Ori are examined and usedto put constraints on the water snowline location in V883 Oriand make predictions for future high resolution observations ofHCO + in the disk around this outbursting star.Based on our models the following conclusions can bedrawn: • The HCO + abundance jumps two orders of magnitudearound the water snowline. • In addition to the water snowline, the optical depth of thecontinuum emission and molecular excitation e ff ects for thelow J -levels contribute significantly to the decrease in theH CO + and HCO + flux in the inner parts of the disk andresult in ring shaped emission. Therefore the e ff ect of thecontinuum optical depth needs to be checked observationallyand the e ff ects of the molecular excitation and HCO + abun-dance need to be modelled in detail. Outbursting sources arethe best targets, as the snowline is shifted to larger radii,where the dust optical depth is lower. • HCO + and H CO + are equally good as a tracers of the watersnowline but the main isotopologue HCO + is more readilyobservable. • For both HCO + and H CO + , the J = − ff ects of the continuum optical depth. Moreover, Article number, page 12 of 19. Leemker et al.: Chemically tracing the water snowline in protoplanetary disks with HCO + it is not expected to be blended with emission from complexorganic molecules, unlike the J = − J = − CO + . • Our analysis of the observations of the HCO + J = − + istracing the water snowline in V883 Ori. Based on our mod-els, the snowline is located around 100 AU and is not cor-related with the opacity change in the continuum emissionobserved at 42 AU.Our results thus show that HCO + and H CO + can be usedto trace the water snowline in warm protoplanetary disks suchas those found around luminous stars. Determining the snowlinelocation in these sources is important to understand the processof planet formation and composition. Acknowledgements.
We thank the referee and the editor for the construc-tive comments, and Catherine Walsh and Jeong-Eun Lee for useful discus-sions. Astrochemistry in Leiden is supported by the Netherlands ResearchSchool for Astronomy (NOVA). M.L.R.H acknowledges support from a Huy-gens fellowship from Leiden University, and from the Michigan Society ofFellows. L.T. is supported by NWO grant 614.001.352. This paper makesuse of the following ALMA data: ADS / JAO.ALMA / JAO.ALMA / NRAO and NAOJ.
References
Adams, N. G., Smith, D., & Grief, D. 1978, International Journal of Mass Spec-trometry and Ion Processes, 26, 405Agúndez, M., Roue ff , E., Le Petit, F., & Le Bourlot, J. 2018, A&A, 616, A19Aikawa, Y., Furuya, K., Nomura, H., & Qi, C. 2015, ApJ, 807, 120Andrews, S. M. 2020, arXiv e-prints, arXiv:2001.05007Andrews, S. M., Wilner, D. J., Espaillat, C., et al. 2011, ApJ, 732, 42Andrews, S. M., Wilner, D. J., Hughes, A. M., Qi, C., & Dullemond, C. P. 2010,ApJ, 723, 1241Anicich, V. G., Futrell, J. H., Huntress, Wesley T., J., & Kim, J. K. 1975, Inter-national Journal of Mass Spectrometry and Ion Processes, 18, 63Banzatti, A., Pinilla, P., Ricci, L., et al. 2015, ApJ, 815, L15Banzatti, A., Pontoppidan, K. M., Salyk, C., et al. 2017, ApJ, 834, 152Bergin, E. A., Hogerheijde, M. R., Brinch, C., et al. 2010, A&A, 521, L33Bergin, E. A., Melnick, G. J., & Neufeld, D. A. 1998, ApJ, 499, 777Blevins, S. M., Pontoppidan, K. M., Banzatti, A., et al. 2016, ApJ, 818, 22Booth, A. S. & Ilee, J. D. 2020, MNRAS, 493, L108Booth, A. S., Walsh, C., Ilee, J. D., et al. 2019, ApJ, 882, L31Bosman, A. D., Tielens, A. G. G. M., & van Dishoeck, E. F. 2018, A&A, 611,A80Botschwina, P., Horn, M., Flügge, J., & Seeger, S. 1993, J. Chem. Soc., FaradayTrans., 89, 2219Bruderer, S. 2013, A&A, 559, A46Bruderer, S., Doty, S. D., & Benz, A. O. 2009, ApJS, 183, 179Bruderer, S., van Dishoeck, E. F., Doty, S. D., & Herczeg, G. J. 2012, A&A, 541,A91Bryden, G., Chen, X., Lin, D. N. C., Nelson, R. P., & Papaloizou, J. C. B. 1999,ApJ, 514, 344Carney, M. T., Fedele, D., Hogerheijde, M. R., et al. 2018, A&A, 614, A106Carr, J. S. & Najita, J. R. 2008, Science, 319, 1504Caselli, P., Walmsley, C. M., Terzieva, R., & Herbst, E. 1998, ApJ, 499, 234Cazzoletti, P., van Dishoeck, E. F., Visser, R., Facchini, S., & Bruderer, S. 2018,A&A, 609, A93Cieza, L. A., Casassus, S., Tobin, J., et al. 2016, Nature, 535, 258D’Alessio, P., Calvet, N., Hartmann, L., Franco-Hernández, R., & Servín, H.2006, ApJ, 638, 314Dong, R., Liu, S.-y., Eisner, J., et al. 2018, ApJ, 860, 124Dr˛a˙zkowska, J. & Alibert, Y. 2017, A&A, 608, A92Du, F., Bergin, E. A., Hogerheijde, M., et al. 2017, ApJ, 842, 98Eistrup, C., Walsh, C., & van Dishoeck, E. F. 2016, A&A, 595, A83Eistrup, C., Walsh, C., & van Dishoeck, E. F. 2018, A&A, 613, A14Fedele, D., Carney, M., Hogerheijde, M. R., et al. 2017, A&A, 600, A72Flower, D. R. 1999, MNRAS, 305, 651Fraser, H. J., Collings, M. P., McCoustra, M. R. S., & Williams, D. A. 2001,MNRAS, 327, 1165 Furlan, E., Fischer, W. J., Ali, B., et al. 2016, ApJS, 224, 5Harsono, D., Bruderer, S., & van Dishoeck, E. F. 2015, A&A, 582, A41Harsono, D., Persson, M. V., Ramos, A., et al. 2020, A&A, 636, A26Hartmann, L., Calvet, N., Gullbring, E., & D’Alessio, P. 1998, ApJ, 495, 385Hasegawa, T. I., Herbst, E., & Leung, C. M. 1992, ApJS, 82, 167Hildebrand, R. H. 1983, QJRAS, 24, 267Hogerheijde, M. R., Bergin, E. A., Brinch, C., et al. 2011, Science, 334, 338Hsieh, T.-H., Murillo, N. M., Belloche, A., et al. 2019, ApJ, 884, 149Huang, J., Andrews, S. M., Dullemond, C. P., et al. 2018, ApJ, 869, L42Isella, A., Guidi, G., Testi, L., et al. 2016, Phys. Rev. Lett., 117, 251101Jørgensen, J. K., Visser, R., Sakai, N., et al. 2013, ApJ, 779, L22Kama, M., Bruderer, S., van Dishoeck, E. F., et al. 2016, A&A, 592, A83Kama, M., Trapman, L., Fedele, D., et al. 2020, A&A, 634, A88Kim, J. K., Theard, L. P., & Huntress, W. T., J. 1974, International Journal ofMass Spectrometry and Ion Processes, 15, 223Kimura, H., Wada, K., Kobayashi, H., et al. 2020, Monthly Notices of the RoyalAstronomical Society, staa2467Klippenstein, S. J., Georgievskii, Y., & McCall, B. J. 2010, Journal of PhysicalChemistry A, 114, 278Kounkel, M., Hartmann, L., Loinard, L., et al. 2017, ApJ, 834, 142Lee, J.-E., Lee, S., Baek, G., et al. 2019, Nature Astronomy, 3, 314Lepp, S., Dalgarno, A., & Sternberg, A. 1987, ApJ, 321, 383Long, F., Pinilla, P., Herczeg, G. J., et al. 2018, ApJ, 869, 17Lynden-Bell, D. & Pringle, J. E. 1974, MNRAS, 168, 603Mandell, A. M., Bast, J., van Dishoeck, E. F., et al. 2012, ApJ, 747, 92Mathews, G. S., Klaassen, P. D., Juhász, A., et al. 2013, A&A, 557, A132Mathis, J. S., Rumpl, W., & Nordsieck, K. H. 1977, ApJ, 217, 425McCall, B. J., Huneycutt, A. J., Saykally, R. J., et al. 2004, Phys. Rev. A, 70,052716McElroy, D., Walsh, C., Markwick, A. J., et al. 2013, A&A, 550, A36Milam, S. N., Savage, C., Brewster, M. A., Ziurys, L. M., & Wycko ff , S. 2005,ApJ, 634, 1126Mitchell, J. B. A. 1990, Phys. Rep., 186, 215Notsu, S., Akiyama, E., Booth, A., et al. 2019, ApJ, 875, 96Novotný, O., Buhr, H., Stützel, J., et al. 2010, Journal of Physical Chemistry A,114, 4870Öberg, K. I., Murray-Clay, R., & Bergin, E. A. 2011, ApJ, 743, L16Okuzumi, S., Momose, M., Sirono, S.-i., Kobayashi, H., & Tanaka, H. 2016,ApJ, 821, 82Phillips, T. G., van Dishoeck, E. F., & Keene, J. 1992, ApJ, 399, 533Pickering, E. C. 1890, Annals of Harvard College Observatory, 18, 1Pickett, H. M., Poynter, R. L., Cohen, E. A., et al. 1998,J. Quant. Spectr. Rad. Transf., 60, 883Pinilla, P., Pohl, A., Stammler, S. M., & Birnstiel, T. 2017, ApJ, 845, 68Qi, C., Öberg, K. I., Andrews, S. M., et al. 2015, ApJ, 813, 128Qi, C., Öberg, K. I., Espaillat, C. C., et al. 2019, ApJ, 882, 160Qi, C., Öberg, K. I., Wilner, D. J., et al. 2013, Science, 341, 630Ruíz-Rodríguez, D., Cieza, L. A., Williams, J. P., et al. 2017, MNRAS, 468, 3266Salinas, V. N., Hogerheijde, M. R., Bergin, E. A., et al. 2016, A&A, 591, A122Salyk, C., Lacy, J., Richter, M., et al. 2019, ApJ, 874, 24Salyk, C., Pontoppidan, K. M., Blake, G. A., et al. 2008, ApJ, 676, L49Sandell, G. & Weintraub, D. A. 2001, ApJS, 134, 115Schöier, F. L., van der Tak, F. F. S., van Dishoeck, E. F., & Black, J. H. 2005,A&A, 432, 369Schoonenberg, D. & Ormel, C. W. 2017, A&A, 602, A21Stevenson, D. J. & Lunine, J. I. 1988, Icarus, 75, 146Strom, K. M. & Strom, S. E. 1993, ApJ, 412, L63Teague, R., Guilloteau, S., Semenov, D., et al. 2016, A&A, 592, A49Teague, R., Semenov, D., Guilloteau, S., et al. 2015, A&A, 574, A137Theard, L. P. & Huntress, W. T. 1974, J. Chem. Phys., 60, 2840Tobin, J. J., Sheehan, P. D., Megeath, S. T., et al. 2020, ApJ, 890, 130Tychoniec, Ł., Manara, C. F., Rosotti, G. P., et al. 2020, A&A, 640, A19van der Marel, N., van Dishoeck, E. F., Bruderer, S., Pérez, L., & Isella, A. 2015,A&A, 579, A106van Gelder, M. L., Tabone, B., Tychoniec, Ł., et al. 2020, A&A, 639, A87van ’t Ho ff , M. L. R., Persson, M. V., Harsono, D., et al. 2018a, A&A, 613, A29van ’t Ho ff , M. L. R., Tobin, J. J., Trapman, L., et al. 2018b, ApJ, 864, L23van ’t Ho ff , M. L. R., Walsh, C., Kama, M., Facchini, S., & van Dishoeck, E. F.2017, A&A, 599, A101van Terwisga, S. E., van Dishoeck, E. F., Ansdell, M., et al. 2018, A&A, 616,A88van Zadelho ff , G. J., Aikawa, Y., Hogerheijde, M. R., & van Dishoeck, E. F.2003, A&A, 397, 789Visser, R. & Bergin, E. A. 2012, ApJ, 754, L18Visser, R., Bergin, E. A., & Jørgensen, J. K. 2015, A&A, 577, A102Vorobyov, E. I., Bara ff e, I., Harries, T., & Chabrier, G. 2013, A&A, 557, A35Walsh, C., Millar, T. J., & Nomura, H. 2013, ApJ, 766, L23Walsh, C., Nomura, H., Millar, T. J., & Aikawa, Y. 2012, ApJ, 747, 114Williams, J. P., Bergin, E. A., Caselli, P., Myers, P. C., & Plume, R. 1998, ApJ,503, 689Woitke, P., Min, M., Pinte, C., et al. 2016, A&A, 586, A103Yen, H.-W., Koch, P. M., Liu, H. B., et al. 2016, ApJ, 832, 204Yen, H.-W., Koch, P. M., Manara, C. F., Miotello, A., & Testi, L. 2018, A&A,616, A100Zhang, K., Blake, G. A., & Bergin, E. A. 2015, ApJ, 806, L7Zhang, K., Pontoppidan, K. M., Salyk, C., & Blake, G. A. 2013, ApJ, 766, 82Zhu, Z., Stone, J. M., Rafikov, R. R., & Bai, X.-n. 2014, ApJ, 785, 122 Article number, page 13 of 19 & A proofs: manuscript no. Leemker_et_al_accepted
Appendix A: Freeze-out and desorption coefficients
The reaction coe ffi cients of the freeze-out and desorption of wa-ter depend on the local conditions in the disk. The rate coe ffi cientof the freeze-out of water, k f can be expressed as: k f = (cid:104) v (cid:105) π a n (grain) S (A.1) = (cid:114) k B T gas m π a n (grain) S , (A.2)where (cid:104) v (cid:105) = (cid:112) k B T gas / m is the thermal velocity of gas-phasewater molecules with mass m in a gas at temperature T gas , and k B is the Boltzmann constant. The DALI models use a distribu-tion of grain sizes, but for the chemistry, only the surface areaof the grains matters. Therefore a grain size of a grain = . µ mand a number density, n (grain) = − × n (H ) with respect tomolecule hydrogen is used. Finally a sticking coe ffi cient, S , of 1is assumed.The reverse process, the desorption of water ice, is describedby the desorption rate coe ffi cient: k d = ν e − E b / k B T dust (A.3) = (cid:114) n s E b π m e − E b / k B T dust , (A.4)with ν = (cid:112) n s E b /π m the characteristic vibrational frequencyof water ice on a grain. Here n s = . × cm − is the numberdensity of surface sites where water can bind (Hasegawa et al.1992). The binding energy, E b , assumed in this work is 5775 K,corresponding to an amorphous water ice substrate (Fraser et al.2001). Appendix B: Chemical network
Appendix B.1: Analytical approximation for network NW
Chemical network NW can be used to derive an analytical ex-pression for the HCO + abundance because of its simplicity. Thetime evolution of the HCO + number density is given by the sumof the formation and destruction rates: dn (HCO + ) dt = k n (CO) n (H + ) − k e − n (HCO + ) n (e − ) , (B.1)with k and k e − the reaction rates as listed in Table 2. This can berewritten to n (HCO + ) = (cid:114) k k e − n (CO) n (H + ) (B.2)under the assumption of steady state. Furthermore, it is assumedthat metallic ions can be neglected as electron donors and thatHCO + is the main electron donor. HCO + has been found to bethe dominant molecular ion protoplanetary disks (Teague et al.2015) and the main charge carrier in starless cores (e.g. Caselliet al. 1998; Williams et al. 1998).Following the approach of Lepp et al. (1987), the numberdensity of H + , which is mainly governed by the ionization, canbe found in a similar way. The three reactions that regulate theionization in the disk: ionization by cosmic rays (reaction ζ c . r . ),the ion-molecule reaction to form H + (reaction k ) and dissocia-tive recombination of H + (reaction k ), see also Table. 2. Thefirst two reactions can be approximated as:H + c . r . → H + , (B.3) as the formation of H + is limited by the cosmic ray ionizationrate. Assuming CO is the main destroyer of H + and thus neglect-ing the dissociative recombination of H + , the time evolution andsteady state abundance of H + can be expressed as: dn (H + ) dt = ζ c . r . n (H ) − k n (CO) n (H + ) (time-dependent), and(B.4) n (H + ) = ζ c . r . n (H ) k n (CO) (steady state) . (B.5)Combining Eq. B.2 and Eq. B.5 gives the analytical approx-imation of the HCO + abundance in the disk: x (HCO + ) = (cid:114) ζ c . r . k e − n (H ) . (B.6)The comparison of chemical network NW and W in Section 3.1and Fig. 3 shows that the HCO + abundance in network NW isvery similar to the HCO + abundance outside the water snowsurface in network W. Therefore, the derived expression for theHCO + abundance can also be used in this disk region in chemicalnetwork W. Appendix B.2: Initial conditions
The e ff ects of the initial conditions on the abundance predictedby chemical network W were discussed in Section 3.1.2 andexpected to be of little importance for HCO + ’s ability to tracethe water snowline. Here, the e ff ects on the corresponding ra-dial emission profiles of the J = − + andH CO + are discussed and shown in Fig. B.1.Previously, it was derived that the column density of HCO + scales with the square root of the cosmic ray ionization rate (seeEq. 9). Similarly, the H CO + emission scales with the squareroot of the cosmic ray ionization rate because it is optically thin,see Fig. B.1. On the other hand, the HCO + emission does notscale with the cosmic ray ionization rate as it is optically thick.As the column density of HCO + decreases, the HCO + emissionbecomes less optically thick and approaches the scaling for thecolumn density.In Section 3.1.2, it was found that the HCO + abundance orcolumn density does not depend strongly on the initial abun-dance of CO or H O. This is is also seen in the radial emissionprofiles both for HCO + and H CO + , because the drop in theHCO + emission in the center is dominated by the e ff ect of theoptical depth of the continuum emission and molecular excita-tion. A small dependence of the HCO + emission on the initialabundance of gas-phase water is seen, but the abundance of gas-phase water seems to be low in the outer regions of protoplane-tary disks (Bergin et al. 2010; Du et al. 2017; Notsu et al. 2019;Harsono et al. 2020). Moreover, the abundance of gas-phase wa-ter in the inner disk can be as high as 10 − (Bosman et al. 2018). Appendix C: HCO + and H CO + J = − , J = − and J = − transitions Radial emission profiles for the J = − J = − J = − + and H CO + are presented in Fig. C.1. Appendix D: V883 Ori
An overview of the model parameters used for the representative
DALI model for V883 Ori is given in Table D.1. Predictions for
Article number, page 14 of 19. Leemker et al.: Chemically tracing the water snowline in protoplanetary disks with HCO + Fig. B.1: HCO + column densities for di ff erent initial conditions in chemical network W (top row) and the corresponding radialemission profiles for the J = − + (middle row) and H CO + (bottom row). The left-hand column shows theresults for a cosmic ray ionization rate of 10 − s − , 10 − s − and 10 − s − . The middle column shows the corresponding models foran initial CO abundance of 10 − , 10 − and 10 − . The right-hand column shows the corresponding models for an initial abundanceof gas-phase water of 3 × − , 3 × − and 3 . × − . The fiducial model is indicated with the light blue line in each panel anduses an initial abundance of 3 . × − for gas-phase water, 10 − for gas-phase CO and 10 − s − for the cosmic ray ionisation rate.The water snowline is indicated with a dashed black line and the position of the star is indicated by the symbol of a black star. Theradial emission profiles are convolved with a 0 (cid:48)(cid:48) .
05 beam.the corresponding emission of the J = − + and H CO + are shown in Fig. D.1 and Fig. D.2. Article number, page 15 of 19 & A proofs: manuscript no. Leemker_et_al_accepted
Fig. C.1: Same as the middle and bottom panel of Fig. 5, but then for the J = − J = − J = − + (top) and H CO + (bottom).Table B.1: Initial conditions for chemical network NW (no wa-ter) and W (water).Model parameter Initial valueNW Initial valueW Ref. ζ c . r . −
17) 1( − T gas DALI model
DALI model T dust DALI model
DALI model n (H ) DALI model
DALI model x (H + ) 0 0 x (H + ) 0 0 x (H) 0 0 x (e − ) 0 0 x (CO) 1( −
4) 1( − x (HCO + ) 0 0 x (H O + ) 0 0 x (H O) 0 3 . −
7) (a) x (H O(ice)) 0 0
Notes. a ( b ) represents a × b . Abundances are defined with respect tomolecular hydrogen. References: (a) McElroy et al. 2013.Article number, page 16 of 19. Leemker et al.: Chemically tracing the water snowline in protoplanetary disks with HCO + Fig. D.1: Same as Fig. 11 but then for the J = − + . The 0 (cid:48)(cid:48) . Article number, page 17 of 19 & A proofs: manuscript no. Leemker_et_al_accepted
Fig. D.2: Same as Fig. 11 but then for the J = − CO + . The 0 (cid:48)(cid:48) . Article number, page 18 of 19. Leemker et al.: Chemically tracing the water snowline in protoplanetary disks with HCO + Table D.1:
DALI model parameters for the representative modelfor V883 Ori.Model parameter Value Ref.
Physical structureR subl R c
75 AU Σ c
35 g cm − M disk .
25 M (cid:12) γ h c ψ Dust properties χ f ls ∆ gas / dust Stellar spectrum (1)
Type Outbursting L (cid:63) + acc (cid:12) L X − T e ff T acc T X ζ c . r . − Stellar properties (2) ˙ M −
5) M (cid:12) yr − M (cid:63) . (cid:12) (a) R (cid:63) . (cid:12) Observational geometryi ◦ (a)P.A. 32 ◦ (a) d
400 pc (b)
Notes. a ( b ) represents a × b . (1) The stellar spectrum is obtained bythe sum of the accretion luminosity L acc and an artificially high stellarluminosity L (cid:63) to shift the snowline to 47, 76 and 119 AU. (2) R (cid:63) is chosento obtain an accretion luminosity of 4 × L (cid:12)(cid:12)