© N. M. Murillo et al. 2022

1 Introduction

The physical structure of a cloud core determines the type of protostellar system that forms and how it evolves. The chemical structure within a cloud core, while an excellent tool to probe the physical structure, is also relevant in star and planet formation. The distribution and phase change of molecules as well as the resulting elemental composition of gas and ice decide the chemical complexity of the planets, comets, and meteorites that may form as part of the protostellar system. The region where a molecule undergoes a phase change to and from gas and ice is referred to as its snowline. Hence understanding how the physical structure determines the distribution of chemical content, including the snowline location, is relevant to understanding both the dynamical processes that occur during star formation as well as the chemical composition of planets, comets, and meteorites.

Within a star-forming cloud core, the protostar is the main source of luminosity and heat due to the release of gravitational energy from contraction and material accretion. The amount of protostellar heating dictates the temperature structure of the cloud core. In an idealized spherical scenario, the temperature alone would dictate the snowline locations within the cloud core. However, star-forming cloud cores are not spherically symmetric. The outflow cavity, flattened structures around the protostar (e.g., pseudo-disks and rotationally supported disks; hereafter referred to as disks for simplicity), and variations within the envelope density can all impact how heat is distributed within the cloud core. Studies have shown that heating mainly escapes through the outflow cavity in deeply embedded sources (van Kempen et al. 2009; Yıldız et al. 2015; Murillo et al. 2018a). Thus the extent of chemical richness in the outflow cavity provides insight into the luminosity of the protostar and the physical conditions of the envelope (e.g., Drozdovskaya et al. 2015; Murillo et al. 2018b; Tychoniec et al. 2019, 2020). Observations of embedded protostars (so-called Class 0 and I systems) have shown the presence of disks, both as flattened dust continuum structures (e.g., J∅rgensen et al. 2009; Enoch et al. 2011; Persson et al. 2016; Segura-Cox et al. 2018; Tobin et al. 2020) and rotationally supported disks traced in molecular gas (e.g., Murillo et al. 2013; Harsono et al. 2014; Yen et al. 2015, 2017; Maret et al. 2020). These disks show a wide range of geometries and radii ranging from a few 10 AU up to ~200 AU. Additional studies have shown that the presence of a disk can alter the temperature profile along the equator (disk mid-plane) of the cloud core (e.g., Murillo et al. 2015, 2018b; van’t Hoff et al. 2018b; Hsieh et al. 2019a). Multiplicity, that is two or more protostars within a single cloud core, can produce further asymmetries due to differences in luminosity from the multiple components and their locations with respect to each other (e.g., Chen et al. 2009; Koumpia et al. 2016; Murillo et al. 2016, 2018b). At an early evolutionary stage, protostellar luminosity is dominated by accretion, that is accretion luminosity (Hartmann & Kenyon 1996). Variability in protostellar luminosity has been detected toward several targets (V1647 Ori: Ábrahám et al. 2004; Andrews et al. 2004; Acosta-Pulido et al. 2007; Fedele et al. 2007; Aspin & Reipurth 2009; OO Serpentis: Kóspál et al. 2007; CTF93 216-2 Caratti o Garatti et al. 2011; VSX J205126.1: Covey et al. 2011; Kóspál et al. 2011; HOPS383 Safron et al. 2015; S255IR-SMA1 Safron et al. 2015; Liu et al. 2018; EC53 Herczeg et al. 2017; Yoo et al. 2017). Such variability is considered to be a product of the nonuniform accretion of material with a variable amount and frequency onto the protostar, that is to say episodic accretion (Audard et al. 2014).

This variability can then change the chemical structure of the cloud core, leading to a dynamic chemical evolution (Taquet et al. 2016; Molyarova et al. 2018). Molecular snowlines provide a way to characterize one aspect of the chemistry within the cloud core by indicating where a particular molecule undergoes a phase change (from gas to ice and vice versa). Hence, the snowline of a molecular species is defined as the radius at which the species is half frozen onto the dust grains (Hayashi 1981; van’t Hoff et al. 2017). Snowline locations are expected to shift back and forth during the star formation process due to the variable accretion luminosity of protostars (e.g., Hsieh et al. 2019b). Thus, measuring snowline locations can provide insight into protostellar evolution during the early embedded phase. In addition, given the likelihood of planets forming during the early stages of star formation (e.g., Tychoniec et al. 2020), these shifts likely also impact the formation and composition of planets. Whether planets form inside or outside of particular snowlines, such as the water snowline, influences the atmospheric C/O ratios and core compositions of planets (e.g., Öberg et al. 2011; Walsh et al. 2015; Eistrup et al. 2016; Bosman et al. 2019).

Observational line emission images of protostars have contributions from different processes (infall, outflow, rotation, internal, and external heating) and structures (envelope, outflow cavity wall, outflow, jet, and disk), in addition to geometrical effects imposed by the viewing angle. It is then difficult to determine what factors are producing the observed chemical structure. Consequently, deducing the snowline location of different molecular species can be challenging. Hence, simple physico-chemical models that explore and compare the impact of individual parameters can provide insight as to what factors affect the observed chemical structure. In addition, because observed chemical signatures may reflect the thermal history of the protostellar system, the comparison of models and observations can reveal the occurrences of different processes. For example, some chemical signatures are believed to trace a previous outburst of mass accretion (Lee 2007; J∅rgensen et al. 2013; Taquet et al. 2016; Molyarova et al. 2018; Wiebe et al. 2019).

Several models studying the location of snowlines have been previously published. Some of them focus on the snowline location of specific molecular species (Rab et al. 2017; van’t Hoff et al. 2017; Frimann et al. 2017), model a particular physical structure (Harsono et al. 2015; Bjerkeli et al. 2016; Rab et al. 2017), have limited chemical modeling (Harsono et al. 2015; Frimann et al. 2017; Owen 2020), or consider more molecular species but use a spherically symmetric envelope without the addition of a disk or outflow cavity (Visser et al. 2015). The effects of gas dispersal and dust evolution on the CO snowline location has been studied in disks around Herbig stars (Panic & Min 2017). The effect of dust grain sizes on molecular gas distributions and the H₂O snowline location has been modeled for a T Tauri disk (Gavino et al. 2021). These models all provide important insight into the factors that dictate the chemical structure set by snowlines, but they are unable to more generally constrain the impact of outflow cavities and disk-like structures on the location of snowlines within the cloud core.

This paper explores the physical conditions that affect the location of CO and H₂O snowlines, as well as quantifies the resulting abundance distribution and emergent line emission from their respective tracers, N₂H⁺ and HCO⁺. Since CO destroys N₂H⁺ via an ion-molecule reaction (N₂H⁺ + CO → HCO⁺ + N₂), N₂H⁺ increases in abundance when CO freezes out onto the dust grains, thus tracing the CO snowline location. A similar process occurs between H₂O and HCO⁺ (HCO⁺ + H₂O → H₃O⁺ + CO). To achieve this, a 2D (cylindrical symmetric) physical model was coupled with a reduced chemical network that includes key formation and destruction for each species. With cylindrical symmetry, the disk, outflow cavity, and envelope structures were included. The conditions of the physical model were varied within a range of parameters in order to study the effects of different physical structures within the cloud core. The range of parameters was chosen to reflect the variety of observed conditions in protostars. Section 2 describes the physical model, chemical network, and how simulated emission line maps were generated. Section 3 describes the results from the molecular distribution models and simulated emission maps, highlighting which parameters affect the snowline location and observability. The discussion and conclusions are presented in Sects. 4 and 5, respectively.

2 Methods

2.1 Physical models

Our physical models consist of a disk, outflow cavity, envelope, and a central heating source (also referred to as the protostar). The structure is illustrated in Fig. 1. The density profile is built starting from a rotating flattened envelope introduced by Ulrich (1976) $${\rho }_{\text{env}}\left(r,\mu \right)={\rho }_{\text{env}}{}_{,0}{\left(\frac{r}{r_{\text{cen}}}\right)}^{-1.5}{\left(1+\frac{\mu }{{\mu }_0}\right)}^{-1/2}{\left(\frac{\mu }{2{\mu }_0}+\frac{r_{\text{cen}}}{r}{\mu }_0^2\right)}^{-1},$$(1)

where r and θ are the spherical coordinates, ρ_env,0 is the density parameter used to scale the density profiles, r_cen is the centrifugal radius, and µ = cos(θ). The parameter µ₀ = cos(θ₀) is a solution of $${\mu }_0^3+{\mu }_0\left(r/r_{\text{cen}}-1\right)-\mu \left(r/r_{\text{cen}}-1\right)=0$$(2)

and sets the location and direction of streamlines. The centrifugal radius r_cen is arbitrarily set to be 50 au in the models presented in this work, representative of embedded small disks.

We adopted a density profile of a flared disk (Chiang & Goldreich 1997; Hartmann et al. 1998) as $${\rho }_{\text{disk}}\left(R,z\right)=\frac{\sum \left(R\right)}{\sqrt{2\pi }H\left(R\right)}\exp \left[-\frac{1}{2}{\left(\frac{z}{H\left(R\right)}\right)}^2\right]$$(3)

with Σ(R) defined as $$\sum \left(R\right)={\sum }_0{\left(\frac{R}{R_{\text{disk}}}\right)}^{-1}\exp \left[-\frac{R}{R_{\text{disk}}}\right]$$(4)

where R and z are the cylindrical coordinates, R_disk is the disk radius in au, and Σ₀ is the disk surface density at R_disk. The disk scale height H(R) is defined as a function of radius R $$H\left(R\right)=H_0{\left(\frac{R}{R_{\text{disk}}}\right)}^{1.3},$$(5)

(Chiang & Goldreich 1997), for which H₀ is the pressure scale height at the radius R_disk, and is tuned as a ratio normalized to R_disk.

The final gas density at a specific pixel is taken as the larger value between ρ_env and ρ_disk. To convert the H₂ density (cm⁻³) to dust density (g cm⁻³), we assumed a gas-to-dust mass ratio of 100 and a mean molecular weight for the gas of 2.3. It is important to note that the dust evolution is not taken into account, and we assume that the dust and gas are well mixed. The dust opacity including absorption and scattering is obtained from the DIANA Opacity Tool¹ (Woitke et al. 2016). This tool computes fast models given a grain size distribution. A power-law index of 3.5 (dn/da ∝ a^−3.5, where n is the number of grains with radius a) is used for both envelope and disk grain size distribution, with a maximum grain size of 1 µm, and 100 µm for the envelope, and disk, respectively. So whilst the full effects of dust evolution are not being taken into account, the assumption that dust grains in the disk are on average bigger than those in the envelope is adopted.

An outflow cavity is carved out from the envelope defined by $$z=z_0\times {\left(\frac{R}{z_0\tan \left({\text{Θ}}_{\text{outflow}}/2\right)}\right)}^{1.5}$$(6)

where Θ_outflow is the outflow cavity opening angle. We adopt z₀ ≡ 50 au.

Two-dimensional density and temperature profiles are generated with the above described structures. The dust temperature is calculated with RADMC3D² (Dullemond et al. 2012), which uses the Monte Carlo method, for a given dust opacity and density with a central heating source. An effective temperature T_eff and central stellar luminosity L_cen are set for the central heating source. In practice, the stellar radius is tuned to the luminosity using the Stefan-Boltzmann law with a given T_eff.

Fig. 1 Cartoon showing the structure of the physical models, along with the illustrated definitions of the disk radius Rdisk, disk scale height H0 at Rdisk, and outflow cavity opening angle Θoutflow.

2.1.1 Model grid setup

A grid of models is generated in order to study the effect of physical conditions on the locations of the CO and H₂O snowlines within the cloud core of embedded protostellar sources. The parameters which characterize the protostar, outflow, disk, and envelope are treated as free parameters. For each model one parameter is changed with the purpose to recognize the factors that influence the location of snowlines. Table 1 lists the parameters along with the range used in the grid of models.

