Physical and Chemical Structure of the Disk and Envelope of the Class 0/I Protostar L1527

To make the modeling effort tractable, the software package we have assembled, RadChemT, breaks the calculation into three parts with several key simplifying assumptions. Step 1 carries out MCRT using Hochunk3d, as described in Section 2.2, to calculate dust temperature based on the current-time luminosity L_int and 2D input physical model. Namely, for a given snapshot in time, both density ρ(r, θ, ϕ) and velocity v (r, θ, ϕ) are specified for the main components during the dynamical collapse of the class 0/I system: outflow, rotationally flattened envelope, and disk. The physical model is tuned to L1527, based on fitting modeling and imaging data from Tobin et al. (2008, 2010, 2012).

Step 2 performs a full time-dependent chemical abundance calculation in a point model, i.e., locally for each cell as described in Section 2.3, using the density and temperature from Step 1. As a simplification, this first version of RadChemT does not track material along a gas streamline as a function of time. Instead, the abundance calculations at each spatial grid point are time-dependent but assume fixed density and temperature. The assumption is that current local conditions (ρ, T) sufficiently describe earlier times. The assumption is most valid in the outer envelope, where density and temperature change the least for a moving gas parcel. There is an initial (brief) transient behavior in the chemistry, from starting with elements in atomic form, until molecular cloud abundances are achieved. After this, the chemical age should be the same as the protostar age. The age of an individual protostar can be estimated using ${t}_{\mathrm{age}}={M}_{* }/{\dot{M}}_{\mathrm{env}}$, but the age of an individual object is not well known during this early stage. We therefore compare chemical abundance and morphology at two time stamps, 10⁴ and 10⁵ yr, that roughly span the range of relevant ages for the Class 0/I phase. Predicted abundances that are relatively constant over that time span give confidence that conclusions based on those abundances will be robust.

Step 3 visualizes the protostar in a PV cube, making use of the Doppler effect in spectral line emission to sense the 3D velocity structure. This investigation uses RADMC-3D, as described in Section 2.4, to perform spectral line radiative transfer assuming LTE. The main inputs are the density, temperature, and velocity grids from Step 1, combined with an abundance grid from Step 2 for a selected molecule at t = 10⁵ yr. From the inputs, RADMC-3D creates a synthetic model PV cube in FITS format having units of Jy pixel⁻¹. The primary and synthesized beams are then applied to the model spectral line cube, in order to compare the model with ALMA and CARMA submillimeter observations.

The physical structure of our representative Class 0/I model is based on the Hochunk3d code (Whitney et al. 2013). This model uses an MCRT calculation to obtain the temperature of the dust and the gas density profile. The 2D geometry considers the system in spherical coordinates with axial symmetry and a logarithmically increasing grid in radius. The radial domain goes from the dust destruction radius, which is 10.04R_*, to the outer edge of the envelope, R_env = 12,500 au. The meridional domain covers the full 180°. The model includes a description of the density structure of the disk component and of the envelope with an outflow cavity. Figure 1 shows a schematic of the TSC geometry. Both the TSC and UCM models have similar parameters to tune; Table 1 tabulates adopted values for L1527 to be described more specifically in Sections 2.2.4 and 2.2.5.

Table 1. Physical Parameters

Parameters	Description	TSC	UCM
R* (R⊙)	Stellar radius	1.70	:
T* (K)	Stellar temperature	3300	:
M* (M⊙)	Stellar mass	0.22	:
Mdisk (M⊙)	Mass of the disk	0.011	0.006
Rdisk (au)	Disk outer radius	75	:
${\dot{M}}_{{\rm{d}}{\rm{i}}{\rm{s}}{\rm{k}}}\,({M}_{\odot }\,{{\rm{y}}{\rm{r}}}^{-1})$	Disk accretion rate	6.6 × 10−7	:
${\dot{M}}_{{\rm{e}}{\rm{n}}{\rm{v}}}\,({M}_{\odot }\,{{\rm{y}}{\rm{r}}}^{-1})$	Envelope infall rate	3.0 × 10−6	5.0 × 10−6
θ1 (deg)	Opening angle of the inner cavity surface	15	:
z (au)	z-intercept, inner cavity surface at ω = 0	75	:
θ2 (deg)	Opening angle of the outer cavity surface	6	:
LISRF (L⊙)	Luminosity due to ISRF	0.49	:
Quantities shown below are derived from input parameters above
Menv (M⊙)	Mass of the envelope	1.77	1.04
cs (km s−1)	Thermal sound speed using ${\dot{M}}_{\mathrm{env}}=0.975{c}_{s}^{3}/G$	0.23	0.27
Rcol (au)	Inside-out collapse radius using Rcol = cs tage	3800	n/a
L* (L⊙)	Stellar luminosity	0.31	:
Lacc,star (L⊙)	Stellar hot-spot accretion luminosity	2.14	:
Lacc,disk (L⊙)	Disk accretion luminosity	0.29	:
Lint (L⊙)	Internal luminosity	2.74	:

Note. The symbol : means the UCM values are the same as the TSC values.

Download table as:  ASCIITypeset image

Zoom In Zoom Out Reset image size

