Modeling snowline locations in protostars: The impact of the structure of protostellar cloud cores

Abridged Context: Snowlines during star and disk formation are responsible for a range of effects during the evolution of protostars, such as setting the chemical composition of the envelope and disk. This in turn influences the formation of planets by changing the elemental compositions of solids and affecting the collisional properties and outcomes of dust grains. Snowlines can also reveal accretion bursts, providing insight into the formation process of stars. Methods: A numerical chemical network coupled with a grid of cylindrical-symmetric physical models was used to identify what physical parameters alter the CO and H$_2$O snowline locations. The investigated parameters are the initial molecular abundances, binding energies of CO and H$_2$O, heating source, cloud core density, outflow cavity opening angle, and disk geometry. Simulated molecular line emission maps were used to quantify the change in the snowline location with each parameter. Conclusions: The models presented in this work show that the CO and H$_2$O snowline locations do not occur at a single, well-defined temperature as is commonly assumed. Instead, the snowline position depends on luminosity, cloud core density, and whether a disk is present or not. Inclination and spatial resolution affect the observability and successful measurement of snowline locations. We note that N$_2$H$^+$ and HCO$^+$ emission serve as good observational tracers of CO and H$_2$O snowline locations. However, constraints on whether or not a disk is present, the observation of additional molecular tracers, and estimating envelope density will help in accurately determining the cause of the observed snowline position. Plots of the N$_2$H$^+$ and HCO$^+$ peak emission radius versus luminosity are provided to compare the models with observations of deeply embedded protostars aiming to measure the CO and H$_2$O snowline locations.


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 dy-namical 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. 2019Tychoniec et al. , 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. 2015Yen et al. , 2017Maret 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. 2015Murillo et al. , 2018bvan '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. 2016Murillo et al. , 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. 2015a;S255IR-SMA1 Safron et al. 2015b;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 physicochemical models that explore and compare the impact of indi-vidual 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 Bjerkeli et al. 2016;Rab et al. 2017), have limited chemical modeling 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 . The effects of gas dispersal and dust evolution on the CO snowline location has been studied in disks around Herbig stars (Panić & Min 2017). The effect of dust grain sizes on molecular gas distributions and the H 2 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 2 O snowlines, as well as quantifies the resulting abundance distribution and emergent line emission from their respective tracers, N 2 H + and HCO + . Since CO destroys N 2 H + via an ion-molecule reaction (N 2 H + + CO − −− → HCO + + N 2 ), N 2 H + increases in abundance when CO freezes out onto the dust grains, thus tracing the CO snowline location. A similar process occurs between H 2 O and HCO + (HCO + + H 2 O − −− → H 3 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 Sections 4 and 5, respectively.

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 Figure 1. The density profile is built starting from a rotating flattened envelope introduced by Ulrich (1976) ρ env (r, µ) = ρ env,0 r r cen where r and θ are the spherical coordinates, ρ env,0 is the density parameter used to scale the density profiles, r cen is the centrifugal Murillo, Hsieh & Walsh: Modeling snowline locations in protostars  Anderl et al. 2016;E b,CO , Collings et al. 2003, Collings et al. 2004, Bisschop et al. 2006, and Noble et al. 2012andE b,H 2 O , Sandford &Allamandola 1993 andFraser et al. 2001.
radius, and µ = cos(θ). The parameter µ 0 = cos(θ 0 ) is a solution of µ 3 0 + µ 0 (r/r cen − 1) − µ(r/r cen ) = 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 with Σ(R) defined as where R and z are the cylindrical coordinates, R disk is the disk radius in au, and Σ 0 is the disk surface density at R disk . The disk scale height H(R) is defined as a function of radius R (Chiang & Goldreich 1997), for which H 0 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 2 density (cm −3 ) to dust density (g cm −3 ), 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 1 (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, A&A proofs: manuscript no. 42982corr 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 where Θ outflow is the outflow cavity opening angle. We adopt z 0 ≡ 50 au.
Two-dimensional density and temperature profiles are generated with the above described structures. The dust temperature is calculated with RADMC3D 2 (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 .

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 2 O snowlines within the cloud core of embedded protostellar sources. The parameters which characterize the protostar, outflow, disk, and en-2 http://www.ita.uni-heidelberg.de/ dullemond/software/radmc-3d/ velope 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 solartype 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 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. and the outflow opening angle can be different depending on the tracers used to make the measurement. The outflow opening angle is defined in Equation 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;Testi et al. 2016;Yen et al. 2017;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 0 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).

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 6 cm −3 , 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. 2015Yen et al. , 2017Maury 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 = 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.

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 2 H + , and associated deuterated ions such as H 2 D + , DCO + , and N 2 D + (see Table 3). As we are interested in HCO + as an accessible tracer of the H 2 O snowline location, and N 2 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 2 O, CO and N 2 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. 2013a,b, 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 3 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 2 O + HCO + . The dissociative recombination of H 3 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 2 , and H 2 O (Table 1). Initial abundances are set for CO, N 2 , H 2 O, HD, and neutral grains, and naturally, for H 2 and H. The total initial abundances of CO, N 2 , H 2 O, HD and neutral grains relative to H 2 are set to 2×10 −4 , 1×10 −4 , 2×10 −4 , 1.6×10 −5 , and 1.3×10 −12 , respectively. The initial abundances for CO, N 2 and H 2 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 2 is kept constant for all models, and is set at E b,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(Collings et al. , 2004Bisschop et al. 2006;Noble et al. 2012); E b,H 2 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 8 years to ensure that steady state has been reached to produce molecular distri-butions 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 8 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.

Simulated line emission maps
In order to examine whether the effects of physical and chemical structure on the CO and H 2 O snowline locations are observable, simulated emission maps are generated from the molecular distributions. The molecules N 2 H + and HCO + are typically used in observations to determine the CO and H 2 O snowline positions, respectively (van 't Hoff et al. 2017;Frimann et al. 2017;Hsieh et al. 2018Hsieh et al. , 2019b. Simulated emission maps of CO (three isotopologues: 12 CO, 13 CO, and C 18 O), N 2 H + , DCO + , H 2 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 2 H + and HCO + for all the models discussed in Section 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 • (faceon), 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 (17580 × 17580 AU) is used for N 2 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 where G is the gravitational constant, M is the central source mass, r is the radius, θ is the angle in polar coordinates, and θ 0 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) For the radial and meridional velocities, a dimensionless factor is applied to the above equations as a function of θ and r, ( r r cen ) 2θ/π . 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 2 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 E up were used, one at 183 GHz with E up = 205 K (ALMA Band 5), and a second at 557 GHz with E up = 61 K (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.

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 Section 2.3. The effect of inclination and spatial resolution on the observed CO and H 2 O snowline locations will then be considered. Finally, plots of peak radius of N 2 H + and HCO + versus luminosity for the explored parameter space are provided and described.
A sample of the models from the full grid is shown in Figures 4 and 6, as well as in Appendix A. The full grid of molecu-lar distributions and simulated line emission maps are available online 3 .

CO snowline location
The two fiducial models are used to explore the robustness of N 2 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 2 H + and DCO + . From the z = 0 slices, N 2 H + begins to increase at a radius where CO is between 60% and 70% of its maximum brightness. The N 2 H + emission then peaks where CO is at ∼1% of its maximum brightness. At low envelope densities (e.g., 10 5 cm −3 ), N 2 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 2 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 2 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 2 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 2 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 Section 4.2. Figure A.1 shows how the emission from the CO isotopologs, 13 CO and C 18 O, does not coincide well with the snowline location traced by the N 2 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 ∼21 K 3 The full model grid can be found at https://starformation.space/models/snowline-models/ and are also archived at CDS.
A&A proofs: manuscript no. 42982corr High spatial resolution is needed to detect the shift, and to differentiate it from other effects. (b) The presence of a disk in the cloud core, and its geometry, constrain the extent of the region where water is in the gas phase, particularly in the radial direction. In the molecular distribution models the presence of the disk does not significantly change the water snowline location. However, the geometry of the disk affects the observability and measurability of the H 2 O snowline location. 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 6 cm −3 ) the CO snowline location moves outward to larger radii, whereas at higher densities (ρ env > 10 6 cm −3 ) 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 2 (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 2 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 2 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 0 . 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 2 H + due to disk mass M disk and radius R disk (Fig. 4 bottom rows). Interestingly, changing H 0 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 = 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.

H 2 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 2 O and p−H 2 O emission are examined. HCO + emission begins to increase at the H 2 O half maximum, and peaks 20 to 40 AU further out. Similar to N 2 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 Section 3.1, the limited chemical network in this work does not account for all the destruction and formation pathways for HCO + and H 2 O, in addition to the simple treatment of water in the network. The H 2 O snowline location is altered by the binding energy E b,H 2 O , luminosity, and envelope density. Increasing the binding energy of water (E b,H 2 O = 4820 → 5700 K) causes the H 2 O snowline position to shift inward between 20 to 50 AU (Fig. B.2). The simulated emission line maps show that the effect of binding energy on the H 2 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 2 O snowline as well. This is expected, given how much more luminous pro-tostellar sources have relatively larger warm ( 100 K) regions than less luminous sources. Changing the envelope density has the same effect on the H 2 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 2 O snowline location shifts proportional to luminosity, and inversely proportional to envelope density ρ env . Overall, the H 2 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 0 has no apparent effect on the actual location of the H 2 O snowline (Fig. 6). In a similar manner, the disk density does not seem to have a significant effect on the H 2 O snowline location in the molecular distributions. However, the disk geometry does help to concentrate the H 2 O gas extent between the outflow cavity wall and disk surface ( Fig. 6 fifth and sixth columns). Hence, a very flared disk (H 0 = 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 2 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

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 2 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 2 H + emission decreases with increasing inclination (Fig. B.6). This trend is more evident when a disk is present (Figs. 8 to 14). For i = 0 • to i = 90 • , the peaks of N 2 H + emission shift less than 500 AU in the case without a disk. In the case with a disk, the N 2 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 2 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 • (Fig. 13). 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 2 O is affected by inclination, more so when a disk is present. Consequently, the H 2 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 sim-  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 R disk (second column), disk mass M disk (fourth column), scale height H 0 (fifth column), and disk density (sixth column). Fractional abundances of all molecular species are relative to total number density of H 2 (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 . ∼600, 70, and 30 AU, respectively. The peak of N 2 H + emission in the case without a disk does not shift with different spatial resolutions. In the case with a disk, the N 2 H + peak moves inward (<100 AU) with increasing beam size. Lower spatial resolutions are better suited to trace the CO snowline location with N 2 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 (Section 3), the H 2 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 2 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 ≥ 10 L , but would benefit from high signalto-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 lumi-nous sources, or those that have undergone a strong luminosity burst in the recent past.

N 2 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 2 H + and HCO + line emission peak radius are provided in Figures 8 to 14. Since inclination impacts the measured peak of N 2 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 2 H + and HCO + . Figure 9 shows the three spatial resolutions (beam sizes) versus luminosity. The three spatial resolutions are described in Section 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.
The effect of low envelope density on N 2 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 Section 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 2 H + curve at ρ env = 10 5 cm −3 in Figure 12 is not indicative of the CO snowline location. The curve drops out at luminosities >50 L because the N 2 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 2 H + and CO (ID 13 in Table 3) is no longer a relevant destruction path for N 2 H + . Hence N 2 H + is present in the gas phase in a larger region of the envelope and not correlated to the distribution of CO. The N 2 H + curve at ρ env = 10 5 cm −3 in Figure 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 2 O snowline position. This effect is somewhat dependent on luminosity and density. For luminosities < 10 L , the H 2 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 2 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 midplane 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.

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 2 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 2 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.

Robustness of N 2 H + and HCO + as snowline position tracers
Observational studies of CO and H 2 O snowline positions typically use the peak emission of N 2 H + and HCO + to trace the respective snowline locations (e.g., Qi et al. 2015;Hsieh et al. 2018;van 't Hoff et al. 2018a;Hsieh et al. 2019b;Qi et al. 2019). The simulated line emission maps discussed in Section 3 examined the robustness of N 2 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 2 H + and DCO + are robust tracers of the CO snowline location, and HCO + of the H 2 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 2 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 2 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 2 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 ar- gued to be a good tracer of the H 2 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 2 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 2 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 2 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 2 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-4 transition (Murillo et al. 2018b). However, combining N 2 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 2 H + (or N 2 D + ) is significantly less abundant and difficult to detect because of the depletion of gas-phase N 2 onto dust grains (e.g., VLA 1623-2417: Murillo et al. 2015). Similar to N 2 H + , HCO + emission peaks some distance away from the H 2 O snowline location. However, HCO + can still be a good observational proxy and provide an outer limit to the H 2 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 2 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 2 O snowline. However, similar to HCO + , their success- ful detection could be dependent on source inclination, spatial resolution, and presence of a disk.

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., Martín-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 (Fig. B.1 and B.2), showing that the CO and H 2 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 2 H + and HCO + emission varies with luminosity (Fig. 10).

Heating sources
The effective temperature T eff of the protostar has a minor effect on the positions of the CO and H 2 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 = 1500 K. 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 2 O snowlines. Comparing the luminosity of a protostar with the location of the N 2 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 Visser et al. 2015;van 't Hoff et al. 2017;Frimann et al. 2017;Hsieh et al. 2018Hsieh et al. , 2019bJø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. , 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 N 2 H + , 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 • (edge-on).

Cloud core structure
The envelope density ρ env is a key parameter in setting the CO and H 2 O snowline locations (Fig. 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 2 number density. As input, the H 2 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 2 O are a stronger function of density (∝ n 2 ) 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 2 O snowlines will change with the cloud core density (Fig. 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. 2018Hsieh et al. , 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 2 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 2 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 2 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. disk. The Panić & Min (2017) models explore minimum dust grain sizes between 0.01 to 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 2 M , L cen = 35 L and T e f f = 10 000 K. The Gavino et al. (2021) models present a range of dust grains between 5 nm to 1 mm, 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 e f f = 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 e f f = 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 2 O snowline locations. Panić & 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 & 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 2 O snowline location is regulated by the penetration of UV and the amount of H 2 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 Panić & 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 Panić & 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.

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 2 O snowline positions are located at smaller radii relative to the case without a disk (Fig. 14). The presence of the disk causes N 2 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 2 O emission, the presence of a disk regulates the radial extent of their emission and thus the H 2 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 2 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 2 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;Persson et al. 2016;Yen et al. 2017;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., NGC1333 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(Tychoniec et al. , 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 phys- ical conditions in the embedded phase that set their chemical composition.

Conclusions
This paper presents a grid of cylindrical symmetric (2D) steadystate physico-chemical models aiming to study the conditions and source parameters that affect the CO and H 2 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 2 H + , HCO + , and DCO + . A simplified treatment of water, which does not consider the formation of H 2 O but only its freeze-out and desorption from dust grains, was used. A range of molecular binding energies for CO and H 2 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 2 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 2 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 emis-sion 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 2 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 2 O snowline position shifts radially outward as luminosity increases. In contrast, when a disk is present, the radial shift of the H 2 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.  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 sev-eral 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.    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 BE CO = 1150 K and BE H 2 O = 4820 K. The column to the left of each fiducial model has BE H 2 O = 5700 K, while the right column has BE CO = 1307 K. Fractional abundances of all molecular species are relative to total number density of H 2 (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 .  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 = 10 6 cm −3 . First and fourth columns show densities of ρ env = 10 5 cm −3 , while the third and sixth columns show ρ env = 10 7 cm −3 . Fractional abundances of all molecular species are relative to total number density of H 2 (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 .
Article number, page 28 of 29 Murillo, Hsieh & Walsh: Modeling snowline locations in protostars   Figures 8 to 14. We note that the peak position of N 2 H + does not show a clear trend with inclination, while the peak position of HCO + decreases with increasing i, highlighting the importance of inclination in measuring an accurate water snowline location.