The luminosity of the central protostar L_cen (L_⊙) is altered given the observed luminosity range of protostars and the variability produced by accretion processes. In addition, the effective temperature of the star T_eff (K) is also kept as a free parameter, and can reflect the mass of the protostar, and its evolutionary stage. For example, T_eff = 1500 K is suitable for a proto-brown dwarf, whereas T_eff = 5000 K is more representative of a solar-type protostar. As the uv-visible light from the star is absorbed and reradiated back out by the circumstellar dust at a longer wavelength, a higher effective temperature for the same luminosity can heat up the disk-envelope more efficiently.

Photons from the protostar escape mainly through the outflow cavity, heating the gas and kickstarting chemical processes within the cavity walls (e.g., Drozdovskaya et al. 2015; Yıldız et al. 2015; van Kempen et al. 2009; Murillo et al. 2018b,a). It is of interest to test whether the outflow cavity opening angle has an effect on the location of snowlines within the cloud core. Observations suggest that the outflow cavity opening angle Θ_outflow (degrees) increases as the protostar evolves (Arce & Sargent 2006; Velusamy et al. 2014; Hsieh et al. 2017; Heimsoth et al. 2022), while a recent study finds no evidence of outflow cavity opening angle growth during the Class I phase (Habel et al. 2021). However, these studies suffer from a lack of good evolutionary indicators (Whitney et al. 2003; Offner et al. 2009), and the outflow opening angle can be different depending on the tracers used to make the measurement. The outflow opening angle is defined in Eq. (6). It should be noted, however, that the models presented here only include simple molecular species. Modeling of complex organic molecules along the outflow cavity walls has been previously studied using more comprehensive gas-grain chemistry in, for example, Drozdovskaya et al. (2015).

As material is accreted from the cloud core onto the protostar (and the disk), the envelope density ρ_env decreases with time, thus changing this parameter emulates the evolution of the cloud core. The envelope density ρ_env is defined in Eq. (1).

Studies in the last decade have revealed both the presence of and lack of disks at all stages of protostellar evolution, with a wide range of physical characteristics (Reviews: Williams & Cieza 2011; Belloche 2013; Li et al. 2014; Surveys: Yen et al. 2015, 2017; Testi et al. 2016; Maury et al. 2019; Tobin et al. 2020; Individual sources: Murillo et al. 2013; Aso et al. 2017; Lee et al. 2017; Hsieh et al. 2019a). Hence, in our physical model, the outer disk radii R_disk (AU), total disk mass M_disk (M_⊙), and disk scale height H₀ are treated as free parameters. The disk density is constrained by these three parameters. Thus, for two disks with the same mass, one with a large radius and scale height will be less dense than one with a small radius and scale height. Models without disks are also generated in order to further explore the effect of these disk structures on the location of snowlines. In the case of no disk, the density structure is simply described by the flattened envelope profile (see Eq. (1)).

2.1.2 Fiducial models

As a point of reference, two fiducial models are generated: one model without a disk, and one with a disk. The density and temperature distributions of the fiducial models are shown in Fig. 2, and all parameters for the fiducial models are listed in Table 2.

Both fiducial models have an effective stellar temperature of T_eff = 5000 K, and envelope density of ρ_env = 10⁶ cm⁻³, representative of solar type stars. A relatively high luminosity of L_cen = 10 L_⊙ is chosen to better highlight features for comparison. Although not all protostellar sources exhibit such high luminosities, some embedded sources have been observed to have bolometric luminosities on the order of 10 L_⊙ (e.g., Murillo et al. 2016). The outflow cavity opening angle is set to Θ_outflow = 50° based on observations of embedded protostars (Arce & Sargent 2006).

Observations have shown the presence of large (≥ 100 AU) disk(-like) structures around protostars (e.g., Murillo et al. 2013; Harsono et al. 2014; Persson et al. 2016; Tobin et al. 2018), while further surveys show that smaller disk(-like) structures are common (e.g., Harsono et al. 2014; Yen et al. 2015, 2017; Maury et al. 2019; Maret et al. 2020). The fiducial model with a disk is then chosen to have a disk with R_disk = 150 AU, M_disk = 0.05 M_⊙, and H₀ = 0.05. The disk mass is chosen based on an assumed central protostar mass of 0.5 M_⊙. This is reasonable assuming the central protostar is a solar type star (T_eff = 5000 K) that has accreted at least half its mass in the embedded phase.

2.2 Chemical network

Given the grid size of each physical model, and the total number of physical models, for speed and simplicity, we use a reduced chemical network based on that compiled by the UMIST Database for Astrochemistry (Woodall et al. 2007; McElroy et al. 2013). Our reduced network captures the main gas-phase reactions important for the formation of HCO⁺ and N₂H⁺, and associated deuterated ions such as H₂D⁺, DCO⁺, and N₂D⁺ (see Table 3). As we are interested in HCO⁺ as an accessible tracer of the H₂O snowline location, and N₂H⁺ and DCO⁺ as tracers of the CO snowline position, we only consider singly deuterated forms, and neglect spin-state chemistry for simplicity. We allow the freeze-out of molecules onto dust grain surfaces as well as thermal desorption, and photodesorption (induced by both cosmic rays and stellar photons). We assume that H₂O, CO and N₂ have already formed. It is assumed that water is already present in the cloud core, in either ice or gas forms (van Dishoeck et al. 2013,?, and references therein). Additional reactions include the destruction of water by reaction with HCO⁺, and the formation of water through the dissociative recombination of H₃O⁺. These two reactions are used to model the water abundance in a simple way, with the key reaction destroying gas-phase water being H₂O + HCO⁺. The dissociative recombination of H₃O⁺ is included to reform the gas-phase water and is necessary to close the network. We do not include grain-surface chemistry, but we do allow the recombination of gas-phase cations with negatively charged grains. We note that this approach will not accurately model the chemistry occurring in the outflow cavity walls; however, we are interested in the locations of snowlines, which are present only in well-shielded gas. The reactions included in the network, along with the parameters to calculate the rate coefficients of each reaction are listed in Table 3. Our treatment of the chemistry represents an intermediate choice between a parametric model (e.g., Yıldız et al. 2010) and a full chemical model (e.g., Drozdovskaya et al. 2015; Notsu et al. 2021). Future work will explore the impact of a more complex chemical network on the abundances and distribution of key snowline tracers.

Several parameters of the chemical network can be changed in order to determine if a resulting outcome is due to a chemical or physical effect. These parameters include the initial abundances of the species in the network, and the binding energies E_b of CO, N₂, and H₂O (Table 1). Initial abundances are set for CO, N₂, H₂O, HD, and neutral grains, and naturally, for H₂ and H. The total initial abundances of CO, N₂, H₂O, HD and neutral grains relative to H₂ are set to 2 × 10⁻⁴, 1 × 10⁻⁴, 2 × 10⁻⁴, 1.6 × 10⁻⁵, and 1.3 × 10⁻¹², respectively. The initial abundances for CO, N₂ and H₂O are set to all ice. The binding energy value reflects the conditions under which a molecule is bound to grain surfaces. The binding energy of N₂ is kept constant for all models, and is set at $E_{\text{b},{\text{N}}_2}=955$ K (Collings et al. 2004; Garrod & Herbst 2006). To probe molecular binding under different ice environments and determine how binding energy affects the CO and water snowline locations, two binding energies are used for each molecule: E_b,_CO = 1150 K and 1307 K for CO (Collings et al. 2003, 2004; Bisschop et al. 2006; Noble et al. 2012); $E_{\text{b},{\text{H}}_2\text{O}}=4820$ K and 5700 K for water (Sandford & Allamandola 1993; Fraser et al. 2001). To simulate the absorption and desorption processes, we use an average dust grain size of 0.1 µm. Thus, absorption and desorption rates will be determined by the temperature and gas density (with a fixed gas-to-dust ratio) at a given binding energy (e.g., Anderl et al. 2016).

The chemical network is evolved up to 10⁸ years to ensure that steady state has been reached to produce molecular distributions for each set of conditions. Given that timescales for protostellar evolution and disk formation are shorter than this, we checked the validity of this assumption by comparing the results at 10⁸ years with those at earlier times. We can confirm that steady state in our reduced network is already reached by 1 Myr timescales. Because of this, there is no difference if the initial abundances are set to all ice, all gas, or a mix. While the network has the capability of time-dependent chemical models, in this paper the steady state output is used since the aim is to study the effect of cloud core physical structures on the locations of the CO and water snowlines. Additional time-dependent models to study the effect of episodic accretion bursts will be considered in a subsequent paper.

Fig. 2 Density (top row) and temperature (bottom row) distribution of the two fiducial models: without a disk (left) and with a disk (center and right). The right column shows a zoom in of the fiducial model with the disk marked with an orange box in the center column.

2.3 Simulated line emission maps

In order to examine whether the effects of physical and chemical structure on the CO and H₂O snowline locations are observable, simulated emission maps are generated from the molecular distributions. The molecules N₂H⁺ and HCO⁺ are typically used in observations to determine the CO and H₂O snowline positions, respectively (van’t Hoff et al. 2017; Frimann et al. 2017; Hsieh et al. 2018, 2019b). Simulated emission maps of CO (three isotopologues: ¹²CO, ¹³CO, and C¹⁸O), N₂H⁺, DCO⁺, H₂O, and HCO⁺ for the two fiducial models are produced in order to examine the robustness of the emission from these species in tracing snowline locations. Simulated emission maps of N₂H⁺ and HCO⁺ for all the models discussed in Sect. 3 are produced. These are used to determine which effects have an observable, and distinguishable, impact on the observation of snowline locations.

The simulated line emission channel maps are calculated using RADMC3D under the assumption of local thermal equilibrium (LTE). Seven different inclination angles are used: 0° (face-on), 15°, 25°, 45°, 65°, 75°, and 90° (edge-on). The distance to Perseus of 293 pc (Ortiz-León et al. 2018; Zucker et al. 2018) is adopted for ray-tracing and a map size of 60″ × 60″ (17 580 × 17 580 AU) is used for N₂H⁺, while a smaller map size of 6″ × 6″ (1758 × 1758 AU) is used for HCO⁺. Table 4 lists the transitions modeled. These transitions are selected because they are commonly observed, and observations of these transitions toward protostars have been reported in literature.

Observational parameters such as bolometric luminosity change with the inclination of the cloud core relative to the line of sight (e.g., Whitney et al. 2003; Crapsi et al. 2008). Thus, the spectral energy distribution (SED) of each model is also generated with RADMC3D for each inclination angle. The modeled SEDs are used to calculate the bolometric luminosity after radiative transfer.

To calculate the emergent line emission maps, a velocity field is required. The velocity structures from infalling-rotating cores (Ulrich 1976; Terebey et al. 1984) are used, and defined as $$V_r=-{\left(\frac{GM}{r}\right)}^{1/2}{\left(1+\frac{\cos \theta }{\cos {\theta }_0}\right)}^{1/2},$$(7) $$V_{\theta }={\left(\frac{GM}{r}\right)}^{1/2}\left(\cos {\theta }_0-\cos {\theta }_0\right){\left(\frac{\cos \theta +\cos \theta }{\cos {\theta }_0{\sin }^2\theta }\right)}^{1/2},$$(8) $$V_{\varphi }={\left(\frac{GM}{r}\right)}^{1/2}{\left(1-\frac{\cos {\theta }_0}{\cos {\theta }_0}\right)}^{1/2}\frac{\sin {\theta }_0}{\sin \theta }$$(9)

where G is the gravitational constant, M is the central source mass, r is the radius, θ is the angle in polar coordinates, and θ₀ is the angle of the initial velocity in polar coordinates.