We adopt a dust model that is appropriate for protostellar environments, based on a comparison of models with observations (Huard & Terebey 2017). Specifically we adopt the thinly ice-mantled, coagulated dust of Ossenkopf & Henning (1994) (see their Table 1, fifth column), often referred to as "OH5" grains in the literature (e.g., Evans et al. 2001; Shirley et al. 2005), augmented by the opacities of Pollack et al. (1994) at wavelengths shorter than 1.25 μm, as described in Dunham et al. (2010). In the next subsections we explain how the main components—envelope (Section 2.2.1), disk (Section 2.2.2), and outflow cavities (Section 2.2.3)—are structured in the model and the physical parameters (Sections 2.2.4 and 2.2.5) we chose to model L1527.

2.2.1. Envelope Prescription

2.2.2. Disk–Envelope Prescription

The boundary between the envelope and the disk is not a well explored area, although it could be important during the Class 0/I phase if infalling material goes through a shock when it impacts the disk. In this first paper, we do not include shocks; however, as a first step toward including shocks we modified slightly the shape in which the disk–envelope boundary is defined in Hochunk3d by including a ram pressure boundary condition for motion perpendicular to the disk midplane. The assumption is akin to that of infalling material hitting a "brick wall" in the disk midplane. If the pressure is expressed in terms of the thermal sound speed as $P=\rho {c}_{s}^{2}$, and the ram pressure as P + ρ v², then the disk–envelope boundary is defined by matching the disk and envelope ram pressures:

$$\begin{eqnarray}&&{\rho }_{\mathrm{disk}}{c}_{{s}_{\mathrm{disk}}}^{2}={\rho }_{\mathrm{env}}({c}_{{s}_{\mathrm{env}}}^{2}+{v}_{{\perp }_{\mathrm{env}}}^{2}),\end{eqnarray} \tag{ 1 }$$

where v_⊥ represents the velocity component that is perpendicular to the disk midplane. The thermal sound speed is calculated assuming ${c}_{s}={({kT}/\mu {m}_{H})}^{1/2}$ and assuming that the dust and gas temperatures are equal, which is a reasonable assumption for protostellar densities. The disk shape that is found using this boundary condition is no longer that of a flared disk, but that of a disk with fairly constant opening angle out to its edge, as described in Section 4.

2.2.3. Outflow Cavities and Outflow Jet Prescription

2.2.4. Parameters

2.2.5. Luminosity

The internal luminosity is fixed at L_int = 2.74 L_⊙ to match the internal luminosity of 2.74 L_⊙ found by Tobin et al. (2008). The internal luminosity includes contributions from the star L_*, plus the accretion luminosity from material falling onto the star L_acc,star, and the accretion luminosity generated within the disk L_acc,disk, as described in Whitney et al. (2013). The mass of a Class 0/I protostar should be small and most of its luminosity should be due to accretion. This led to choosing a stellar radius of 1.70 R_⊙ with an effective temperature of 3300 K, which gives a stellar luminosity of 0.31 L_⊙ (see Table 1). One additional term, L_ISRF, contributes to the total luminosity of the system, L_tot = L_int + L_ISRF = 3.23 L_⊙; the term L_ISRF is due to external illumination by the interstellar radiation field (ISRF) and is based on the Galactic value computed near our solar system; see the description in Huard & Terebey (2017).

The largest luminosity term is due to the stellar hot-spot accretion luminosity, which is defined as

$$\begin{eqnarray}&&{L}_{\mathrm{acc},\mathrm{star}}={{GM}}_{* }{\dot{M}}_{\mathrm{disk}}\left(\displaystyle \frac{1}{{R}_{* }}-\displaystyle \frac{1}{{R}_{\mathrm{trunc}}}\right)\end{eqnarray} \tag{ 2 }$$

where R_trunc is the truncation radius where the disk is truncated by the stellar magnetospheric field, and its value is approximately the same as in Tobin et al. (2008). The disk accretion rate ${\dot{M}}_{\mathrm{disk}}$ is an important parameter that leads to the stellar accretion luminosity, which encompasses about 66% of the total luminosity of the system. The code determines that the value ${\dot{M}}_{\mathrm{disk}}=6.6\times {10}^{-7}$ M_⊙ yr⁻¹ leads to an accretion luminosity onto the star of 2.14 L_⊙.

We also include the disk accretion luminosity in terms of the energy dissipated at the inner boundary of the disk (see Shakura & Sunyaev 1973; Lynden-Bell & Pringle 1974; Kenyon & Hartmann 1987; Whitney et al. 2003b, for more details). Table 1 summarizes the values of the different luminosity terms.

We model the local chemical evolution to determine the abundances of the molecules that have been observed. The chemical models are time-dependent, and we focus on the results at two epochs, 10⁴ and 10⁵ yr, that span the range of possible ages for L1527.

Our chemical network is taken from the UMIST database, RATE12 (McElroy et al. 2013). The reactions of the carbon and oxygen isotopes have been added such that each reaction involving an atom of the major isotopes will have an equivalent reaction involving the minor isotopes (see Willacy & Woods 2009, for more details). Fractionation of the oxygen and carbon isotopes can occur via the reactions listed in Table 2. For reactions involving ¹⁷O the values of ΔE are taken to be the same as for the equivalent reaction of ¹⁸O, and the pre-exponential part of the rate calculation is scaled by the reduced mass (Young 2007).

Table 2. Isotopic Fractionation Reactions Used in the Model

13C+	+	CO	$\rightleftharpoons$	C+	+	13CO	ΔE = 35 K
13C+	+	C18O	$\rightleftharpoons$	C+	+	13C18O	ΔE = 36 K
HCO+	+	13CO	$\rightleftharpoons$	CO	+	H13CO+	ΔE = 9 K
HCO+	+	C18O	$\rightleftharpoons$	HC18O	+	CO	ΔE = 14 K
HCO+	+	13C18O	$\rightleftharpoons$	H13C18O+	+	CO	ΔE = 22 K
H13CO+	+	C18O	$\rightleftharpoons$	HC18O+	+	13CO	ΔE = 5 K
H13CO+	+	13C18O	$\rightleftharpoons$	H13C18O+	+	13CO	ΔE = 13 K
HC18O+	+	13C18O	$\rightleftharpoons$	H13C18O+	+	C18O

Note. ΔE values are taken from Langer et al. (1984). ΔE values for reactions involving ¹⁷O are assumed to be the same as for the equivalent reaction of ¹⁸O, with the pre-exponential factor in the rate calculation scaled by the reduced mass (Young 2007).

Download table as:  ASCIITypeset image

The network also includes gas–grain reactions, i.e., freezeout on the grain (Hasegawa & Herbst 1993), thermal desorption (Hasegawa et al. 1992), as well as photodesorption, and desorption by heating of grains by cosmic rays (Öberg et al. 2009a, 2009b). The freezeout and desorption reactions are also described in Woods & Willacy (2009). Freezeout is assumed to occur with a sticking coefficient of 1.0 (Bisschop et al. 2006) for all species. For desorption processes, the binding energies required are taken from UMIST12. Cosmic-ray heating rates are given by

$$\begin{eqnarray}&&{k}_{{\rm{crh}}}=3.16\times {10}^{-19}\times \exp (-{E}_{b}/70\,{\rm{K}})\end{eqnarray} \tag{ 3 }$$

(Hasegawa & Herbst 1993), where E_b is the binding energy of the accreted molecule. Photodesorption rates are given by

$$\begin{eqnarray}&&{k}_{\mathrm{photd}}={FY}\langle \pi {a}^{2}{n}_{g}\rangle {{\rm{\Theta }}}_{X}\end{eqnarray} \tag{ 4 }$$

(Willacy & Langer 2000) where F is the UV field (the total of the stellar, interstellar, and cosmic-ray-induced fields in units of G₀), Y is the yield per photon, which is taken to be 10⁻³ for all species (Westley et al. 1995), except for oxygen atoms (Y = 10⁻⁴) and H₂O (Y = 10⁻³) and Y = 2 × 10⁻³ for desorption as OH (Hollenbach et al. 2009). The dust is assumed to be well mixed with the gas, n_g is the number density of dust grains (n_g = 10⁻¹² n_H ), and the average 〈π a² n_g 〉 = 2.1 × 10⁻²¹ n_H (standard interstellar value). Θ_X is the surface coverage of species X (= n_s (X)/Σn_s (Y), where n_s (X) is the abundance of X in the ices and Σn_s (Y) is the total abundance of ices). The stellar UV field is assumed to have a typical T Tauri value of 500 G₀ at 100 au (unextincted) from the star, where G₀ is the standard ISRF (Bergin et al. 2003). The dense envelope generates significant extinction. To account for this, the local UV field is decreased to take into account the extinction calculated along the line of sight to the star.

Grain surface reactions are included using the approach of Garrod & Pauly 2011. Only atoms, H₂, and simple hydrides (OH, CH, NH, and their isotopologues) are assumed to be mobile on the grain surface.

Initially we assume that all elements are in their atomic form, except for carbon which is ionic and hydrogen which is 95% molecular. The initial abundances used are given in Table 3. We assume a cosmic-ray ionization rate of ζ_CR = 1.3 × 10⁻¹⁷ s⁻¹.

For the photodissociation of CO (and its isotopologues) and H₂ we use the self-shielding coefficients provided by Visser et al. (2009) assuming a Doppler width of 0.3 km s⁻¹ and isotopic ratios of ¹²C/¹³C = 89, ¹⁶O/¹⁸O = 498, and ¹⁶O/¹⁷O = 1988, which are taken to be the same as values in the local interstellar medium (Wouterloot et al. 2008). Abundances here are prescribed as X, the fractional abundances relative to total hydrogen, ${n}_{X}/({n}_{H}+2{n}_{{H}_{2}})$.

Starting from t = 0, the molecular abundances grow rapidly until t = 10⁴ yr, when the abundance of ¹²CO reaches its maximum achievable value of 7.22 × 10⁻⁵, a value that is set by the assumed carbon abundance (Table 3) and is based on Taurus observations. Therefore we chose t = 10⁴ yr as our starting reference time for the chemistry. The grid in the chemical model follows the same structure as in the physical model (Section 2.2) but, in order to speed up the computation, the chemistry was only computed on every fifth grid point in the spatial grid. The solution of the chemical reaction network reaches good convergence throughout most of the spatial grid, including the disk, envelope, and outflow shell. However, we exclude from consideration the low-density (n = 4000 cm⁻³) outflow cavity due to reduced convergence in this mostly atomic region. There is little impact on our study because the outflow cavity contributes little to the molecular emission that is the focus of this investigation. However, we do capture the chemistry that happens based on our prescription of the outflow shell (see Figure 9) when we add constant velocity in this region between the envelope and outflow; see Section 2.4 for more details.