Keplerian rotation is adopted for the inner disk region with R < r_cen. The gravitational force decreases as the orthogonal distance from the disk mid-plane increases such that (Pinte et al. 2018) $$V_{\varphi }={\left(\frac{GM}{r}\right)}^{1/2}\sin \left(\theta \right),r<r_{\text{cen}}.$$(10)

For the radial and meridional velocities, a dimensionless factor is applied to the above equations as a function of θ and r, ${\left(\frac{r}{r_{\text{cen}}}\right)}^{2\theta /\pi }$. This factor gradually decreases the rotation velocity along the mid-plane, and, with a fixed r_cen, the velocity field is scaled by the squared root of the central mass M that is set to be 0.5 M_⊙.

The modeled transitions and convolving beams from the Atacama Large Millimeter/submillimeter Array (ALMA) observations of Hsieh et al. (2019b) are adopted. For N₂H⁺ 1−0 (ALMA Band 3), the convolving beam is 2″, while for HCO⁺ 3−2 (ALMA Band 6) the convolving beam is 0.25″. For water, two transitions with low upper energy Eup were used, one at 183 GHz with E_up = 205 K (ALMA Band 5), and a second at 557 GHz with E_up = 61K (Herschel Space Observatory). The latter transition has been observed and reported in Kristensen et al. (2012) for embedded Class 0 and I low-mass protostars. Continuum subtraction is not performed on the simulated line emission channel maps.

3 Results

The outcome of the models is described in this section. First, the robustness of each snowline location tracer is discussed based on the two fiducial models. Then each parameter explored in the models and their effect on each snowline position will be described. In this way which parameters produce an observable effect can be examined. Table 5 lists the parameters studied and which snowline location is affected by a change in a parameter. Details on how the simulated line emission maps are generated are presented in Sect. 2.3. The effect of inclination and spatial resolution on the observed CO and H₂O snowline locations will then be considered. Finally, plots of peak radius of N₂H⁺ and HCO⁺ versus luminosity for the explored parameter space are provided and described.

A sample of the models from the full grid is discussed in Sects. 3.1 and 3.2, as well as in Appendix A. The full grid of molecular distributions and simulated line emission maps are available online³.

3.1 CO snowline location

The two fiducial models are used to explore the robustness of N₂H⁺ and DCO⁺ as tracers of the CO snowline location. Figure 3 shows the integrated intensity (moment 0) maps of the simulated line emission, and slices extracted along z = 0. The slices are normalized to their respective peaks for comparison purposes, since CO is several orders of magnitude brighter than N₂H⁺ and DCO⁺. From the z = 0 slices, N₂H⁺ begins to increase at a radius where CO is between 60% and 70% of its maximum brightness. The N₂H⁺ emission then peaks where CO is at ~1% of its maximum brightness. At low envelope densities (e.g., 10⁵ cm⁻³), N₂H⁺ is no longer a robust tracer of the CO snowline location (see Sect. 3.4). In contrast, the DCO⁺ emission is located between the CO and N₂H⁺ peaks. From the z = 0 slice the DCO⁺ emission peaks before CO has decreased to half its maximum brightness. This is clearly visible in the case without a disk. In the case with a disk, the DCO⁺ emission is less bright, but the trend is the same as evidenced by the DCO⁺ “wings” at ~20% brightness (Fig. 3). It should be noted that the central DCO⁺ peak in the case with a disk is a contribution from the continuum. Both behaviors, that of N₂H⁺ and DCO⁺ are expected from the formation and destruction pathways of each molecule (see Table 3, reactions 11 and 13). These results suggest that N₂H⁺ and DCO⁺ can provide observational constraints on the CO snowline location. It should be noted, however, that the limited chemical network used in this work does not consider additional direct or indirect reactions leading to the formation or destruction of N₂H⁺ and DCO⁺. For example, the warm formation pathway for DCO⁺ (e.g., Favre et al. 2015; Murillo et al. 2018b) is not included in our current network. This is further discussed in Sect. 4.2. Figure A.1 shows how the emission from the CO isotopologs, ¹³CO and C¹⁸O, does not coincide well with the snowline location traced by the N₂H⁺ emission.

A radial shift of the CO snowline location occurs throughout the envelope with changes in the binding energy E_b,CO, luminosity, and envelope density. These parameters change the CO snowline position whether a disk is present or not. The location of the CO snowline moves to smaller radii with an increase in the binding energy (E_b,co = 1150 1307 K, Fig. B.1). These binding energies result in sublimation temperatures of ~21K and ~25 K, respectively. The simulated emission line maps show that the snowline position shifts by at least a few 100 AU when changing the binding energy, whether there is a disk present or not. Hence, the effect of E_b,CO on the CO snowline location is potentially observable and measurable, but needs to be differentiated from other effects. With increasing luminosity L_cen, the 20 K gas temperature contour moves to larger radii, and consequently, so does the CO snowline position. This is expected, and is widely used as a test for occurrence of episodic accretion in embedded protostellar systems. The gas density of the cloud core plays a major role in the distribution of different molecules. Thus, at low densities (ρ_env < 10⁶ cm⁻³) the CO snowline location moves outward to larger radii, whereas at higher densities (ρ_env > 10⁶ cm⁻³) the snowline location moves inward to smaller radii (Fig. B.3). In either case, the snowline position shift does not coincide with the change in temperature structure. This is due to the balance between the thermal desorption and freeze-out rates. The thermal desorption rate varies with density as n^b where b is a value between 0 and 1, whereas the rate of freeze-out goes as n² (Cuppen et al. 2017). Thus, at low densities, the freezeout of molecules is slow causing snowline locations to occur at lower temperatures, and vice versa. Consequently, the snowline position does not occur at a single, well defined temperature as commonly assumed, but is also dependent upon density. In turn, density also affects the robustness of snowline tracers. The simulated line emission maps confirm that the change of cloud core density produces an observable and measurable effect on the snowline location.

The presence of a disk structure produces an inward shift in the CO snowline position along the disk mid-plane. From the molecular distributions (Fig. 4), it can be noted that “rings” of N₂H⁺ and DCO⁺ are generated around the disk edge, in contrast to the more spherical distributions of both molecules in the case without a disk. For N₂H⁺, the effect is also noticeable in the simulated emission maps. For DCO⁺ the effect is seen in the simulated emission maps at i = 0° (face-on), but less noticeable at i > 0° due to the contribution from the continuum component in the simulated emission maps. However, DCO⁺ rings have been observed toward embedded protostars (e.g., Murillo et al. 2018b) and disks (e.g., Mathews et al. 2013; Salinas et al. 2018). The effect on the CO snowline location from each parameter that defines the disk geometry is explored here. This is done by comparing the fiducial models with models where only one parameter is changed, either disk radius R_disk, disk mass M_disk, or scale height H₀. Figure 4 shows the two fiducial models in addition to one model for each disk parameter that is changed. The molecular distributions show that the 20 K gas temperature contour for the fiducial model without a disk is located at ~ 1400 AU (Fig. 4 first column). While a disk with R_disk = 50 and 150 AU cause the 20 K contour to move inward to 1000 and 800 AU along the disk mid-plane, respectively (Fig. 4 second and third columns). Changing M_disk by a factor of a few or even an order of magnitude only shifts the CO snowline location by ~ 100 AU or less (Fig. 4 third, fourth and sixth columns). In contrast to the molecular distributions, the simulated line emission maps do not show a particularly strong effect on the CO snowline position as traced by N₂H⁺ due to disk mass M_disk and radius Rdisk (Fig. 4 bottom rows). Interestingly, changing H₀ not only alters the CO snowline position along the disk mid-plane, but also vertically (Fig. 4 third and fifth columns). This effect is evident in the simulated line emission images, where a flared disk (e.g., H₀ = 0.3) efficiently shadows the envelope both along the disk midplane and vertically. Thus, a flared disk causes a larger radial and vertical shift inward of the CO snowline position (between 100 to 1000 AU) than that caused by a flatter disk. The effect of disk density on the CO snowline location is implicit in the models (Fig. 4 fifth and sixth columns). In the simulated line emission maps, disk density does not produce a significant snowline position shift, similar to disk mass and radius.

Two parameters in our models are found to have no effect on the CO snowline location: effective temperature of the protostar, and outflow cavity opening angle. The effective temperature of the central protostar generates a more noticeable effect along the outflow cavity than within the cloud core. Thus, only a slight increase of the presence of CO in the gas phase along the outflow cavity at z > 1500 AU is apparent in the models. The outflow cavity opening angle Θ_out only alters the CO gas distribution at r, z ≳ 1500 AU.

Fig. 3 Robustness of N2H+ and DCO+ as CO snowline location tracers. Integrated intensity maps of the simulated line emission show CO (top row), DCO+ (second row) and N2H+ (third row) emission. Contours in all three rows are for CO at 10, 30, 50, and 100 mJy beam−1 km s−1. Images that have been scaled for better comparison have the scaling factor noted on the top right corner. Different inclinations are shown in order to determine if inclination affects the robustness of N2H+ and DCO+ as snowline position tracers. Inclinations shown are from face-on (i = 0°) to edge-on (i = 90°). The first three columns show the fiducial model without disk, while the other three columns show the fiducial model with disk. Slices extracted along z = 0 from the simulated line emission maps are shown in the bottom row normalized to the peak of each profile.

Fig. 4 Effect of disk geometry parameters on the CO snowline location. Molecular distributions are shown in the top four rows. The fifth row shows the maps of the integrated intensity of the simulated line emission for N2H+ 1−0. The sixth row shows the corresponding slice extracted along z = 0 from the N2H+ simulated emission maps. The average distance of the peaks from the center (i.e., peak radius) is indicated in AU in the top left corner. These positions are shown with orange vertical lines. The first and third columns show the fiducial models without and with disk, respectively. The additional columns show the effect of changing disk radius (second column), disk mass (fourth column), scale height H0 (fifth column), and disk density (sixth column). Fractional abundances of all molecular species are relative to total number density of H2 (top row). The black and gray contours show gas temperature and density, respectively. The emission maps are shown at i = 45° and convolved to a beam of 2″.

3.2 H₂O snowline

The robustness of HCO⁺ emission as a water snowline location tracer is examined by using the two fiducial models (Fig. 5). Models of both o–H₂O and p–H₂O emission are examined. HCO⁺ emission begins to increase at the H₂O half maximum, and peaks 20–40 AU further out. Similar to N₂H⁺ for the CO snowline, the HCO⁺ emission peak is not at the exact location of the water snowline, but provides an indication of where water begins to freeze-out onto the dust grains. As noted in Sect. 3.1, the limited chemical network in this work does not account for all the destruction and formation pathways for HCO⁺ and H₂O, in addition to the simple treatment of water in the network.

The H₂O snowline location is altered by the binding energy $E_{\text{b},{\text{H}}_2\text{O}}$, luminosity, and envelope density. Increasing the binding energy of water ($E_{\text{b},{\text{H}}_2\text{O}}$ = 4820 → 5700 K) causes the H₂O snowline position to shift inward between 20 and 50AU (Fig. B.2). The simulated emission line maps show that the effect of binding energy on the H₂O snowline location is observable. With sufficient spatial resolution and additional constraints, the effect can potentially be measured. Increasing the protostellar luminosity moves the 100 K gas temperature contour outward to larger radii, in turn shifting the location of the H₂O snowline as well. This is expected, given how much more luminous protostellar sources have relatively larger warm (≳100 K) regions than less luminous sources. Changing the envelope density has the same effect on the H₂O snowline location as it does on the CO snowline position (Fig. B.4). In the simulated emission line maps without a disk, the H₂O snowline location shifts proportional to luminosity, and inversely proportional to envelope density ρ_env. Overall, the H₂O snowline location shifts on the order of ~100 AU or less.

Based on the molecular distributions, changing the disk radius R_disk, disk mass M_disk, and scale height H₀ has no apparent effect on the actual location of the H₂O snowline (Fig. 6). In a similar manner, the disk density does not seem to have a significant effect on the H₂O snowline location in the molecular distributions. However, the disk geometry does help to concentrate the H₂O gas extent between the outflow cavity wall and disk surface (Fig. 6 fifth and sixth columns). Hence, a very flared disk (H₀ = 0.3) causes the gas-phase water to be located in a narrow region between the outflow cavity wall and the disk surface. Similarly, changing the outflow cavity opening angle will affect the vertical distribution of water gas but not the snowline location. This effect becomes relevant in the face-on inclination for the case without a disk (see Sect. 3.3). The simulated line emission maps provide further insight into the H₂O snowline location when a disk is present (Fig. 6 bottom rows). The presence of a disk concentrates the gas-phase water by limiting snowline location from shifting outward. This is somewhat dependent on luminosity and density (Sect. 3.4), but occurs regardless of disk radius R_disk, mass M_disk, scale height H₀.

The effective temperature of the central source only has a slight effect on the presence of gas-phase water in the z direction along the outflow cavity walls. The simulated line emission images further support this result.

3.3 Inclination and spatial resolution

To properly compare the results presented in this work with observations, the inclination of the protostellar cloud core and the spatial resolution of the observations must be considered. The inclination i of a cloud core, the orientation with respect to the line of sight, can change the observed source SED (Fig. B.5), derived bolometric luminosity, and determine which features of the circumstellar material at scales < 100 AU can be observed (e.g., Crapsi et al. 2008). The spatial resolution of observations is relevant for characterizing snowline locations, in particular those located within a few hundred AU.

Figures 3 and 5 show how inclination affects the peak location of N₂H⁺ and HCO⁺ line emission, respectively, for the two fiducial models. In this work, i = 0° is defined as face-on, whereas i = 90° is edge-on. The peak position of simulated N₂H⁺ emission decreases with increasing inclination (Fig. B.6). This trend is more evident when a disk is present (see Sect. 3.4). For i = 0° to i = 90°, the peaks of N₂H⁺ emission shift less than 500 AU in the case without a disk. In the case with a disk, the N₂H⁺ emission peaks move inward from i = 0° to i = 90° by 25%. Hence, the inclination of the cloud core along the line of sight does not affect the robustness of N₂H⁺ and DCO⁺ as CO snowline location tracers.

The simulated HCO⁺ emission peak radius decreases with increasing inclination (Fig. B.6). In the case without a disk, the outflow cavity opening angle Θ_out causes the HCO⁺ peak position to be overestimated for low inclinations and in particular i = 0° (see Sect. 3.4). This is because Θ_out alters the vertical distribution of the HCO⁺ gas, and looking down along the outflow cavity does not allow the radial and vertical distributions to be distinguished. When a disk is present, the disk itself is the main constraining factor and the HCO⁺ emission peak position is not overestimated at low inclinations due to Θ_out or other factors. For the fiducial models, the peak moves from 210 (i = 0°) to 38 AU (i = 90°) for the case without a disk, and from 90 (i = 0°) to 9 AU (i = 75°) in the case with a disk (Fig. B.6). At 0° < i < 45°, the HCO⁺ gap can be more clearly seen than at higher inclinations. At i = 45°, the HCO⁺ gap is only marginally seen, and may be difficult to observe depending on the disk geometry when a disk is present. This is most likely caused by the optical depth, that is, column density, of the HCO⁺ emission (see also Hsieh et al. 2019b; van’t Hoff et al. 2021). Thus, the robustness of HCO⁺ as a snowline position tracer of H₂O is affected by inclination, more so when a disk is present. Consequently, the H₂O snowline location can be better traced with HCO⁺ emission at inclinations close to face-on.

To explore how spatial resolution affects the measurement of different snowlines positions, simulated line emission maps are convolved with three spatial resolutions (beams): 2″, 0.25″, and 0.1″ (Fig. 7). These spatial resolutions are representative of ALMA configurations C-1, C-5, and C-7 in Band 6, and C-2, C-6, and C-8 for Band 3, respectively. The distance for the simulated line emission maps is set to 293 pc (Perseus, Ortiz-León et al. 2018), hence the three beams represent physical scales of ~600, 70, and 30 AU, respectively. The peak of N₂H⁺ emission in the case without a disk does not shift with different spatial resolutions. In the case with a disk, the N₂H⁺ peak moves inward (<100 AU) with increasing beam size. Lower spatial resolutions are better suited to trace the CO snowline location with N₂H⁺ line emission. The reason is two-fold: the CO snowline position is usually located at radii ≳ 1000 AU, and higher spatial resolutions can recover the peak but not the extended emission which is necessary to properly characterize the CO snowline location from observations. Based on the physico-chemical models (Sect. 3), the H₂O snowline location will extend out to a few 100 AU at most. In addition, the HCO⁺ gap is only marginally seen at 15° < i ≤ 45°. Thus, a spatial resolution of 2” cannot resolve the H₂O snowline location with HCO⁺ line emission. A spatial resolution of 0.1″ would be able to unambiguously detect the HCO⁺ gap for a protostar with luminosity of 10 L_⊙ (the fiducial model) at the distance of Perseus. A beam of 0.25″ can also detect the HCO⁺ gap for L_cen ≥ 10L_⊙, but would benefit from high signal-to-noise observations. The HCO⁺ gap may not be observable for sources with lower luminosities unless the spatial resolution can trace few AU scales. In contrast, the HCO⁺ gap will be easily observable with lower spatial resolutions for much more luminous sources, or those that have undergone a strong luminosity burst in the recent past.

Fig. 5 Robustness of HCO+ as a H2O snowline location tracer. Simulated emission maps show HCO+ (top row), p–H2O (second row) and o–H2O (third row) emission. Contours in all three rows are H2O at 0.4, 0.6, 0.9, 1.6, 1.9, and 2.1 mJy beam−1 km s−1. Images that have been scaled for better comparison have the scaling factor on the top right corner. Different inclination are shown in order to determine if it affects the robustness of HCO+ as a snowline location tracer. Inclinations shown are from face-on (i = 0°) to edge-on (i = 90°). The first three columns show the fiducial model without disk, while the other three columns show the fiducial model with disk. Slices extracted along z = 0 from the simulated line emission maps are shown in the bottom row normalized to the peak of each profile.

Fig. 6 Effect of disk geometry parameters on the H2O snowline location. Molecular distributions are shown in the top three rows. The fourth row shows the simulated intensity integrated line emission maps for HCO+ 3−2. Images that have been scaled for better comparison have the scaling factor on the top right corner. The fifth row shows the corresponding slice extracted along z = 0 from the HCO+ simulated emission maps. The average distance of the peaks from the center (i.e., peak radius) is indicated in AU in the top left corner. These positions are shown with orange vertical lines. No scaling has been applied to the HCO+ profiles. The first and third columns show the fiducial models without and with disk, respectively. The additional columns show the effect of changing disk radius Rdisk (second column), disk mass (fourth column), scale height H0 (fifth column), and disk density (sixth column). Fractional abundances of all molecular species are relative to total number density of H2 (top row). The black and gray contours show gas temperature and density, respectively. The simulated line emission maps are shown at i = 45° and convolved to a beam of 0.25″.

3.4 N₂H⁺ and HCO⁺ peak radius versus luminosity

In order to provide a quick way to compare the models presented here with observations, plots of the N₂H⁺ and HCO⁺ line emission peak radius are provided in Figs. 8 to 14. Since inclination impacts the measured peak of N₂H⁺ and HCO⁺, the range of modeled inclinations are shown as shaded areas, with a curve showing the trend for i = 45°. Section 3.3 demonstrated the relevance of spatial resolution in measuring the emission peaks of N₂H⁺ and HCO⁺. Figure 9 shows the three spatial resolutions (beam sizes) versus luminosity. The three spatial resolutions are described in Sect. 3.3 and the simulated emission maps are shown in Fig. 7.

The luminosity in the plots is that of the central protostar (L_cen). Thus, Fig. 8 shows a comparison of L_cen with the bolometric luminosity L_bol derived from the simulated line emission maps. While there is a difference between L_cen and L_bol, the trends are similar. Figures 10 and 11 show the effects of binding energy and effective temperature, respectively, on the snowline location.

The effect of low envelope density on N₂H⁺ as a tracer of the CO snowline location can be better visualized when comparing peak radius versus luminosity for different envelope densities (Fig. 12). As noted in Sect. 3.1, the freeze-out of molecules at low densities is much slower, and thus the snowline locations are no longer clearly defined. The N₂H⁺ curve at ρ_env = 10⁵ cm⁻³ in Fig. 12 is not indicative of the CO snowline location. The curve drops out at luminosities >50 L_⊙ because the N₂H⁺ emission has become so diffuse and extended that a clear peak is no longer discernible. This is because the envelope density has become so low that the reaction between N₂H⁺ and CO (ID 13 in Table 3) is no longer a relevant destruction path for N₂H⁺. Hence N₂H⁺ is present in the gas phase in a larger region of the envelope and not correlated to the distribution of CO. The N₂H⁺ curve at ρ_env = 10⁵ cm⁻³ in Fig. 12 is kept as a reference, but caution must be taken when comparing this model with observations. This result illustrates that the snowline location of molecules also depends on density, and consequently so does the robustness of a molecular tracer for a particular snowline location.

An interesting result worth noting is how the presence of a disk limits the outward shift of the H₂O snowline position. This effect is somewhat dependent on luminosity and density. For luminosities <10 L_⊙, the H₂O snowline position increases sharply up to about 50 AU, regardless of which parameter is changed. For luminosities ≥ 10 L_⊙, the water snowline position is limited to within radii below 150 AU. The H₂O snowline position moves out to larger radii with decreasing disk density. Thus, for a disk with the same radius and scale height, a lower disk mass will allow the snowline position to move out further than a higher disk mass (Fig. B.3) In contrast, when a disk is present the CO snowline location shifts inward by ~1000 AU along the disk mid-plane relative to the models without a disk. However, the CO snowline location still shifts to larger radii with increasing luminosity when a disk is present.

Fig. 7 Simulated line emission maps showing the effect of spatial resolution (0.1″, 0.25″, and 2″) on the measured emission peak positions of N2H+ and HCO+. All models have an inclination of 45°. The first three columns show the fiducial model without disk. The other three columns show the fiducial model with disk. The second and third rows show corresponding slices extracted along the z = 0, with the number on the top left corner indicating the distance in AU of the peak from the center (i.e., peak radius). These positions are shown with orange vertical lines.

4 Discussion

4.1 Caveats

Several caveats to the models presented in this work should be noted. The chemical network only takes into account a limited set of simple molecules (three atoms or less per molecule) and a reduced number of reactions. The reduced chemical network in this work does not allow the additional chemical processes to be studied in the context of the CO and H₂O snowline locations. Examples of such processes are alternative formation and destruction paths for molecular species in the network (e.g., warm DCO⁺, Favre et al. 2015; Murillo et al. 2018b), extent of other simple molecules (e.g., SO, CN), and the relation between complex molecular species and snowline position. Since water is given a simple treatment in our network, additional effects on the H₂O snowline location from water chemistry cannot be characterized or discussed here. It is known from previous observational studies with Herschel (e.g., Kristensen et al. 2012) that gas-phase water is also found along the outflow cavity. We do not include processes such as photodissociation and photoionization in our models, and so the chemistry in the outflow cavity is not accurately modeled. Hence, the calculated abundances within the outflow cavity are not shown in the molecular distributions and simulated line emission maps. This could alter the distribution and measured peak radius of the HCO⁺ emission, and consequently the inferred water snowline location. It is also possible that wider cavities could lead to more of the envelope being heated, leading to more water in the envelope, and a different distribution of HCO⁺ emission. To accurately simulate the chemical composition of outflow cavity walls, photodissociation and hot gas-phase chemistry are needed (e.g., Drozdovskaya et al. 2015).