Step 3 of RadChemT calculates the protostellar environment in a PV cube. Because Hochunk3d did not include spectral line emission we selected RADMC-3D, version 0.41 (Dullemond et al. 2012), to generate synthetic spectral line emission assuming LTE. LTE improves computational speed and holds in locations where the gas density is above the critical density. The generation of synthetic line images takes into consideration the density, temperature, and velocity grids from Step 1, and the abundance grid from Step 2 at t = 10⁴ yr and t = 10⁵ yr time stamps, all of which are translated into the RADMC-3D file format.

Finally, for each model density distribution and epoch, we carry out line-of-sight radiative transfer calculations using RADMC-3D to construct synthetic line and continuum observations, and compare these against the multi-telescope data set for L1527.

Additional inputs for L1527 in RADMC-3D were source inclination i = 85°, distance d = 140 pc, and system velocity v_sys = 6.0 km s⁻¹. For consistency the OH5 dust opacity law is the same as that used for the Step 1 continuum radiative transfer. Standard and reasonable assumptions for the protostellar environment are T_dust = T_gas and 100 for the gas-to-dust ratio. The microturbulent velocity was fixed at 0.1 km s⁻¹, a value that is required for LTE and results in smoothing the line profile. A larger microturbulent velocity value could affect the hyperfine structures of the molecular spectrum by causing line blending. However, data from CARMA CO observations indicate that somewhat larger values for the microturbulence should be used for the outflow shell region. Based on the ¹²CO data presented in Section 4.5.1 and in Appendix C, we chose a very simple outflow velocity prescription having a constant outward radial motion of 3 km s⁻¹ in the outflow cavity and also in the outflow shell. However, a careful treatment of the outflow shell lies outside the scope of this paper. From the inputs RADMC-3D creates a synthetic model PV cube in FITS format that has units of Jy pixel⁻¹.

Comparison with observations also requires applying telescope-specific parameters. The velocity channel width and spatial pixel size were specified as inputs to RADMC-3D. The synthetic images are sampled in real space in order to be able to see more clearly the variations in the chemical abundances, and therefore the dynamical range of the system, which in our models are greater than the interferometric data. To compare with millimeter interferometer data, the synthetic model images are multiplied by a Gaussian primary beam (peak value normalized to unity), and then convolved with a circular Gaussian synthesized beam (area normalized to unity). See Section 3 for ALMA and CARMA observational parameters.

L1527 was observed in 2009 August and November using CARMA, which is located at an altitude of 2200 m on the Eastern California Inyo mountains. Observations were obtained in D- and C-array configurations, which provide a uniform uv-coverage between 9 and 371 m. The CARMA correlator was set to observe the ¹²CO (1–0) emission line (ν = 115.271 GHz) in the upper side band, and the ¹³CO (1–0) (ν = 110.201 GHz) line in the lower side band. These two lines were observed at a velocity resolution of 0.34 km s⁻¹ in two 8 MHz spectral windows. The dust continuum emission was observed in two 1 GHz windows separated by 3.6 GHz and centered at the mean frequency of 112.73625 GHz (λ = 2.66 mm). Here we do not use ¹³CO (1–0) since it has higher optical depth than C¹⁸O, making it more difficult to see the emission from the disk. The bandpass shape was calibrated by observing 3C 84, and the flux calibration was set by observing Uranus. The quasar 3C 11 was observed every 12 minutes to correct for atmospheric and instrumental effects. The data reduction and the image reconstruction were obtained using the MIRIAD software package. For N₂H⁺ data and model comparison, we convolve the model using a geometric mean FWHM of 9'', namely, $b=\sqrt{{b}_{max}\times {b}_{min}}=\sqrt{{10.97}^{{\rm{{\prime} }}{\rm{{\prime} }}}\times {8.73}^{{\rm{{\prime} }}{\rm{{\prime} }}}}={9}^{{\rm{{\prime} }}{\rm{{\prime} }}}$. The observational parameters are listed in Table 4.

Table 4. CARMA Observational Parameters

	12CO (1–0)	N2H+ (1–0)
Target, Date	L1527 Infrared Source, 2009 Aug and Nov
Coordinate center	R.A.(J2000) = 4h39m539000
	Decl.(J2000) = 26°03'10000
Frequency	115.271 GHz	93.17378 GHz
Synthesized beam	332 × 295	1097 × 873
Primary beam	54''	67''
Velocity resolution	0.34 km s−1	0.26 km s−1
Noise level (detected channel)	0.187 Jy beam−1	0.200 Jy beam−1

Download table as:  ASCIITypeset image

The observational data presented for C¹⁸O (2–1), which have a spatial resolution of 0.96'' × 0.73'' (that is the geometric beam) and a maximum recoverable scale of ∼15'', are based on data from the ALMA archive that were taken during cycle 0 on 2012 August 26 (Project code: 2011.0.00210.S). ALMA data with higher spatial resolution exist, but they resolve out much of the ∼10'' emission that is relevant to this study. We convolve the model using the same geometric mean formulation as in CARMA, giving b = 08 for the ALMA spatial resolution. The model is also multiplied by the ALMA primary beam, which was taken to be a smooth 28'' FWHM Gaussian image. Table 5 summarizes observational parameters. More information about the observations and calibration is given in Table 1 in Ohashi et al. (2014).

Table 5. ALMA Observational Parameters

	C18O (2–1)
Target, Date	L1527 Infrared Source, 2012 Aug 26
Coordinate Center	R.A.(J2000) = 4h39m539000
	Decl.(J2000) = 26°03'10000