Opacity also plays a role in the simulated line emission maps and the measured emission peak location, especially for HCO⁺. In our models, optical depths are mainly determined by the density and kinematic structure. If the line emission is optically thick, the intensity map cannot properly reflect the molecular spatial distribution (e.g., HCO⁺, van’t Hoff et al. 2021). It is also worth noting that dust grain growth can occur in the early evolutionary stages and can block line emission from the region immediately around a protostar (e.g., Harsono et al. 2018). This could result in a ring structure after continuum subtraction (Lee et al. 2019). No continuum subtraction is performed on our simulated line emission images to avoid this artificial chemical substructure.

Fig. 8 Peak radius of N2H+ (left column), and HCO+ (right column) simulated emission versus luminosity for the fiducial models without disk (top row), and with disk (bottom row). The black hatched area shows the central protostellar luminosity, while the colored area shows the bolometric luminosity. The solid curve shows the relation for i = 45°. For N2H+, the shaded area indicates the range of peak radii for inclinations between 0° (face-on) and 90° (edge-on). For HCO+, the shaded area indicates the range of peak radii for inclinations between 15° and 90°, while the dotted line shows the peak radius versus luminosity for i = 0° (face-on).

4.2 Robustness of N₂H⁺ and HCO⁺ as snowline position tracers

Observational studies of CO and H₂O snowline positions typically use the peak emission of N₂H⁺ and HCO⁺ to trace the respective snowline locations (e.g., Qi et al. 2015, 2019; Hsieh et al. 2018, 2019b; van’t Hoff et al. 2018a). The simulated line emission maps discussed in Sect. 3 examined the robustness of N₂H⁺ and HCO⁺ as snowline position tracers, as well as DCO⁺, within the context of the limited chemical network used here. To fully assess whether N₂H⁺ and DCO⁺ are robust tracers of the CO snowline location, and HCO⁺ of the H₂O snowline location, a complete chemical network which includes all formation and destruction pathways for the relevant species is needed. The current work aims to find the conditions that produce significant and observable impacts on the CO and H₂O snowline locations. Thus, complexity is sacrificed for computational speed. The use of a more complete network informed by the results in this work will be investigated in the future.

Previous studies of the robustness of N₂H⁺ as a tracer of the CO snowline location focused on isolated protoplanetary disks (few 100 AU; van’t Hoff et al. 2017). This work examines the robustness of N₂H⁺ in the envelope (few 1000 AU) during the early embedded protostellar phase. The robustness of HCO⁺ was studied observationally in van’t Hoff et al. (2018a), and was argued to be a good tracer of the H₂O snowline position. However, the optical depth of the HCO⁺ emission can affect the real peak radius related to the snowline location (van’t Hoff et al. 2020). In our results, we find that by changing the optical depth (i.e., ρ_env), the HCO⁺ emission might suggest that the snowline is located at a different latitude if there is no disk present. We note that this process is also affected by the velocity structure and inclination which determines the optical depth at a specific velocity range. A possible solution might be to observe a range of isotopologues with different abundances that are less optically thick (e.g., van’t Hoff et al. 2021).

Although the peak of the N₂H⁺ emission does not trace the exact position of the CO snowline, it does provide an accessible observational proxy, except at low envelope densities where N₂H⁺ no longer traces the CO snowline location (see Sect. 3.4 and Fig. 12). van’t Hoff et al. (2017) also found that N₂H⁺ peaks further out than the CO snowline position in protoplanetary disks, which is consistent with the results presented here. While van’t Hoff et al. (2017) suggests chemical modeling is necessary to derive a robust location of the CO snowline using N₂H⁺ observations, the results could still be degenerate, and require more observational constraints. Constraining additional parameters, such as the envelope density, and the presence of a disk would prove to be helpful, but it would require multiwavelength observations with a range of spatial resolutions and sufficient spectral resolution to probe the kinematics. Additional molecular tracers would also provide constraints regarding the CO snowline location. By itself, DCO⁺ is not a robust tracer of the CO snowline location since it will also strongly depend on density and chemistry (e.g., Qi et al. 2015). DCO⁺ also has a second, warm gas formation pathway (Favre et al. 2015) not included in the chemical network used here. This warm formation pathway can enhance DCO⁺ in the disk region, and can be easily observed with the DCO⁺ 5–1 transition (Murillo et al. 2018b). However, combining N₂H⁺ and DCO⁺ observations can provide a better estimate of the CO snowline position (Fig. 3). Furthermore, having a second CO snowline location tracer can prove useful in protostellar systems where the outer envelope is so cold that N₂H⁺ (or N₂D⁺) is significantly less abundant and difficult to detect because of the depletion of gas-phase N₂ onto dust grains (e.g., VLA 1623-2417: Murillo et al. 2015).

Similar to N₂H⁺, HCO⁺ emission peaks some distance away from the H₂O snowline location. However, HCO⁺ can still be a good observational proxy and provide an outer limit to the H₂O snowline location, as suggested by observational studies (Hsieh et al. 2019b; van’t Hoff et al. 2018a), as long as HCO⁺ emission can be spatially resolved. Given that H₂O is in the gas phase, where the gas temperature is ≥100 K, other molecular species that sublimate at similar temperatures, such as methanol, could also be used as independent tracers of the approximate location of the H₂O snowline. However, similar to HCO⁺, their successful detection could be dependent on source inclination, spatial resolution, and presence of a disk.

Fig. 9 Peak radius of N2H+ (left column), and HCO+ (right column) simulated emission versus central protostellar luminosity for the fiducial models (black hatch) without disk (top row), and with disk (bottom row). The different shaded areas show how the measured peak radius is altered by the spatial resolution used for observations. The solid curve shows the relation for i = 45°. For N2H+, the shaded area indicates the range of peak radii for inclinations between 0° (face-on) and 90° (edge-on). For HCO+, the shaded area indicates the range of peak radii for inclinations between 15° and 90°, while the dotted line shows the peak radius versus luminosity for i = 0° (edge-on).

Fig. 10 Peak radius of N2H+ (left column), and HCO+ (right column) simulated emission versus central protostellar luminosity for the fiducial models (black hatch) without disk (top row), and with disk (bottom row). The shaded areas show the effect of CO and H2O binding energy on the peak radius of N2H+ and HCO+. The solid curve shows the relation for i = 45°. For N2H+, the shaded area indicates the range of peak radii for inclinations between 0° (face-on) and 90° (edge-on). For HCO+, the shaded area indicates the range of peak radii for inclinations between 15° and 90°, while the dotted line shows the peak radius versus luminosity for i = 0° (face-on).

4.3 Chemical effects

The binding energy of molecules onto dust grains depends on the conditions under which the molecule sticks to the grain. Binding energies of molecules in “pure” ices are typically lower than when the same species are in mixtures, for example when mixed with water ice (e.g., Martm-Doménech et al. 2014) This has an impact on the sublimation temperature of a particular molecule, which in turn can alter its snowline location. Both molecular distributions and simulated line emission maps reflect this behavior (Figs. B.1 and B.2), showing that the CO and H₂O snowline locations move inward when their respective binding energies are increased. This effect is further highlighted when considering how the peak radius of N₂H⁺ and HCO⁺ emission varies with luminosity (Fig. 10).

4.4 Heating sources

The effective temperature T_eff of the protostar has a minor effect on the positions of the CO and H₂O snowlines (Fig. 11). In fact, the most notable effect of T_eff is along the outflow cavity. For the same bolometric luminosity, the models with T_eff = 5000 K heat their outflow cavity to 50 K out to a distance 1.5 times that of the models with T_eff = 1500K. This is because the central source with a higher T_eff produces more photons at short wavelengths, heating up the surrounding material more efficiently. This is consistent with observations that show that the majority of heating from the protostar escapes through the outflow cavity (van Kempen et al. 2009; Yıldız et al. 2015; Murillo et al. 2018b,a).

The protostellar luminosity L_cen has a more significant effect on the temperature of the cloud core than T_eff. Because of this, the protostellar luminosity has a much larger role in determining the positions of the CO and H₂O snowlines. Comparing the luminosity of a protostar with the location of the N₂H⁺ and HCO⁺ emission peaks is commonly used, in both observations and models, to determine whether a protostar has undergone an accretion burst or not (Harsono et al. 2015; Visser et al. 2015; van’t Hoff et al. 2017; Frimann et al. 2017; Hsieh et al. 2018, 2019b; J∅rgensen et al. 2020). The models presented here are indeed consistent with this practice. While luminosity is relevant in setting the location of snowlines, the models presented here suggest that the cloud core conditions need to be constrained in more detail before determining if the snowline location is due solely to the luminosity of the protostar. Figures 10 to 14 show that the presence or lack of a disk, and parameters such as envelope density and binding energy, can significantly change the snowline position even if the luminosity is the same.

Fig. 11 Peak radius of N2H+ (left column), and HCO+ (right column) simulated emission versus central protostellar luminosity for the fiducial models (black hatch) without disk (top row), and with disk (bottom row). The shaded areas show the impact of effective temperature on the peak radius of N2H+ and HCO+. The solid curve shows the relation for i = 45°. For N2H+, the shaded area indicates the range of peak radii for inclinations between 0° (face-on) and 90° (edge-on). For HCO+, the shaded area indicates the range of peak radii for inclinations between 15° and 90°, while the dotted line shows the peak radius versus luminosity for i = 0° (face-on).

4.5 Cloud core structure

The envelope density ρ_env is a key parameter in setting the CO and H₂O snowline locations (Figs. B.3 and B.4). Physically, the density structures regulate the radiation field by redirecting the photons (Fig. B.3, top row), resulting in distinguishable thermal structures in each model. Changes in density are, in a way, also a chemical effect for two reasons. First, the fractional abundances of all molecules in the chemical network are with respect to the H₂ number density. As input, the H₂ number density is set equal to the density profile of the physical model. Thus changing the density profile changes the absolute abundances. Second, the accretion rates of CO and H₂O are a stronger function of density (∝ n²) than the thermal desorption rates. Hence, the density, as well as the temperature, sets the snowline location. The lower the density, the longer the freezeout times. At sufficiently low densities, the freezeout timescale will be so long that the molecule will stay in the gas phase regardless of the temperature. Consequently, the location of the CO and H₂O snowlines will change with the cloud core density (Figs. B.3 and B.4), and snowline positions do not occur at a single, well defined temperature as commonly assumed. The freeze-out timescale is commonly used as a chemical clock to estimate the burst interval (Frimann et al. 2017; Hsieh et al. 2018, 2019b). In such cases, the cloud core density proportional to the desorption rate would need to be constrained in order to use the freeze-out timescale as a reliable clock. Comparing the peak emission radius of N₂H⁺ and HCO⁺ with luminosity when considering different envelope densities (Fig. 12) shows that for the same luminosity, the envelope density can drastically alter the peak emission radius, unless there is a disk present.

As the protostar evolves, the envelope is expected to disperse either due to accretion onto the protostar or by material being pushed out by the outflow. In turn the outflow cavity is expected to widen as the protostar evolves (Arce & Sargent 2006). Because heating from the protostar mainly escapes through the outflow cavity, wider outflow cavities allow more of the envelope to be heated by the protostar. Despite this, the emission peaks of N₂H⁺ and HCO⁺ in our model present no significant shift with increasing outflow cavity opening angle (Fig. 13). Thus, the chemical model results and simulated line emission maps do not suggest that the outflow cavity angle has a strong impact on the location of the CO and H₂O snowlines. However, outflow cavity wall chemistry needs to be added to accurately study how the outflow cavity opening angle affects the water snowline position.