Frequency	219.5603 GHz
Synthesized beam	096 × 073 (+11°)
Primary beam	286
Velocity resolution	0.17 km s−1
Noise level (detected channel)	8.0 mJy beam−1
Minimum baseline	18 m
Maximum recoverable scale	157
Flux calibrator	Callisto
Gain calibrator	J0510+180

Download table as:  ASCIITypeset image

We present a comparison of the flux density of the simulated L1527 for both UCM and TSC models. The simulated L1527 SEDs are shown in Figure 2 plotted at 85° edge-on (pink), the appropriate inclination of L1527 (Oya et al. 2015). Panels (a) and (b) show the model SEDs computed for a 1000 au (714) aperture and a larger 10,000 au (714) aperture, respectively. We also include a plot to illustrate the strong effect of source inclination on the SED model curves. We re-generate SEDs for L1527 by including the flux density values from Tobin et al. (2008) plus we added flux density values from Herschel (see Tables 6 and 7), which come mostly from the thermal radiation of the envelope. The Herschel data were downloaded from the IRSA/IPAC archive. Since the Herschel CDF spectrum (Green et al. 2016) and the flux density points in the HPPSC catalog (Marton et al. 2017) have apertures of 6–14'', i.e., they lie in between the model 714 and 714 apertures, we chose to show the Herschel data on both SED aperture plots. Although we do not intend to redo the best parameter fit for L1527, here we revisit the outskirts of the envelope, where the emission is also important for the chemistry.

Table 6. Photometry at 714 Aperture Size

Wavelength	Fλ	Aperture a	Reference
(μm)	(mJy)	(arcsec)
2.16	0.594 ± 0.16	7.14	1
3.6	6.936 ± 0.69	7.14	1
4.5	22.75 ± 2.28	7.14	1
5.8	29.93 ± 2.99	7.14	1
8.0	18.83 ± 3.80	7.14	1
24	660.6 ± 66	13	1
70–160 b	22,000–64,000 ± 500–7000	14	2
70 b	16,746.0 ± 49.0	6	3
100 b	28,942.0 ± 653.0	6	3
160 b	47,011.0 ± 16,291.0	12	3
1300	375.0 ± 75.0	6	4
2700	47 ± 5.6	3.2 c	5, 6

Note.

^a Radius. ^b The apertures of the Herschel HPPSC catalog and CDF archive spectrum generally lie between the 714 and 714 aperture values, thus the same Herschel data are shown on both SED plots. ^c Radius of the cited beam; either the FWHM size or geometric mean, divided by two.

References. 1, Tobin et al. (2008); 2, Green et al. (2016); 3, Marton et al. (2017); 4, Motte & Andrè (2001); 5, Ohashi et al. (1997); 6, Terebey et al. (1993).

Download table as:  ASCIITypeset image

Table 7. Photometry at 714 Aperture Size

Wavelength	Fλ	Aperture a	Reference
(μm)	(mJy)	(arcsec)
2.16	35.2 ± 16.2	71.4	1
3.6	141.8 ± 16.2	71.4	1
4.5	225.1 ± 16.3	71.4	1
5.8	149.5 ± 45.0	71.4	1
8.0	54.5 ± 25.0	71.4	1
25	743.6 ± 70.0	23 × 150 c	6
60	17,770.0 ± 1600.0	45 × 150 c	6
70–160 b	22,000–64,000 ± 500–7000	14	2
70 b	16,746.0 ± 49.0	6	3
70	24,170.0 ± 4834.0	75	1
100	73,260.0 ± 11,700.0	90 × 150 c	6
100 b	28,942.0 ± 653.0	6	3
160 b	47,011.0 ± 16,291.0	12	3
160	94,000.0 ± 38,000.0	60	7
350	44,000.0 d ± 20,000.0 d	45–60e	1
450	33,125.0 d ± 20,900.0 d	40–120 e	1
750	8400.0 ± 1100.0	45	8
800	1400.0 ± 560.0	60	7
850	6167.0 d ± 480.0 d	40–120 e	1
1300	1110.0 d ± 110.0 d	30–40 e	1
2700	47.0 ± 5.6	3.2 c	4, 5

Notes.

^a Radius. ^b The apertures of the Herschel HPPSC catalog and CDF archive spectrum generally lie between the 714 and 714 aperture values, thus the same Herschel data are shown on both SED plots. ^c Radius of the cited beam; either the FWHM size or geometric mean, divided by two. ^d For a given wavelength >200 μm, the average of values that are listed in Tobin et al. (2008) and the largest error fluxes is adopted. ^e Range of contributing apertures.

References. 1, Tobin et al. (2008); 2, Green et al. (2016); 3, Marton et al. (2017); 4, Ohashi et al. (1997); 5, Terebey et al. (1993); 6, Beichman et al. (1988); 7, Ladd et al. (1991); 8, Chandler & Richer (2000).

Download table as:  ASCIITypeset image

Zoom In Zoom Out Reset image size

Overall, both the UCM and TSC models provide reasonable SED fits, with the TSC model providing a better fit for the 10,000 au aperture that is consistent with the large spatial extent of the system. Another regime that supports our later statement is the region where the disk emits, which is between 2.16 and 8.0 μm. We find that the TSC is a closer fit to the flux points, although not perfectly due the different dust opacity, population, and density that exist here, which are not necessarily in accordance with the dust opacity model we use for the envelope. The PACS data are also of particular interest because they occur near the peak of the SED. The Herschel PACS HPPSC data have an aperture of 6'' at 70 μm and 100 μm, and 12'' at 160 μm. The PACS data at 70 μm and 100 μm have a spatial resolution that is comparable with that of the 714 aperture model (Figure 2(a)), and moreover, the flux density values do not deviate much from the modeled ones. However, for PACS > 100 μm the resolution is slightly offset since the discrepancy in aperture size is greater. The Herschel spectrum does not match the Herschel photometry points either, so the Herschel data are not consistent with each other at around 100 μm and longer wavelengths. Some possible explanations are that (a) the source is extended, which might mean the data calibration is off/incorrect (data issue), and (b) the source is varying in luminosity (source variability).

The properties of the central region were chosen carefully (see Section 2.2.4) to best describe the physical structure of L1527 as constrained by imaging data and the SED photometry. We estimated these stellar properties properly describing L1527 based on Tobin et al. (2008, 2010, 2013). We also adjusted the description of the outflow cavity to transition from jet to wide-angle outflow based on ALMA CO data (Figure 8) at 75 au (see Section 2.2.4). For the fixed protostellar parameters listed in Table 1 and described in Sections 2.2.4 and 2.2.5, we varied the mass infall rate and disk mass, finding the mass infall rate, ${\dot{M}}_{\mathrm{env}}$, to be 5 × 10⁻⁶ M_⊙ yr⁻¹ and 3 × 10⁻⁶ M_⊙ yr⁻¹ for the UCM and TSC models, respectively. Based on the magnitude of ${\dot{M}}_{\mathrm{env}}$ chosen, the age of the simulated L1527 is estimated to be $\tfrac{0.22\,{M}_{\odot }}{3\times {10}^{-6}\,{M}_{\odot }\,{\mathrm{yr}}^{-1}}$ ∼ 7 × 10⁴ yr. We adopt that L1527 is in the protostar collapse phase with an age of t ≃ 10⁵ yr.

To further assess the self-consistency of the model, we note that in every evolutionary stage of the star formation process, the mass budget of the infalling envelope and accreting disk have to be consistent with the feedback of the outflows and winds. In the absence of any outflow, ${\dot{M}}_{\mathrm{env}}={\dot{M}}_{\mathrm{disk}}$, all of which falls onto the central star. As discussed in Section 2.2.5 and shown by Equation (2), the disk accretion ${\dot{M}}_{\mathrm{disk}}$ leads to generous L_acc,star accretion luminosity, which is fit by the modeling procedure. Therefore a mismatch in the two accretion rates is related to outflow feedback. The infall efficiency is given by the ratio of the accretion rates; for L1527 our model value of ${\dot{M}}_{{\rm{disk}}}/{\dot{M}}_{{\rm{env}}}\,=6.6\times {10}^{-7}\,{M}_{\odot }\,{{\rm{yr}}}^{-1}$/3 × 10⁻⁶ M_⊙ yr⁻¹ = 0.22. We point out that this value is similar to the value of 0.25 estimated for the protostar TMC-1 (Terebey et al. 2006). These values differ from unity, and provide interesting constraints on theoretical discussions of star formation efficiency. See Matzner & McKee (2000), and in the context of high-mass star formation, Zhang et al. (2014) for an extensive treatment of outflow feedback and its relation to star formation efficiency.

Next, for comparison with van 't Hoff et al. (2018), we produced a second UCM run by increasing the stellar mass to 0.45 M_⊙ and the disk radius to 125 au. The untuned model produced SED fits that were worse but still reasonable. No significant difference could be discerned in the density, temperature, or resultant chemistry. However, the velocities increased by a factor of $\sqrt{2}$ due to higher gravity from the larger mass. The effect was clearly seen in model spectral line profiles and PV diagrams (see further in Section 4.5.2), thus making it a promising way to investigate the dynamical mass of L1527 in future studies.

Figure 3 presents the dust temperature (top panels) and gas density (bottom panels), computed using the RadChemT package, as meridional cuts through the envelope + disk + bipolar outflow cavity. Views on three different scales are shown from left to right. The central star is located at the origin, the outflow cavities extend horizontally, and the edge-on disk (vertical) is oriented 90° with respect to the outflow cavities. The dust temperature is calculated from the radiative equilibrium solution using Hochunk3d and its distribution varies as the power r^−0.5 in optically thick regions, and as the power r^−0.33 in optically thin regions (Kenyon et al. 1993). Near the protostar, the temperature reaches the sublimation temperature of ∼1600 K inside the outflow cavities, and the small amount of dust sublimates due to exposure to stellar radiation in this region, making it feasible for the atomic gas to be present at vibrational energy levels. The radius of the collapsing region (i.e., expansion wave) is 3800 au, and grows larger in time at the sound speed. Outside the collapse radius, the distribution of the envelope is a power-law function of the radius, ρ ∼ r⁻². Inside the collapse radius the density distribution becomes flat, transitioning at smaller radius to the freefall zone where the density behaves as ρ ∼ r^−3/2 well outside the disk (Shu 1977; Terebey et al. 1984). Contour lines in white show T = 25 K (top panels), our definition of the evaporation temperature of CO (e.g., Qi et al. 2015, 2019; Wiebe et al. 2019), and number density n = 10⁹ cm⁻³ (bottom panels). The shape of the dual outflow cavity is seen most clearly in the middle density panel in blue; near 75 au the narrow outflow/jet opens into a wide-angle outflow. The circle in the lower right shows the (small) 1000 au radius aperture (714) that is used for aperture photometry. Note that this aperture covers a small portion of L1527, making it necessary to compare PACS data points with 714 (10,000 au), which covers a much larger region.