Observational studies have shown that dust grain growth can occur in embedded protostellar sources (e.g., Harsono et al. 2018). The effect of different dust grain size distributions, dust growth, and dynamical processes of dust are not treated in this work. There are, however, previous papers which have studied the effect of dust grain size on snowline locations. For example, Panic & Min (2017) have studied the effect of gas and dust evolution in the disk mid-plane around typical Herbig Ae stars, while Gavino et al. (2021) have studied the effect of dust grain size distribution, grain-size dependent temperature and dust settling on the gas-grain chemistry in the context of a representative T Tauri disk. The Panic & Min (2017) models explore minimum dust grain sizes between 0.01 and 100 µm, and maximum dust grain sizes in the range of 1 mm to 1 km, with a power-law index of 3.5. The modeled region spans r = 0.24−500 AU for an A-type star with 2M_⊙, L_cen = 35 L_⊙ and T_eff = 10 000 K. The Gavino et al. (2021) models present a range of dust grains between 5 nm and 1mm, with a power-law index of 3.5 within a region up to r = 30–250 AU and z/r = 0–0.5, for a pre-main sequence star with L_cen = 0.75 L_⊙ and T_eff = 3900 K. Meanwhile, our models of the cloud core of embedded protostars (r, z = 8790 AU) include a grain-size distribution of 0.5–100 µm for the disk, and 0.5–1 µm for the envelope, both with a power-law index of 3.5, L_cen = 0.01–200 L_⊙ and T_eff = 1500 or 5000 K. Despite the differences, the key results of all model grids can be compared.

The lack of grains ≤0.01 µm in our models avoids wide temperature fluctuations such as those found in Gavino et al. (2021) which would affect the CO and H₂O snowline locations. Panic & Min (2017) note that the temperature profile in the disk depends on both the total gas mass and the size of the smallest grains. In particular, dust growth leads to a decreased temperature given dust settling and reduced τ = 1 surfaces in the gas. Our models confirm that the effect of reducing the envelope density parameter, and thus the total gas mass, leads to lower temperatures in disk, in particular along the mid-plane, even when the dust grain size does not change (Fig. B.4). In contrast, the envelope shows an increase in temperature when the total gas mass is reduced (Figs. B.3 and B.4). The different effects in the disk mid-plane and envelope are most likely product of the dust grain sizes in the disk (0.5–100 µm) and envelope (0.5–1 µm).

While we do not consider the vertical water snowline location, and Gavino et al. (2021) cannot provide constraints on the radial water snowline location (≪30 AU), both models agree on two aspects regarding the water snowline. The H₂O snowline location is regulated by the penetration of UV and the amount of H₂ in the disk; and water in the gas phase is mainly present in the warm regions (upper layers) of the disk. These aspects are shown in our models by the fact that the presence of the disk, and its characteristics, limits the water snowline location, and the gas-phase water is located between the outflow cavity wall and the disk surface. When a disk is not present, and thus the largest grains are 1 µm, the water snowline location shifts outward with increasing luminosity. Gavino et al. (2021) points out that dust temperature impacts chemistry. Our models demonstrate that the presence of a disk, with larger dust grains, significantly changes the temperature profile of the cloud core, and consequently the snowline location of different molecular species. We further highlight that gas density also plays a key role in the distribution of different molecular species.

Given that this work models embedded protostars with significant envelope, while Panic & Min (2017) and Gavino et al. (2021) model disks around Ae Herbig and T Tauri sources, the CO snowline location cannot be directly compared. However, the UV shielding and grain temperature relation to freeze-out effects for CO reported by Gavino et al. (2021) are seen in our models as well. Similarly, the importance of gas mass and smallest grain sizes on the CO snowline location found in Panic & Min (2017) is also reproduced. Changes in envelope gas mass alter the CO snowline, whether a disk is present or not. The presence of a disk, which has a larger range of dust grain sizes, changes the CO snowline location and distribution of CO in the gas phase throughout the cloud core relative to the case without a disk. Indeed, a flared disk further produces a change in the CO gas phase distribution at the edge of the disk along the mid-plane (Fig. 4). In observational studies of embedded protostars this effect is often referred to as disk shadowing (e.g., Murillo et al. 2015), which is basically UV shielding of the envelope by the dust disk causing change in the dust temperature in the envelope along the disk plane. While this comparison provides some insight into the effect of dust grain size distributions on molecular gas, it does not provide the full picture for the case of embedded protostars and snowline locations. For example, it is not clear how different ratios of large to small dust grains within the envelope and disk (i.e., the power-law index), and a wider range of dust grain sizes affect the molecular gas distribution. This topic is left for future work.

Fig. 12 Peak radius of N2H+ (left column), and HCO+ (right column) simulated emission versus central protostellar luminosity for the fiducial models (black hatch) without disk (top row), and with disk (bottom row). The shaded areas show how the measured peak radius shifts with changes in envelope density. The N2H+ curve for ρenv = 105 cm−3 is not indicative of the CO snowline due to the low density, and thus caution must be taken when comparing with observations. See main text for further discussion. The solid curve shows the relation for i = 45°. For N2H+, the shaded area indicates the range of peak radii for inclinations between 0° (face-on) and 90° (edge-on). For HCO+, the shaded area indicates the range of peak radii for inclinations between 15° and 90°, while the dotted line shows the peak radius versus luminosity for i = 0° (face-on).

Fig. 13 Peak radius of N2H+ (left column), and HCO+ (right column) simulated emission versus central protostellar luminosity for the fiducial models (black hatch) without disk (top row), and with disk (bottom row). The shaded areas show the effect of outflow cavity opening angle ⊙out on the peak radius of N2H+ and HCO+. The solid curve shows the relation for i = 45°. For N2H+, the shaded area indicates the range of peak radii for inclinations between 0° (face-on) and 90° (edge-on). For HCO+, the shaded area indicates the range of peak radii for inclinations between 15° and 90°, while the dotted line shows the peak radius versus luminosity for i = 0° (face-on).

4.6 Effect of disk geometry

In all the molecular distribution and simulated line emission maps, the presence of a disk generates a significant change in the temperature structure of the cloud core. When a disk is present, both CO and H₂O snowline positions are located at smaller radii relative to the case without a disk (Fig. 14). The presence of the disk causes N₂H⁺ emission peak to move inward along the disk mid-plane, producing an hour-glass-like morphology. In contrast, the extent of the HCO⁺ emission becomes more compact when a disk is present (Figs. 6 and 5).

The presence of a disk (Fig. 2), its radius, mass, scale height and density (Figs. 4, 6, 14) alter the temperature and chemical structure throughout the core. The molecular distributions alone would suggest that the disk mainly affects regions beyond the disk edge. However, the simulated line emission maps show that even the inner regions are affected by the presence of a disk, no matter how small.

For HCO⁺ and H₂O emission, the presence of a disk regulates the radial extent of their emission and thus the H₂O snowline location (Figs. 5 and 14). These results suggest that the presence of a disk in embedded cloud cores could help limit the radial extent of warm molecules in protostellar cloud cores to regions within the disk. This of course is true if the only heating source is the central protostar. An external heating source, such as accretion shocks or another protostar, may change the distribution of warm molecules through the disk. The first row of Fig. 14 shows that a disk of R_disk = 50 AU and M_diSk = 0.05 M_⊙ causes the HCO⁺ emission to be very compact, unless the protostar reaches luminosities of 100 L_⊙ or higher. This could provide an argument as to why very bright sources, or massive protostars, appear so chemically rich out to large radii in comparison with low mass protostellar sources even when a disk is present. Given the shape of disks, inclination plays a role in whether regions within the H₂O snowline position can be observed. However, we cannot necessarily know the inclination of a disk in advance, especially if it is small, below the spatial resolution of the data, or there is no kinematic data to infer the presence of a disk. If emission from warm molecules is detected then the disk, if present, is close to face-on (i < 45°) and the snowline location can be measured. On the other hand, if warm molecules are not observed (despite other indicators that they should be, e.g., bolometric luminosity), then it might suggest that a disk, if present, is close to edge-on (i > 45°), and measurement of the H₂O snowline location becomes more challenging.

Thus, the disk behaves like an “umbrella.” The disk shields the cloud core from protostellar heating along the mid-plane, while at the same time it retains the warm (and hot) gas at scales within the disk radius. The presence of disks in embedded sources, both rotationally supported and flattened dust structures, with a range of radii and masses has been confirmed observationally (J∅rgensen et al. 2009; Enoch et al. 2011; Murillo et al. 2013; Harsono et al. 2014; Yen et al. 2015, 2017; Persson et al. 2016; Maret et al. 2020; Tobin et al. 2020). In addition, molecular species tracing warm regions have a tendency to be observed in the inner regions of the protostar and along the disk-like structures or the outflow cavity (e.g., Murillo et al. 2018b; Artur de la Villarmois et al. 2019). Further observational evidence of the disk altering the location of cold molecules in embedded sources has also been reported (Murillo et al. 2015).

The presence of a disk dictates where warm and hot gas are located inside the cloud core. Thus, species such as water and complex organic molecules (COMs) released from the dust grains are very likely present in every protostellar source. Their successful detection is dependent on the presence of a disk, in combination with the protostellar luminosity, inclination angle, and spatial resolution. When a disk is present, water and COMs are most likely to be detected at inclinations i < 45°. The necessary spatial resolution will depend on the luminosity of the protostellar system, with the lowest luminosity objects requiring spatial resolutions equivalent to ~10 AU. At i ~ 45°, the degree of disk flaring will affect how well the warm molecular regions can be detected. The impact of a disk on the chemistry of a cloud core could provide some insight into observations of protostellar systems, a few examples are noted here. The concentration of emission from warm molecules to a small region around the protostars (e.g., IRAS 16293-2422, J∅rgensen et al. 2016; Murillo et al. 2022). The presence of a disk, and its inclination, can obscure the warm molecular regions (e.g., NGC 1333 IRAS4A, De Simone et al. 2020). An aspect that would still need to be explored is how warm disks, like those in Taurus (e.g., van’t Hoff et al. 2020), affect the snowline location.

If the disk does indeed help constrain the radial extent of all warm molecular species, this would have an implication for planet formation. Recent studies show that embedded disks have the dust mass needed to produce planets (Tychoniec et al. 2018, 2020), and possibly the dust size as well (Harsono et al. 2018). Hence, the effect seen in the results in this work would suggest that embedded disks not only produce planets, but it is the physical conditions in the embedded phase that set their chemical composition.

Fig. 14 Peak radius of N2H+ (left column), and HCO+ (right column) simulated emission versus central protostellar luminosity for the fiducial models with disk (black hatch). For comparison, the top row shows the fiducial model without disk in red. The shaded areas show the effect of different disk parameters (Rdisk, Mdisk, and H0) on the peak radius of N2H+ and HCO+. The solid curve shows the relation for i = 45°. For N2H+, the shaded area indicates the range of peak radii for inclinations between 0° (face-on) and 90° (edge-on). For HCO+, the shaded area indicates the range of peak radii for inclinations between 15° and 90°, while the dotted line shows the peak radius vs luminosity for i = 0° (face-on).

5 Conclusions

This paper presents a grid of cylindrical symmetric (2D) steady-state physico-chemical models aiming to study the conditions and source parameters that affect the CO and H₂O snowline locations within protostellar cloud cores. The chemical network included deuterated species and the most important ion-molecule reactions for prescribing the abundances of N₂H⁺, HCO⁺, and DCO⁺. A simplified treatment of water, which does not consider the formation of H₂O but only its freeze-out and desorption from dust grains, was used. A range of molecular binding energies for CO and H₂O were used to simulate different ice environments (pure ices versus mixed ices). For physical parameters, the effective temperature of the protostar, luminosity, cloud core density, outflow cavity opening angle, and disk geometry (radius, mass, scale height, and density) were considered. Two fiducial models, one with a disk and one without, were used as references to understand how each individual parameter affected the CO and H₂O snowline positions. We simulated the molecular line emission from snowline location tracers with the purpose of determining which parameters produce observable effects on the snowline locations. With the simulated line emission maps, the robustness of N₂H⁺ and HCO⁺ as snowline position tracers, the impact of inclination, and the spatial resolution on the emission peak positions as measured from the simulated line emission maps were also addressed. Finally, for the purpose of comparing the models with observations, plots of the peak emission radius versus luminosity for the studied parameter space are presented.