Zoom In Zoom Out Reset image size

The outer disk radius equals the centrifugal radius of the envelope, 75 au, and the disk–envelope boundary is wedge-shaped due to the ram pressure boundary condition (see Section 2.2). The region of highest density is found at the disk midplane with the number density of the order of ∼10¹⁴ cm⁻³. In this region of high density inside 50 au, the temperature in the disk midplane is slightly higher than 30 K. Beyond 75 au, the midplane is shielded from the stellar radiation, leading molecules to gradually freeze out onto dust grains until T = 25 K at ∼225 au, where the CO snowline lies and freeze-out is expected to happen more rapidly. Beyond the T = 25 K contour line, the temperature decreases until it reaches the typical temperature of the envelope, T = 10 K. Since we do not include shocks, there is no increase in the temperature profile due to them. The luminosity that shocks can produce is a small fraction of the stellar luminosity; however, it is been suggested that can be enough to liberate molecules from grains (Sakai et al. 2014). Comparing the UCM and TSC models, the temperature and the density distributions of both show very little difference in terms of structure.

The timescale for freeze-out onto grains depends on both the temperature and density (Hartmann et al. 2004), which vary throughout the cloud. The chemistry at the outer boundary of the cloud is determined by the low temperature, low density, and exposure to the ISRF. Within the cloud envelope, at first the rapid increase of density inward dominates, so that molecules freeze out onto grains in regions of low temperature and high density. Nearer the protostar the temperature rises above the desorption temperature, at which point the behavior changes, and molecules come off the grains and are released into the gas phase. The chemistry is also affected by the UV radiation field from the protostar. Dust extinction shields regions near the midplane from the protostellar UV field.

The dense gas distribution is often traced using isotopologues of CO, such as C¹⁸O, which are less likely to be optically thick than CO itself. For both time steps, t = 10⁴ yr and t = 10⁵ yr, high C¹⁸O abundances (∼10⁻⁷ relative to hydrogen) primarily come from the disk and envelope at radii smaller than 400 au, where the temperature is ≥25 K (see Figure 3) and the chemistry can be dominated by thermal desorption reactions and perhaps photodesorption reactions due to the capture of stellar radiation. Due to the lower abundance and optical depth of C¹⁸O, the effect of photodissociation is visible as a drop of 10× in C¹⁸O abundance already in the outer envelope (∼8000 au) at t = 10⁵ yr. But in terms of spatial extent, the C¹⁸O is similar to ¹²CO, as seen in Figure 4(a), but different for absolute abundance. We found no significant difference between UCM and TSC models for chemical abundances, and therefore limit the discussion to the TSC model abundances.

Zoom In Zoom Out Reset image size

Our chemistry matches expectations in the envelope outside the edge of the disk, where we expect lower concentrations of C¹⁸O (≤10⁻⁸) due to the high-density area that is shielded from stellar radiation, and thus it is cooler in this region. An abundance gap in the midplane beyond 300 au is present for both t = 10⁴ and t = 10⁵ yr, and is consistent with the depletion expected in the region beyond the CO snowline (T ≤ 25 K) in the disk and envelope midplane. At t = 10⁵ yr there is a steep drop in the C¹⁸O abundance between 400 au and 4000 au due to the short freeze-out timescale compared with the outermost part of the cloud, where the freeze-out timescale becomes longer (i.e., Caselli et al. 1999).

In the outermost region in the envelope (r > 4000 au), the temperature drops to T = 10 K; the light enhancement of C¹⁸O abundance (∼10⁻⁸) in this region is influenced by photodesorption reactions such as cosmic rays coming from outside the cloud.

The chemistry results presented here for ¹²CO, C¹⁸O, and also for N₂H⁺ (in Section 4.4) represent a small subset of molecular abundances that are available for comparison with observations. Our RadChemT model includes a chemical network that provides abundance predictions for 292 chemical species; this number includes carbon and oxygen isotopologues and 90 grain surface abundances, a promising area for further studies.

Snowlines, such as those for water and CO, are important because they can influence the efficiency of planet formation within the cold shielded regions of disks. Because CO line emission is optically thick, which makes it difficult to observe the dense interior regions, recent studies by Aikawa et al. (2015) and van 't Hoff et al. (2017) suggest that it is necessary to derive the location of the CO snowline from N₂H⁺ observations. Observations have established that the N₂H⁺ abundance is anticorrelated with CO abundance in systems ranging from starless cores (i.e., Tafalla et al. 2004), to class 0 protostars (i.e., Jørgensen 2004), to protoplanetary disks (i.e., Qi et al. 2013, 2015, 2019; Wiebe et al. 2019).

Our chemistry results show a general predicted trend, that gas-phase ¹²CO and N₂H⁺ are anticorrelated in the envelope. In the upper middle and right panels of Figure 4(a), at t = 10⁴ yr the ¹²CO is depleted from the gas phase, and at t = 10⁵ yr (lower middle and right panels) the ¹²CO depletion region increases outward. The increase in N₂H⁺ concentration, beyond 400 au, is consistent with the ¹²CO depletion to t = 10⁵ yr (Figure 4(b), lower middle and right). Inside 320 au, ¹²CO is enhanced since T > 25 K in this region (see the first two top panels in Figure 3) due to thermal and photodesorption reactions whereas N₂H⁺ is absent. Therefore, the ¹²CO snowline is present and the best anticorrelation takes place at T ∼ 25 K in the midplane.