The results of this study are as follows:

1. The CO and H₂O snowline locations are mainly dictated by luminosity and cloud core density. Increasing luminosity shifts the snowline location outward to larger radii. In contrast, increasing the protostellar cloud core density causes the snowline location to shift inward to smaller radii, regardless of the protostellar luminosity.

2. The CO snowline location shifts radially outward or inward in all directions when there is no disk present. When a disk is present, the CO snowline position shifts inward along the disk mid-plane. Vertical shifting of the CO snowline position also occurs if the disk is flared.

3. When no disk is present, the H₂O snowline position shifts radially outward as luminosity increases. In contrast, when a disk is present, the radial shift of the H₂O snowline position along the disk mid-plane is limited to radii below the disk radius, concentrating the gas-phase water to small regions around the protostar. The exception to this trend occurs with small disks and high luminosities (R_disk = 50 AU and L_cen > 150 L_⊙ in this work). This effect would also concentrate all warm (and hot) gas-phase molecules in the cloud core to small regions around the protostar.

4. Both inclination of the protostellar cloud core along the line of sight and the spatial resolution of the data affect the observability and measurement of snowline locations, in particular when disks are present. The H₂O snowline location is easier to measure with inclinations i ≲ 45° and spatial resolutions that can resolve ~ 10 AU scales. At higher inclinations i > 45° and lower spatial resolutions, the H₂O snowline location is difficult to measure from HCO⁺ emission. At inclinations i ~ 45°, the degree of disk flaring affects how well the H₂O snowline location can be measured. Measurement of the CO snowline location is not significantly affected by inclination or spatial resolution.

5. Finally, N₂H⁺ and HCO⁺ emission serve as good observational tracers of the CO and H₂O snowline locations. Determining whether a disk is present, additional molecular tracers (e.g., DCO⁺ and line emission from warm or hot molecules, such as methanol), and envelope density (even an upper limit) would provide additional constraints that would help in accurately determining the factors that influence the observed CO and H₂O snowline positions.

The models presented in this work show that the snowline position is not related to a single physical parameter, such as luminosity, but rather the snowline location is dependent on several factors. The physical structure of the protostellar cloud core plays a key role in determining the location of snowlines, and how well these can be observed and measured.

Acknowledgements

N.M.M. acknowledges support from the RIKEN Special Postdoctoral Researcher Program (Fellowships), and a pioneering project in RIKEN (Evolution of Matter in the Universe). C.W. acknowledges financial support from the University of Leeds, the Science and Technology Facilities Council, and UK Research and Innovation (grant numbers ST/T000287/1 and MR/T040726/1).

Appendix A CO isotopologues

Section 3.1 and Figure 3 show the robustness of N₂H⁺ as a tracer of the CO snowline location. In this appendix, N₂H⁺ emission is compared with CO isotopologs to further explore the robustness of N₂H⁺ as a CO snowline tracer. The ¹³CO and C¹⁸O emission are less extended than the ¹²CO emission, which is expected from observations. As the column densities of ¹³CO and C¹⁸O decrease with radius, the optical depth of the lines also decrease, and thus the emission strength declines sharply as well. Hence, r(¹²CO) >> r(¹³CO) >> r(C¹⁸O), where r(X) is the snowline location of X species. It is possible, as well, that C¹⁸O is optically thin throughout the observed region, whereas ¹²CO and ¹³CO will remain optically thick out to some distance.

Fig. A.1 Robustness of N2H+ as a CO snowline location tracer compared to CO isotopolog emission. Simulated emission maps show 13CO (top row), C18O (second row) and N2H+ (third row) emission. Contours in all three rows are N2H+ at 0.30, 0.6, 0.9, 1.2, 1.5 mJy beam−1 km s−1. Different inclinations are shown in order to determine if inclination affects the robustness of N2H+ as a CO snowline position tracers. Inclinations shown are from face-on (i = 0°) to edge-on (i = 90°). The first three columns show the fiducial model without disk, while the other three columns show the fiducial model with disk. Slices extracted along z = 0 from the simulated line emission maps are shown in the bottom row normalized to the peak of each profile.

Appendix B Sample of model grid

The full model grid is available online at https://starformation.space/models/snowline-models/ and are also archived at CDS.. The model grid includes the calculated molecular distributions of the species included in the chemical network used in this work (Table 3). The molecular distributions are provided as ascii tables. Simulated emission maps of N₂H⁺ and HCO⁺ in fits format are also included for seven inclinations: i = 0°, 15°, 25°, 45°, 65°, 75°, and 90°. In addition, the spectral energy distribution for each inclination is provided. Finally, the peak positions of N₂H⁺ and HCO⁺ versus luminosity are provided for each set of models with the same parameters. These are the data used to generate the plots in Figures 8 to 14.

In this appendix a sample of the data available from the model grid is shown. Figures B.1 to B.4 show comparisons of the models without and with disk for changes in binding energy and envelope density. These models are discussed in the main text. Figure B.5 shows the spectral energy distributions (SEDs) for all inclinations used in this work for the two fiducial models. It is interesting to note that the presence of the disk changes the brightness of the SED peak with increasing inclination. Figure B.6 illustrates how the cloud core inclination affects the measured peak position of N₂H⁺ and HCO⁺ with respect to luminosity for the disk, envelope and outflow cavity conditions of the fiducial models. The models presented in Figure B.6 are the same as the black curves in Figures 8 to 14.

Fig. B.1 Effect of binding energy on the CO snowline location. Molecular distributions are shown in the top four rows. The fifth row shows the intensity integrated simulated line emission maps for N2H+ 1–0. The sixth row shows the corresponding slice extracted along z = 0 from the N2H+ simulated emission maps. The average distance of the peaks from the center (i.e., peak radius) is indicated in AU in the top left corner. These positions are shown with orange vertical lines. The second, and fifth columns show the fiducial models without, and with disk, respectively, with Eb,CO = 1150 K and $E_{\text{b},{\text{H}}_2\text{O}}=4820$ K. The column to the left of each fiducial model has $E_{\text{b},{\text{H}}_2\text{O}}=5700$ K, while the right column has Eb,CO = 1307 K. Fractional abundances of all molecular species are relative to total number density of H2 (top row). The black and gray contours show gas temperature and density, respectively. The simulated line emission maps are shown at i = 45° and convolved to a beam of 2″.

Fig. B.2 Effect of binding energy on the H2O snowline location. Molecular distributions are shown in the top three rows. The fourth row shows the intensity integrated simulated line emission maps for HCO+ 3−2. Images that have been scaled for better comparison have the scaling factor on the top right corner. The fifth row shows the corresponding slice extracted along z = 0 from the HCO+ simulated emission maps. The average distance of the peaks from the center (i.e., peak radius) is indicated in AU in the top left corner. These positions are shown with orange vertical lines. No scaling has been applied to the HCO+ profiles. The second, and fifth columns show the fiducial models without, and with disk, respectively, with Eb,CO = 1150 K and $E_{\text{b},{\text{H}}_2\text{O}}=4820$ K. The column to the left of each fiducial model has $E_{\text{b},{\text{H}}_2\text{O}}=5700$ K, while the right column has Eb,CO = 1307 K. Fractional abundances of all molecular species are relative to total number density of H2 (top row). The black and gray contours show gas temperature and density, respectively. The simulated line emission maps are shown at i = 45° and convolved to a beam of 0.25″.

Fig. B.3 Effect of envelope density on the CO snowline location. Molecular distributions are shown in the top four rows. The fifth row shows the intensity integrated simulated line emission maps for N2H+ 1−0. Images that have been scaled for better comparison have the scaling factor on the top right corner. The sixth row shows the corresponding slice extracted along z = 0 from the N2H+ simulated emission maps. The average distance of the peaks from the center (i.e., peak radius) is indicated in AU in the top left corner. These positions are shown with orange vertical lines. No scaling has been applied to the N2H+ profiles. The second and fifth columns show the fiducial model without and with disk, respectively, having a density of ρenv = 106 cm−3. First and fourth columns show densities of ρenv = 105 cm−3, while the third and sixth columns show ρenv = 107 cm−3. Fractional abundances of all molecular species are relative to total number density of H2 (top row). The black and gray contours show gas temperature and density, respectively. The simulated line emission maps are shown at i = 45° and convolved to a beam of 2″.

Fig. B.4 Effect of envelope density on the H2O snowline location. Molecular distributions are shown in the top three rows. The fourth row shows the intensity integrated simulated line emission maps for HCO+ 3−2. Images that have been scaled for better comparison have the scaling factor on the top right corner. The fifth row shows the corresponding slice extracted along z = 0 from the HCO+ simulated emission maps. The average distance of the peaks from the center (i.e., peak radius) is indicated in AU in the top left corner. These positions are shown with orange vertical lines. No scaling has been applied to the HCO+ profiles. The second and fifth columns show the fiducial model without and with disk, respectively, having a density of ρenv = 106 cm−3. First and fourth columns show densities of ρenv = 105 cm−3, while the third and sixth columns show ρenv = 107 cm−3. Fractional abundances of all molecular species are relative to total number density of H2 (top row). The black and gray contours show gas temperature and density, respectively. The simulated line emission maps are shown at i = 45° and convolved to a beam of 0.25″.

Fig. B.5 Calculated spectral energy distributions (SEDs) for the fiducial models, without disk (left) and with disk (right). The different colored lines indicate inclination.

Fig. B.6 Peak position of N2H+ (left) and HCO+ (right) versus luminosity for the disk, envelope and outflow cavity conditions of the fiducial models. The different colors indicate inclination. These curves are combined into the black curve shown in Figures 8 to 14. We note that the peak position of N2H+ does not show a clear trend with inclination, while the peak position of HCO+ decreases with increasing highlighting the importance of inclination in measuring an accurate water snowline location.

References

All Tables

All Figures

	Fig. 1 Cartoon showing the structure of the physical models, along with the illustrated definitions of the disk radius Rdisk, disk scale height H0 at Rdisk, and outflow cavity opening angle Θoutflow.
In the text

	Fig. 2 Density (top row) and temperature (bottom row) distribution of the two fiducial models: without a disk (left) and with a disk (center and right). The right column shows a zoom in of the fiducial model with the disk marked with an orange box in the center column.
In the text

Fig. 3 Robustness of N2H+ and DCO+ as CO snowline location tracers. Integrated intensity maps of the simulated line emission show CO (top row), DCO+ (second row) and N2H+ (third row) emission. Contours in all three rows are for CO at 10, 30, 50, and 100 mJy beam−1 km s−1. Images that have been scaled for better comparison have the scaling factor noted on the top right corner. Different inclinations are shown in order to determine if inclination affects the robustness of N2H+ and DCO+ as snowline position tracers. Inclinations shown are from face-on (i = 0°) to edge-on (i = 90°). The first three columns show the fiducial model without disk, while the other three columns show the fiducial model with disk. Slices extracted along z = 0 from the simulated line emission maps are shown in the bottom row normalized to the peak of each profile.

In the text

Fig. 4 Effect of disk geometry parameters on the CO snowline location. Molecular distributions are shown in the top four rows. The fifth row shows the maps of the integrated intensity of the simulated line emission for N2H+ 1−0. The sixth row shows the corresponding slice extracted along z = 0 from the N2H+ simulated emission maps. The average distance of the peaks from the center (i.e., peak radius) is indicated in AU in the top left corner. These positions are shown with orange vertical lines. The first and third columns show the fiducial models without and with disk, respectively. The additional columns show the effect of changing disk radius (second column), disk mass (fourth column), scale height H0 (fifth column), and disk density (sixth column). Fractional abundances of all molecular species are relative to total number density of H2 (top row). The black and gray contours show gas temperature and density, respectively. The emission maps are shown at i = 45° and convolved to a beam of 2″.