The predicted anticorrelation covers a larger region at t = 10⁵ yr, which is the nominal age of L1527, and shows that the N₂H⁺ concentration grows and extends from 400 au to ∼2000 au (close to the midplane). The N₂H⁺ is present in the outer envelope at larger abundances at t = 10⁵ yr than at t = 10⁴ yr. Thus, our models show a general predicted trend, that gas-phase ¹²CO and N₂H⁺ are anticorrelated in the envelope. Although the CO depletion timescale differs somewhat from the N₂H⁺ growth timescale, the increase in N₂H⁺ abundance should closely follow the slow depletion of ¹²CO since gas-phase timescales are much shorter than freeze-out timescales. The chemistry models are thus consistent with N₂H⁺ being a good tracer for the ¹²CO snowline. Note that the abundances in the outflow cavity are excluded from consideration (see Section 2.3). The jagged boundary between outflow cavity and envelope is due to grid sampling effects mentioned in Section 2.3.

In order to compare our models with observations, we generate synthetic model spectral line images of N₂H⁺ for CARMA and compare their anticorrelation with ¹²CO CARMA data (Section 4.5.1). In order to convey the validation of RadChemT, we generate a synthetic spectral line for C¹⁸O, from which we extract synthetic PV diagrams to compare with ALMA observations (Section 4.5.2). Each model spectral line data cube is based on the density, temperature, velocity grids, and chemical abundances analyzed in previous sections.

4.5.1. N₂H⁺ and ¹²CO CARMA

Figure 5 presents the RadChemT model in the form of an integrated intensity map of N₂H⁺ at t = 10⁵ yr. The best N₂H⁺ and ¹²CO anticorrelation takes place at t = 10⁵ yr, where the greater extent of the N₂H⁺ abundance (see Figure 4(a)) is more consistent with the CARMA data than at t = 10⁴ yr. In Figure 5 we predict the spatial extent of the N₂H⁺ intensity. By not convolving the modeled emission of Figure 5 with the beam, we are able to see in more detail the spatial structure of the N₂H⁺ emission in the envelope, otherwise the X-shape emission would partially disappear and look more like a vertical bar, very similar to how it looks in Figure 6(a). As previously discussed in Section 4.4, the predicted N₂H⁺ emission traces cold dense gas in the envelope between 400 au and ∼3000 au, extending north–south along the midplane. The predicted emission is seen to peak north and south of the protostar position, with less emission at ≤400 au, or about 3'' in radius, where Figure 4 predicts that N₂H⁺ should be absent. Although less visually prominent, there is also extended N₂H⁺ that corresponds to the outer envelope, which is faintly seen in Figure 5 at the level of ∼0.002 Jy pixel⁻¹, extending across the simulated image.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

The N₂H⁺ spectrum is constructed by integrating over the entire 10,000 au image, and this includes and confirms the hyperfine components seen in Figure 6(b). The red solid line (model) and green plus symbols (data) show good correspondence at t = 10⁵ yr. Note that the extended (faint) emission contributes significantly to the predicted spectrum. For the two emission peaks from the inner envelope (two red peaks in Figure 5), there is a discrepancy between the spatial extent of our prediction (Figures 4 and 5) and the observation (Figure 6), which we conjecture might be improved by increasing the collapse age of the system in terms of changing the dynamics so that R_col leads to a larger collapsing region. This is a good motivation to follow a time-dependent evolutionary parcel and thus confirm the spatial extent. Overall, the strength of the observed N₂H⁺ emission is consistent with the model prediction as seen in the spectrum.

Zoom In Zoom Out Reset image size

Figure 6 also presents the integrated intensity maps of ¹²CO (1–0) from CARMA that cover the same 714 (10,000 au) field of view and which show outflowing CO gas that extends east–west, perpendicular to the N₂H⁺ emission. The dust continuum emission (red contours) is centered on the protostar position. In ¹²CO a narrow "jet" extends east–west from the protostar, merging into a wide-angle outflow on both sides of the protostar. The narrow "jet" is also confirmed by ALMA (Figure 8), therefore supporting the inclusion of an inner outflow jet in our model. Figure 9 in Appendix C contains the ¹²CO velocity channel maps, which further show that the emission arises in an (±3 km s⁻¹) outflow shell. None of the emission extends north–south, since the high optical depth of ¹²CO blocks any meaningful view of the low-velocity envelope.

We conclude that the N₂H⁺ emission shows evidence for the predicted anticorrelation with ¹²CO. The N₂H⁺ appears to be missing within 400 au of the protostar, just where full-strength ¹²CO emission is predicted. One difference is that our predicted emission is less extended than the data (Figure 5), which extends out to ∼30'' (∼4200 au). However, the data agree in showing N₂H⁺ enhancement at ±18'' ∼2500 au from the protostar, where the model predicts there is severe depletion for ¹²CO (see Figure 4). A factor that might help to further improve the comparison is to increase R_col and follow the evolutionary process of the chemical parcel spatially.

4.5.2. C¹⁸O ALMA