In the text

Fig. 5 Robustness of HCO+ as a H2O snowline location tracer. Simulated emission maps show HCO+ (top row), p–H2O (second row) and o–H2O (third row) emission. Contours in all three rows are H2O at 0.4, 0.6, 0.9, 1.6, 1.9, and 2.1 mJy beam−1 km s−1. Images that have been scaled for better comparison have the scaling factor on the top right corner. Different inclination are shown in order to determine if it affects the robustness of HCO+ as a snowline location tracer. Inclinations shown are from face-on (i = 0°) to edge-on (i = 90°). The first three columns show the fiducial model without disk, while the other three columns show the fiducial model with disk. Slices extracted along z = 0 from the simulated line emission maps are shown in the bottom row normalized to the peak of each profile.

In the text

Fig. 6 Effect of disk geometry parameters on the H2O snowline location. Molecular distributions are shown in the top three rows. The fourth row shows the simulated intensity integrated line emission maps for HCO+ 3−2. Images that have been scaled for better comparison have the scaling factor on the top right corner. The fifth row shows the corresponding slice extracted along z = 0 from the HCO+ simulated emission maps. The average distance of the peaks from the center (i.e., peak radius) is indicated in AU in the top left corner. These positions are shown with orange vertical lines. No scaling has been applied to the HCO+ profiles. The first and third columns show the fiducial models without and with disk, respectively. The additional columns show the effect of changing disk radius Rdisk (second column), disk mass (fourth column), scale height H0 (fifth column), and disk density (sixth column). Fractional abundances of all molecular species are relative to total number density of H2 (top row). The black and gray contours show gas temperature and density, respectively. The simulated line emission maps are shown at i = 45° and convolved to a beam of 0.25″.

In the text

Fig. 7 Simulated line emission maps showing the effect of spatial resolution (0.1″, 0.25″, and 2″) on the measured emission peak positions of N2H+ and HCO+. All models have an inclination of 45°. The first three columns show the fiducial model without disk. The other three columns show the fiducial model with disk. The second and third rows show corresponding slices extracted along the z = 0, with the number on the top left corner indicating the distance in AU of the peak from the center (i.e., peak radius). These positions are shown with orange vertical lines.

In the text

Fig. 8 Peak radius of N2H+ (left column), and HCO+ (right column) simulated emission versus luminosity for the fiducial models without disk (top row), and with disk (bottom row). The black hatched area shows the central protostellar luminosity, while the colored area shows the bolometric luminosity. The solid curve shows the relation for i = 45°. For N2H+, the shaded area indicates the range of peak radii for inclinations between 0° (face-on) and 90° (edge-on). For HCO+, the shaded area indicates the range of peak radii for inclinations between 15° and 90°, while the dotted line shows the peak radius versus luminosity for i = 0° (face-on).

In the text

Fig. 9 Peak radius of N2H+ (left column), and HCO+ (right column) simulated emission versus central protostellar luminosity for the fiducial models (black hatch) without disk (top row), and with disk (bottom row). The different shaded areas show how the measured peak radius is altered by the spatial resolution used for observations. The solid curve shows the relation for i = 45°. For N2H+, the shaded area indicates the range of peak radii for inclinations between 0° (face-on) and 90° (edge-on). For HCO+, the shaded area indicates the range of peak radii for inclinations between 15° and 90°, while the dotted line shows the peak radius versus luminosity for i = 0° (edge-on).

In the text

Fig. 10 Peak radius of N2H+ (left column), and HCO+ (right column) simulated emission versus central protostellar luminosity for the fiducial models (black hatch) without disk (top row), and with disk (bottom row). The shaded areas show the effect of CO and H2O binding energy on the peak radius of N2H+ and HCO+. The solid curve shows the relation for i = 45°. For N2H+, the shaded area indicates the range of peak radii for inclinations between 0° (face-on) and 90° (edge-on). For HCO+, the shaded area indicates the range of peak radii for inclinations between 15° and 90°, while the dotted line shows the peak radius versus luminosity for i = 0° (face-on).

In the text

Fig. 11 Peak radius of N2H+ (left column), and HCO+ (right column) simulated emission versus central protostellar luminosity for the fiducial models (black hatch) without disk (top row), and with disk (bottom row). The shaded areas show the impact of effective temperature on the peak radius of N2H+ and HCO+. The solid curve shows the relation for i = 45°. For N2H+, the shaded area indicates the range of peak radii for inclinations between 0° (face-on) and 90° (edge-on). For HCO+, the shaded area indicates the range of peak radii for inclinations between 15° and 90°, while the dotted line shows the peak radius versus luminosity for i = 0° (face-on).

In the text

Fig. 12 Peak radius of N2H+ (left column), and HCO+ (right column) simulated emission versus central protostellar luminosity for the fiducial models (black hatch) without disk (top row), and with disk (bottom row). The shaded areas show how the measured peak radius shifts with changes in envelope density. The N2H+ curve for ρenv = 105 cm−3 is not indicative of the CO snowline due to the low density, and thus caution must be taken when comparing with observations. See main text for further discussion. The solid curve shows the relation for i = 45°. For N2H+, the shaded area indicates the range of peak radii for inclinations between 0° (face-on) and 90° (edge-on). For HCO+, the shaded area indicates the range of peak radii for inclinations between 15° and 90°, while the dotted line shows the peak radius versus luminosity for i = 0° (face-on).

In the text

Fig. 13 Peak radius of N2H+ (left column), and HCO+ (right column) simulated emission versus central protostellar luminosity for the fiducial models (black hatch) without disk (top row), and with disk (bottom row). The shaded areas show the effect of outflow cavity opening angle ⊙out on the peak radius of N2H+ and HCO+. The solid curve shows the relation for i = 45°. For N2H+, the shaded area indicates the range of peak radii for inclinations between 0° (face-on) and 90° (edge-on). For HCO+, the shaded area indicates the range of peak radii for inclinations between 15° and 90°, while the dotted line shows the peak radius versus luminosity for i = 0° (face-on).

In the text

Fig. 14 Peak radius of N2H+ (left column), and HCO+ (right column) simulated emission versus central protostellar luminosity for the fiducial models with disk (black hatch). For comparison, the top row shows the fiducial model without disk in red. The shaded areas show the effect of different disk parameters (Rdisk, Mdisk, and H0) on the peak radius of N2H+ and HCO+. The solid curve shows the relation for i = 45°. For N2H+, the shaded area indicates the range of peak radii for inclinations between 0° (face-on) and 90° (edge-on). For HCO+, the shaded area indicates the range of peak radii for inclinations between 15° and 90°, while the dotted line shows the peak radius vs luminosity for i = 0° (face-on).

In the text

Fig. A.1 Robustness of N2H+ as a CO snowline location tracer compared to CO isotopolog emission. Simulated emission maps show 13CO (top row), C18O (second row) and N2H+ (third row) emission. Contours in all three rows are N2H+ at 0.30, 0.6, 0.9, 1.2, 1.5 mJy beam−1 km s−1. Different inclinations are shown in order to determine if inclination affects the robustness of N2H+ as a CO snowline position tracers. Inclinations shown are from face-on (i = 0°) to edge-on (i = 90°). The first three columns show the fiducial model without disk, while the other three columns show the fiducial model with disk. Slices extracted along z = 0 from the simulated line emission maps are shown in the bottom row normalized to the peak of each profile.

In the text

Fig. B.1 Effect of binding energy on the CO snowline location. Molecular distributions are shown in the top four rows. The fifth row shows the intensity integrated simulated line emission maps for N2H+ 1–0. The sixth row shows the corresponding slice extracted along z = 0 from the N2H+ simulated emission maps. The average distance of the peaks from the center (i.e., peak radius) is indicated in AU in the top left corner. These positions are shown with orange vertical lines. The second, and fifth columns show the fiducial models without, and with disk, respectively, with Eb,CO = 1150 K and $E_{\text{b},{\text{H}}_2\text{O}}=4820$ K. The column to the left of each fiducial model has $E_{\text{b},{\text{H}}_2\text{O}}=5700$ K, while the right column has Eb,CO = 1307 K. Fractional abundances of all molecular species are relative to total number density of H2 (top row). The black and gray contours show gas temperature and density, respectively. The simulated line emission maps are shown at i = 45° and convolved to a beam of 2″.

In the text

Fig. B.2 Effect of binding energy on the H2O snowline location. Molecular distributions are shown in the top three rows. The fourth row shows the intensity integrated simulated line emission maps for HCO+ 3−2. Images that have been scaled for better comparison have the scaling factor on the top right corner. The fifth row shows the corresponding slice extracted along z = 0 from the HCO+ simulated emission maps. The average distance of the peaks from the center (i.e., peak radius) is indicated in AU in the top left corner. These positions are shown with orange vertical lines. No scaling has been applied to the HCO+ profiles. The second, and fifth columns show the fiducial models without, and with disk, respectively, with Eb,CO = 1150 K and $E_{\text{b},{\text{H}}_2\text{O}}=4820$ K. The column to the left of each fiducial model has $E_{\text{b},{\text{H}}_2\text{O}}=5700$ K, while the right column has Eb,CO = 1307 K. Fractional abundances of all molecular species are relative to total number density of H2 (top row). The black and gray contours show gas temperature and density, respectively. The simulated line emission maps are shown at i = 45° and convolved to a beam of 0.25″.

In the text

Fig. B.3 Effect of envelope density on the CO snowline location. Molecular distributions are shown in the top four rows. The fifth row shows the intensity integrated simulated line emission maps for N2H+ 1−0. Images that have been scaled for better comparison have the scaling factor on the top right corner. The sixth row shows the corresponding slice extracted along z = 0 from the N2H+ simulated emission maps. The average distance of the peaks from the center (i.e., peak radius) is indicated in AU in the top left corner. These positions are shown with orange vertical lines. No scaling has been applied to the N2H+ profiles. The second and fifth columns show the fiducial model without and with disk, respectively, having a density of ρenv = 106 cm−3. First and fourth columns show densities of ρenv = 105 cm−3, while the third and sixth columns show ρenv = 107 cm−3. Fractional abundances of all molecular species are relative to total number density of H2 (top row). The black and gray contours show gas temperature and density, respectively. The simulated line emission maps are shown at i = 45° and convolved to a beam of 2″.

In the text

Fig. B.4 Effect of envelope density on the H2O snowline location. Molecular distributions are shown in the top three rows. The fourth row shows the intensity integrated simulated line emission maps for HCO+ 3−2. Images that have been scaled for better comparison have the scaling factor on the top right corner. The fifth row shows the corresponding slice extracted along z = 0 from the HCO+ simulated emission maps. The average distance of the peaks from the center (i.e., peak radius) is indicated in AU in the top left corner. These positions are shown with orange vertical lines. No scaling has been applied to the HCO+ profiles. The second and fifth columns show the fiducial model without and with disk, respectively, having a density of ρenv = 106 cm−3. First and fourth columns show densities of ρenv = 105 cm−3, while the third and sixth columns show ρenv = 107 cm−3. Fractional abundances of all molecular species are relative to total number density of H2 (top row). The black and gray contours show gas temperature and density, respectively. The simulated line emission maps are shown at i = 45° and convolved to a beam of 0.25″.

In the text

	Fig. B.5 Calculated spectral energy distributions (SEDs) for the fiducial models, without disk (left) and with disk (right). The different colored lines indicate inclination.
In the text

Fig. B.6 Peak position of N2H+ (left) and HCO+ (right) versus luminosity for the disk, envelope and outflow cavity conditions of the fiducial models. The different colors indicate inclination. These curves are combined into the black curve shown in Figures 8 to 14. We note that the peak position of N2H+ does not show a clear trend with inclination, while the peak position of HCO+ decreases with increasing highlighting the importance of inclination in measuring an accurate water snowline location.

In the text