Free Access
Issue
A&A
Volume 635, March 2020
Article Number A188
Number of page(s) 20
Section Interstellar and circumstellar matter
DOI https://doi.org/10.1051/0004-6361/201936877
Published online 02 April 2020

© ESO 2020

1 Introduction

Starless cores are the sites of future star and planet formation. Although we know gravitational collapse plays a main role during star formation, the details of the collapse process, particularly the timescale, are not yet well understood. There are two major dynamical theories with different dominating mechanisms. One theory is the ambipolar diffusion-controlled scenario, advocating that molecular clouds evolve quasi-statically to form starless cores after 1–10 Myr (Shu et al. 1987; Tassis & Mouschovias 2004; Mouschovias et al. 2006; McKee & Ostriker 2007). The other theory is the supersonic turbulence-induced scenario, in which cores can form faster through gravo-turbulent fragmentation with a free-fall timescale less than 1 Myr (Mac Low et al. 1998; Padoan & Nordlund 1999; Hartmann et al. 2001; Federrath & Klessen 2012). Measuring the age of the cores would therefore allow us to select the correct scenario but the task is difficult. Time-dependent chemical analysis has long been thought to be a solution to this problem but unknown initial conditions too often hamper a precise determination of the elapsed time from the molecular cloud assembly to the formation of the core (see, e.g., Nilsson et al. 2000 for a discussion on the sensitivity of a chemical model to initial conditions). To circumvent the initial condition problem, Pagani et al. (2009a, 2013) proposed an approach based on the N2 D+/N2H+ deuteration profile across the starless core. These works showed that the profile itself is mostly dependent on the time it took to form the starless core, showing a clear distinction between fast collapsing cores with steep profiles and slow collapsing cores with much higher deuteration level and much flatter profiles. This is possible thanks to the slow decrease of the ortho/para-H2 (thought to be originally produced on grains in the 3:1 ratio), which is comparable to the lifetime of starless cores and because a high abundance of ortho-H2 prevents thedeuteration to proceed (Pagani et al. 1992a, 2009a, 2011). Varying the initial conditions did not change the clear-cut separation between the two types of profiles. L 183 was shown to clearly belong to the fast collapse timescale (in a free-fall time; Pagani et al. 2013). More recently, other dynamical studies also based on deuteration (but not on N2 H+ deuteration radial profile evolution) have been published (Kong et al. 2015; Lackington et al. 2016; Körtgen et al. 2017, 2018; Giannetti et al. 2019); but these works have contradictory conclusions concerning the typical lifetime of cores.

The depletion effect on heavy species (e.g., CO, CS, SO) makes these species ineffective tracers of nH_2 and Tkin in starless cores. Undepleted high-density tracers with low excitation energies are needed to reveal the inner structure of starless cores. NH3 and N2 H+ line emissions are commonly found to be spatially anti-correlated with CO line emission among low- and high-mass starless cores, which suggests that light nitrogen-bearing species are better tracers of the inner structure of starless cores (e.g., Tafalla et al. 2002; Fontani et al. 2006; Pagani et al. 2007). The N2 H+ (J = 1–0) rotational line provides a better constraint on both and Tkin in dense cores than the NH3 (J, K) = (1,1) inversion line, since the former has a higher critical density (5 × 105 cm−3) compared with the low critical density (2 × 103 cm−3) of the latter (Pagani et al. 2007). However, when density is higher than ~106 cm−3, N2 H+ and its isotopologs start to decrease in abundance, which suggests that their parent molecule, N2, has frozen out onto dust (e.g., Bergin et al. 2002; Pagani et al. 2005, 2007, 2012). In this case, since the CO- (and N2 -) depleted environment enhances the deuterium fractionation, H+, H, and their isotopologs would be the most abundant ions. Then H2 D+ becomes the best tracer for the innermost region in starless cores owing to the low excitation energy (ΔEk = 18 K) of the ortho-H2D+ ground transition, = 110 –111 (Khersonkii et al. 1987; Pagani et al. 1992b, 2009a; Bergin et al. 2002; van der Tak et al. 2005; Caselli et al. 2008).

Many studies have estimated the depletion level in starless cores from column density maps, which imply a uniform density along the line of sight (e.g., Caselli et al. 1999; Bacmann et al. 2002; Bergin et al. 2002; Kong et al. 2016). However, the depletion level at the center of cores could be underestimated because of line-of-sight integration of their undepleted outer regions. Pagani et al. (2007) performed spherically symmetric “onion shell” modeling of N2 H+ data toward L 183, and showed a volume depletion of N2 H+ by a factor of 6, which is larger than the depletion factor of <2 toward the same source estimated from the column density in their earlier work (Pagani et al. 2005). Similarly, Pagani et al. (2012) show that CO depletion reaches a factor of 2000 instead of the lower limit of 43 they found in an earlier work (Pagani et al. 2005).

With the tools we developed (Pagani et al. 2007, 2009a, 2012, 2013) we can now expand our work, solely based on L 183 until now, to other clouds. We have selected a list of ten clouds spread among several star-forming regions (Auriga, Taurus, Serpens, and ρ Ophiuchus) for which we would like to study the dust content and dust growth (via the coreshine effect, Pagani et al. 2010a; Steinacker et al. 2010), the CO and N2 depletion profiles and the N2 D+/N2H+ deuteration profile. Dust growth, depletion profile, and deuteration profile are three relatively independent ways to measure the age of starless cores and therefore reinforce our confidence in the age determination of each of these cores. This allows us to address the discrepancies revealed by the various works on dynamical deuteration mentioned above. Each core is studied in detail. From the comparison of all the cores, we gain a better understanding of starless core formation and beyond, of low-mass star formation. This first paper focuses on the gas phase of L 1512.

L 1512 is a relatively isolated nearby starless core that is located at a distance of 140 pc near the edge of the Taurus–Auriga molecular cloud complex (Myers et al. 1983; Launhardt et al. 2013). At 1 pc scale, L 1512 is surrounded by CO gas elongated in the north-south direction with a sharp edge in the southern region (Falgarone et al. 1998, 2001). The L 1512 envelope is found to have no significant infall but probably oscillatory (or expansion) motion (Lee et al. 2001; Sohn et al. 2007; Lee & Myers 2011). Kim et al. (2016) further propose that the southern sharp edge is caused by external radiation from an A-type main-sequence star located at 97 pc away from us; this means a distance of ~40 pc between the two if L 1512 is at the Taurus–Auriga distance, which is a separation too large for a significant impact. The L 1512 core has a globular shape with a size of ~0.08 pc, and is quiescent with line widths of ~0.2 km s−1 (Cox et al. 1989; Fuller & Myers 1993; Caselli et al. 1995, 2002). The comparison between the HC3 N and NH3 line widths shows a small nonthermal velocity dispersion of 0.038 km s−1 and an averaged gas kinetic temperature of 11.6 K (Fuller & Myers 1993). N2 H+ observations reveal a nondepleted morphology at 0.1 pc scale, but is found to be depleted in the innermost region of L 1512 with an offset between the N2 H+ and dust emission peaks (Lee et al. 2003).

In this paper, we describe our observations in Sect. 2 and present them in Sect. 3. In Sect. 4, we perform the 1D spherically symmetric nonlocal thermal equilibrium (non-LTE) radiative transfer model and time-dependent chemical model. In Sect. 5, we discuss the lifetime and the possible formation mechanism of L 1512. We summarize our results in Sect. 6.

2 Observations and data reduction

Observational parameters are summarized in Table 1 and the pointings of spectral line observations are shown in Fig. 1.

2.1 IRAM 30 m observations

We observed L 1512 using the IRAM 30 m telescope in December 2013, May and October 2014, and September 2017. Five molecular lines were observed in frequency-switching mode using the dual polarization Eight MIxer Receiver (EMIR), the Versatile Spectral Autocorrelator (VESPA), and the Fourier Transform Spectrometer (FTS); these lines are N2 H+ (1–0), N2 H+ (3–2), N2 D+ (2–1), DCO+ (2–1), and DCO+ (3–2). The frequency resolution was adapted to keep the velocity sampling close to 50 m s−1. Both polarizations are systematically averaged. The spectra were observed on a 12′′ grid in an east-west strip 12′′ south of the core center (RA, Dec)J2000 = (5h 04m07. s5, +32°43′25.′′0) and along a north-south strip across the core center. Pointing was regularly checked every 1.5 h and pointing error was monitored to be less than 3′′. Data were subsequently folded and baseline subtracted with CLASS1. The original cut passes 12′′ south of the true center because when the project started the CFHT and SCUBA-II maps were not yet available and the core center position was not accurately known. Complementary observations in September 2017 (Director Discretionary Time) were performed. Two short east-west strips across the core center and 12′′ above the center, and one short strip along the north-south strip (mentioned above) for N2 H+ (1–0) and DCO+ (2–1) allowed us to check that the modeling of the core is still consistent with the included offset. Figure 1 shows the pointings of spectral line observations overlaid with the SCUBA-II 850 μm map in Fig. 2g, and the observational parameters are summarized in Table 1. We complemented our observations with on-the-fly (OTF) maps of N2 H+ (1–0) and C18O (2–1) obtained by Lippok et al. (2013), and multipointing mosaic maps of C18 O (1–0) and 13CO (1–0) obtained by Falgarone et al. (1998). Lippok et al. (2013) and Falgarone et al. (1998) give observational details. The above data are regridded to supplement two orthogonal strips (Figs. 3 and 4) and measure the spectral center velocity (Fig. 5d) using a Gaussian fit routine in CLASS. We also produce the moment maps of N2 H+ (1–0) and C18O (1–0) (Figs. 2h,i and Figs. 5a,b). Among these, the N2 H+ moment maps include our observations and the OTF data.

thumbnail Fig. 1

Multipointing grids and onion-shell model overlaid with the SCUBA-II 850 μm map. The black and yellow dots in a (ΔRA, ΔDec) = (12′′, 12′′)-spacing grid show the pointings of IRAM 30-m and GBT observations, respectively. The red dots in a (ΔRA, ΔDec) = (15′′, 30′′)-spacing grid show those of JCMT HARP observations. The circles at top and bottom indicate the beam sizes of each spectral observations. The gray concentric circles with center of L1512 represent the onion-shell model with shell widths of 12′′. The 850 μm map is shown in blue-scale, and its white and black contours have the same levels as in Fig. 2g.

2.2 James Clerk Maxwell Telescope (JCMT) observations

The N2 H+ (4–3) and H2 D+ (110 –111) observations were carried out in December 2015 using the JCMT 15 m telescope and the 16 pixels HARP receiver (two pixels of which are non-functioning) in frequency–switching mode. Thanks to the proximity of the two lines (250 MHz), they are routinely observed together with two sub-bands from the Auto Correlation Spectral Imaging System (ACSIS) spectrometer. Data are converted to the CLASS format to be reduced with CLASS (folding and baseline subtraction). Observational parameters are summarized in Table 1, and the pointings of spectral line observations are shown in Fig. 1.

SCUBA-II observations were retrieved from the JCMT archive2. They are part of the M13BC01 project (PI: James di Francesco). As far as we know, these data have not been published previously.

2.3 GBTobservations

We observed L 1512 using the GBT 100 m off-axis telescope in November 2014. We used the MM1 (67–74 GHz) and MM2 (73–80 GHz) W-band dual polarization sub-band receivers in in-band frequency-switching mode with the Versatile GBT Astronomical Spectrometer (VEGAS) backend. Pointing was checked every hour and each session started with a surface shape verification by using Out Of Focus Holography (OOF) on a strong source. Data were preprocessed (calibrated, folded, and baseline subtracted) in the GBTIDL data reduction program and converted to CLASS format for subsequent reduction. Thanks to the very smooth baseline of the receiver (due to the use of large band amplifiers instead of mixers), the preprocessing of the data allows us to treat separately the ON and OFF observations (i.e., total power mode) and gain a factor of on the noise (the method will be presented in a subsequent paper).

2.4 CFHT observations

The CFHT Wide InfraRed CAMera (WIRCAM) was used with the wide filters J, H, and Ks to observe the source on the night of 26 December 2013. To preserve the extended emission (scattered light known as cloudshine; Foster & Goodman 2006), we used a SKY-TARGET-SKY nodding mode to subtract the atmospheric infrared emission. Seeing conditions were typically 0.′′8. Data reduction was performed at the (now closed) TERAPIX center using a specific reduction method to recover the extended emission. The integration times were typically of 0.5 to 1 h on-source per filter, to reach a completed magnitude of 21.5 (J band) to 20 (Ks band). The extended emission will be analyzed in a subsequent paper.

2.5 Spitzer observations

Spitzer observations are collected from the Spitzer Heritage Archive (SHA)3. L 1512 was observed in two programs with Spitzer InfraRed Array Camera (IRAC), Program id. 94 (PI: Charles Lawrence) and Program id. 90 109 (PI: Roberta Paladini); the second program occurred during the warm period of the mission (only 3.6 and 4.5 μm channels still working). The P94 data were already discussed in Stutz et al. (2009), while the warm mission data are analyzed in Steinacker et al. (2015). These two papers provide the presentation of the observations. During the warm mission, deep observations of the cloud were performed and a completed magnitude of 18 in band I2 (4.5 μm) was reached.

3 Results

3.1 Continuum maps

Figure 2 shows the continuum maps of L1512 at near-infrared (NIR), mid-infrared (MIR), and submillimeter (submm) wavelengths. The NIR maps of L1512 at J, H, and Ks bands (Figs. 2a–c) reveal extended emission that is distributed from an annular shape at J band and progressively merges to a concentrated shape at Ks band. These structures are caused by the scattered light of ambient interstellar radiation field by dust grains (agr ~ 0.1 μm) at the periphery of L1512. Similar structures in dark clouds in the Perseus molecular complex have been found and studied by Foster & Goodman (2006) and Juvela et al. (2006). Foster & Goodman (2006) introduced this scattering effect as cloudshine, and proposed to use it as a general tracer of the column density of dense clouds. Since the wavelength of the maximum scattering cross section increases with the size of dust grains, the different spatial spread of the cloudshine of L1512 at the three bands also indicates that the dust grains have the smallest size at the outskirts, and become larger toward the inner region as discussed in Steinacker et al. (2010) and Lefèvre et al. (2014, 2016).

On the other hand, the MIR maps of L1512 from IRAC1 and IRAC2 bands (Figs. 2d,e) show a more compact emission feature, while the IRAC4 map (Fig. 2f) shows absorption instead. The absorption feature is mainly caused by the bright PAH emission from the diffuse medium, the extinction of which becomes superior to the inner scattered light contribution (Lefèvre et al. 2016). The similar spatial distributions between the compact emission and absorption indicate that the emission also comes from the central core. Pagani et al. (2010a) and Steinacker et al. (2010) firstly recognized this phenomenon in their sample of nearby starless and Class 0/I protostellar cores. By analogy to cloudshine, they named this effect coreshine. With the benefit of the smaller beam sizes of IRAC maps compared with those of JCMT SCUBA-II maps, we determined the center of L1512 core at (RA, Dec)J2000 = (5h 04m07. s5, +32°43′25.′′0) using IRAC4 band (8 μm) with a beam size of 2′′, and denote it as a cross shown in each panel of Fig. 2.

thumbnail Fig. 2

L1512 maps of continuum and line integrated intensity. The CFHT NIR maps at (a) J band, (b) H band, and (c) Ks band. Spitzer MIR maps at (d) IRAC1 band, (e) IRAC2 band, and (f) IRAC4 band with its contours at 5.15, 5.25, 5.35, 5.45, 5.55 MJy sr−1. (g) JCMT SCUBA-II map at 850 μm with its contours at 0, 0.05, 0.15, 0.25 mJy arcsec−2. Integrated intensity maps of (h) N2H+ J = 1–0 and (i) C18O J = 1–0 calculated within VLSR = [−3 km s−1, 16 km s−1] and [6.0 km s−1, 7.5 km s−1], respectively. The panels h and i are overlaid with contours from the 850 μm image. The central cross in each panel indicates the center of L1512 determined by the IRAC4 map. The scale bars of 0.05 pc and millimeter-wavelength beam sizes are denoted in the bottom right and bottom left corners, respectively.

3.2 Molecular emission lines

Figure 1 shows the pointings of spectral line observations using IRAM 30 m, GBT, and JCMT HARP. Horizontally, our IRAM/GBT observations at ΔDec = −12′′ and our JCMT observation at ΔDec = −16′′ have the widest coverage across the L1512 core. Vertically, our IRAM/GBT observations at ΔRA = 0′′ have the widest coverage, and the JCMT observation closest to it is at ΔRA = −2.′′2. We used these cuts (hereafter referring to them as the main horizontal cuts and the main vertical cuts) to analyze the physical and chemical properties of L1512 (Sect. 4). Figures 3 and 4 show the spectra of each emission line in the scale (except for C18O (1–0) in the TMB scale) along the main horizontal and vertical cuts, respectively.

Since L1512 is a quiescent starless core, each hyperfine component of N2 H+ (1–0) and N2 D+ (1–0) is well separated with line widths of ~0.2 km s−1. Along the main horizontal cuts, the coverage of the N2 H+ (1–0) spectra includes the whole starless core with the intensity peak toward the central region. Likewise, the spectra of N2 H+ (3–2), N2 D+, and DCO+ also peak toward the central region, which suggests that these four cations trace the dense core region. In addition, our highest observed transition of N2 H+, the J = 4–3 line, shows no signal according tothe 3σ criterion with a rms noise level of ~30 mK because of the low kinetic temperature in the core region and low density compared to the requested critical density for this transition (two orders of magnitude lower). Among the four cations, the ortho-H2D+ (110 –111) line is detected in a smaller spatial region (ΔRA = −17.2′′ to 42.8′′), which indicates that ortho-H2D+ traces the innermost region of the L 1512 core. On the other hand, the intensity of the C18 O (1–0) and (2–1) spectra are relatively constant across the inner core, which is a clear indication that C18 O is depleted across the whole core.

Along the main vertical cuts, Fig. 4 shows that N2 H+, N2 D+, and DCO+ are also depleted in the innermost region. Their spectral peaks drop between the region from ΔDec = −24 to +12′′ for which the intensity decrements range from 0.04 to 0.15 K with ~2–5σ significant levels, while they should increase to follow the increase of column density toward the center if there were no depletion. Similarly, the ortho-H2D+ spectra show a roughly constant intensity across the peak, which is a sure sign of depletion since the line is optically thin and should be roughly proportionalto the H2 column density. The depletion effects are also found in the integrated intensity maps of N2 H+ (1–0) and C18O (1–0) shown in Figs. 2h and i. The C18O emission clearly shows depletion toward the center of the cloud because the peak emission is situated on the outskirts, while an offset between N2 H+ emission peak and dust peak suggests that N2 H+ also depletes in the innermost region.

Table 1

Observational parameters.

3.3 Velocity structure

Figures 5a and b show the mean velocity maps of L 1512, which reveals a velocity gradient in the north-south direction. The constant gradient inside the core itself is compatible with a rigid-like rotation feature or could be due to an oscillation of the whole filament. Figure 5d shows the central velocity of each spectral line along the main vertical cuts fitted with the CLASS Gaussian fitting task. The hyperfine structures of the N2 H+, N2 D+, and DCO+ rotational lines are taken into account for the fitting. We see that the N2 H+ (1–0) emission appears to have a rather uniform velocity gradient of 2.26 ± 0.04 km s−1 pc−1 between ΔDec = −48′′ and +36′′, which is compatible with a rigid body rotation. The N2 H+ (3–2), N2 D+ (2–1), DCO+ (2–1), (and possibly ortho-H2D+ (110 –111) although its velocity measurement is hampered by a low signal-to-noise ratio (S/N)) emissions also have similar dynamical behaviors. The C18 O (1–0) and 13 CO (1–0)data from Falgarone et al. (1998) partially follow the ionic species velocity drift but with an offset of ~50 m s−1 toward the center. Because the 13CO emission traces a more extended region than C18 O, the similarity between the C18 O and 13 CO emissions might be due to C18O depletion in the core region and its emission being dominated by the L 1512 envelope. The CO isotopologs redshift by ~50 m s−1 toward the center could be explained by a contraction of the envelope, if the line is optically thick enough to mask the corresponding blueshifted component that we would expect on the back side of the cloud; however, this seems difficult to justify for the C18 O (1–0) transition. The other possible explanation is that the envelope is experiencing an oscillation in which the core is coupled to it with a delay and is presently toppling over. This kind of oscillation has been proposed to explain the Orion and the G350.54+0.69 filaments (Stutz et al. 2018; Liu et al. 2019) and has also been invoked for L 1512 itself (Lee et al. 2001; Sohn et al. 2007; Lee & Myers 2011).

4 Analysis

Our first step is to analyze the dust absorption maps from 1.2 to 4.5 μm to estimate the total column density of gas+dust all over the cloud. We converted this map into a sphere on first-order approximation to analyze the line emissions. We adopted a correspondingly 1D spherically symmetric non-LTE radiative transfer code (Bernes 1979; Pagani et al. 2007) to reproduce our observed emission line spectra. Since the shape of L 1512 is globular, we can assume that its distributions of volumetric number density (), gas kinetic temperature (Tkin), volumetric relative abundances with respect to H2 of the observed species () and turbulent velocity (Vturb) can be approximated by an onion-shell structure composed of multiple concentric homogeneous layers. We note however that the cloud is slightly dissymmetric and we used two slightly different models to describe the east and west sides of the cloud. In addition, rotational velocity field (Vrot) and radial velocity field (Vrad) can also be applied to each layer. We chose the widths of the layers to be 12′′ (=1680 AU), which is the smaller spacing of the two pointing-grids (Fig. 1). The advantage is that we can determine the physical parameters and chemical abundances at each layer sequentially from the outermost to the innermost layer by sampling sightlines of progressively decreasing radius along a cut across L 1512.

Our procedure to determine , Tkin, Xspecies, and Vturb of each layer is as follows. First, the profile was determined independently from an extinction map. Second, we used N2 H+ data to determine their Tkin, X(N2H+), and Vturb profiles. Third, we assumed that all the other species have the same Tkin and Vturb profiles as those of N2 H+ to determine the abundance profiles of N2 D+, DCO+, and ortho-H2D+. Given the abundance profiles of the various ions, we adopted a pseudo time-dependent chemical model (Pagani et al. 2009a) to estimate the lifetime scale of each layer with our derived and Tkin, and we also derived the abundance profiles of their parent species, CO and N2.

thumbnail Fig. 3

Spectral observations along the main horizontal cuts at ΔDec = −12′′ (IRAM/GBT) or −16′′ (JCMT) compared to our best-fit radiative transfer model. The blue spectra show the observational data, and the red spectra show the models. Each column corresponds to different horizontal offsets from the center of L 1512 according to Fig. 1. Each row shows a spectral line, except that the N2 H+ (1–0), N2D+ (1–0) and (2–1) lines are split into three rows corresponding to their different F1 -transition groups. For C18O (2–1), our data are supplemented with OTF spectra from previous works at large offsets. The green dashed lines indicate the three σ noise level.

thumbnail Fig. 4

Spectral observations along the main vertical cuts at ΔRA = 0′′ (IRAM/GBT) or −2.2′′ (JCMT) compared to our best-fit radiative transfer model. For N2H+ (1–0) and C18O (2–1), our data are supplemented with OTF spectra from previous works at large offsets. Other annotations are the same as in Fig. 3.

4.1 Visual extinction

We adopted the Near-Infrared Color Excess Revisited technique (NICER; Lombardi & Alves 2001) to derive the visual extinction, AV, with color excesses of stars from NIR and MIR bands. Since the NIR data preferentially trace the diffuse region, we used the RV = 3.1 model from Weingartner & Draine (2001) to derive theAV map of the envelope (Fig. 6a) using J, H, and Ks bands. On the other hand, the MIR data preferentially trace the dense region. We alternatively used the RV = 5.5 case B model from Weingartner & Draine (2001) to derive AV map of the core region (Fig. 6b) using H, Ks, and IRAC2 bands. The case B model is an improved variant of the RV = 5.5 model from Weingartner & Draine (2001), but characterizes the denser region better (by including spherical grains up to 10 μm in size), although not perfectly (Cambrésy et al. 2011; Ascenso et al. 2013). Thanks to the deep sensitivity of the CFHT and Spitzer observations, the density of stars across the cloud core is high enough to convolve the reddening data with a 2D Gaussian profile of half-power beam width (HPBW) as small as 30′′ to construct the maps. The dense core map traces the central region better than the diffuse envelope map but could not properly trace the outer low-AV region, and vice versa. In order to combine these two maps, we calculated the azimuthally-averaged AV profiles of the envelope and core maps shown in Fig. 6d from the 1′ horizontal strip across the center (Figs. 6a,b) to avoid the cometary tail. We could see that the two profiles merge at ~2.5 mag and become similar in the outer region. Therefore, an AV threshold of 2.5 mag is used to combine the central region from the core map and the outer region from the envelope map, and the result is shown in Fig. 6c.

thumbnail Fig. 5

Velocity structure. Mean velocity maps of (a) the F1F = 2,3–1,2 component of N2H+ J = 1–0 and (b) C18O J = 1–0 calculated within VLSR = [6.75 km s−1, 7.50 km s−1] and [6.0 km s−1, 7.5 km s−1], respectively. Contours from the 850 μm image are overlaid. The rest of the symbols in the panels a and b have the same meaning as in Fig. 2. (c) Model vertical rotation velocities, Vrot, against radius, Δr, adopted in the radiative transfer model. (d) Central velocities and their error bars of the spectra along the main vertical cuts. Two y-axes are vertical offsets to the center and declination, respectively.

4.2 Density profile

The volume density () profile of L 1512 is modeled with the spherically symmetric Plummer-like profile, (1)

where n0 is the central density, R0 is the characteristic radius, η is the power-law index of profile as rR0, and Redge is set as 320′′ (=44 800 AU) to cover the whole core. Since the contribution from the diffuse envelope is small compared with the dense core, we simply adopted the RV = 5.5B model for the conversion of the extinction, AV, into column density, . Therefore we obtain a relation between AV and (Bohlin et al. 1978), (2)

In order to derive n0, R0 and η in the volume density model of L 1512, we produced an AV model map by (1) integrating a Plummer-like model along the line of sight and convolving it with the same beam size as the combined AV map (Fig. 6c) to obtain an model map, and then (2) derive an AV model map using Eq. (2). Therefore, we could perform a fitting on the azimuthally-averaged profiles of the AV model map and the combined AV map with the Levenberg-Marquardt algorithm included in the Python scipy.optimize.curve_fit function.

Figure 7 shows the fitting results for the east and west sides of L 1512 using the AV computed in the strip shown in Fig. 6c, which spatially covers the main horizontal cuts. In order to ensure that the volume density models of both sides agree with each other at the center, n0 was previously fitted by an averaged profile from both sides. Then n0 was fixed inthe individual fitting for each side. These determined parameters and volume density profiles are shown in the first row in Fig. 8. We see that the horizontal profile is skewed in such a way that it is steeper on the west side than the other side. The volume densities drop to ~104 cm−3 at both edges of our spectral observation coverage.

thumbnail Fig. 6

Visualextinction maps and profiles. (a) Envelope AV map derived from J, H, and Ks bands with the RV = 3.1 model. (b) Core AV map derived from H, Ks, and IRAC2 bands with the RV = 5.5 B model. (c) Combined AV map merged from the previous maps with boundaries at AV = 2.5 mag. Each map is convolved with a Gaussian beam size of 30′′. (d) Comparison between AV profiles inside the strips of panel a and b. The AV values at each pixel are denoted in light blue and pink dots, while their average errors are indicated in the bottom left corner. The AV profiles are averaged in azimuth with radial 10′′ bin and shown in blue triangles and red squares with their error bars. The green dashed lines in the four planes indicate the AV threshold of 2.5 mag.

4.3 Radiative transfer applied to the onion-shell model

We assumed an onion-shell model comprised of nine and six 12′′ layers in thickness on the east and west sides, respectively. Then we reproduced the spectra of N2 H+ (1–0, 3–2, and 4–3) along the main horizontal cuts with the 1D spherically symmetric non-LTE radiative transfer code (Pagani et al. 2007) by varying Tkin, X(N2H+), Vturb, Vrot, and Vrad at each layer. Here we used the para-H2–N2H+ collisional coefficients instead of the original He–N2H+ hyperfine collisional coefficients (Lique et al. 2015). We also used specific para-H2–N2D+ collisional rate coefficients for the deuterated isotopolog (see Appendix B). One of the important features of this code is that the line overlap between hyperfine transitions at close frequencies is taken into account to calculate the correct spectral line shape and excitation. Among the input parameters, Vturb, Vrot, and Vrad are determined from observations. We set both Vrot and Vrad to zero, since the L 1512 core has no obvious infall/expansion motion and only a vertical rotational component, while our model follows the RA (horizontal) structure. For Vturb, we found that a uniform turbulent velocity of 80 m s−1 could reproduce the spectral line widths along the cuts. Afterward, we performed an iterative spectral fitting procedure to obtain Tkin and X(N2H+) profiles, taking care to eliminate nonphysical solutions such as parameter oscillations between layers, even if they provide better χ2.

An important feature of our approach is that instead of converting the data into the TMB scale for direct comparison with the non-LTE model, which is absolutely not valid for extended regions such as starless clouds, we reproduce the observations by computing the telescope coupling to the sky, based on its main beam, error beams, and forward scattering and spillover efficiency (ηfss)4. We then obtained our model predictions directly calibrated in the scale, whereas the TMB scale is overcorrected for low ηMB values such as those found for the IRAM 30 m in the 1.3–0.8 mm range since the error beams pick up a substantial fraction of the signal when the source is extended. This leads to wrong line ratios and overestimates of the temperature and/or density of the cloud and also has an impact on the species abundance.

Once the density, temperature, and kinematic properties of the onion-shell model were constrained by the N2 H+ modeling, we reproduced the spectra of N2 D+ (1–0, and 2–1), DCO+ (1–0, 2–1, and 3–2), and ortho-H2D+ (110 –111) by varying only their abundances and assuming that they share the same physical parameters as N2 H+. Similarly, X(N2D+), X(DCO+), and X(ortho-H2D+) profiles were obtained using the above iteratively spectral fitting procedure. Since we lacked N2 D+ observations at large radii, we set upper limits of X(N2D+) in the outer layers from X(N2H+) multiplied by the minimum observed X(N2D+)/X(N2H+) = 0.05. For the non-detection observations of ortho-H2D+, we also estimated their abundance upper limits.

Figure 8 shows the profiles of Tkin, Vturb, abundances of the above four cations, and the N2 H+ fractionation, X(N2D+)/ X(N2H+). The above best-fit profiles are also numerically listed in Table C.1. The error bar in each layer is determined by the range of quantity, where (3)

We only determined the upper limit of X(ortho-H2D+) toward the innermost layer. A higher S/N of ortho-H2D+ would be needed to derive its lower limit because the volume of the innermost layer is the smallest and its contribution to the emergent intensity is also small.

Finally, the best-fit modeled spectra of the four cations along the main horizontal cuts are shown in red in Fig. 3. The fit with the observed spectra is remarkably good. To see if the observed spectra along the main vertical cut could also be reproduced by our onion-shell model, we simply assumed that the physical parameters and abundance profiles were the same as those of the east side except for Vrot profile that is now not zero. Since the L 1512 core has a vertical motion compatible with a solid-body rotation, we adopted the Vrot profile of N2 H+, N2 D+, and DCO+ shown in Fig. 5c, which averages the central velocities of the vertical spectra (Fig. 5d) between the north and south sides. For ortho-H2D+, we simply adopted the Vrot profile of N2 H+ because of its weak detection. The above Vrot profiles are given at layer interfaces and the rotational velocity inside a layer is linearly interpolated between the two given velocities at the inner and outer sides of that layer. The modeled spectra of the four cations are shown in red in Fig. 4, and are again in good agreement with the observations. Our results suggest that the 1D spherically symmetric assumption is globally valid for L 1512.

thumbnail Fig. 7

Plummer-like density profile fittings. The black dots are the AV values at each pixel in the horizontal strip shown in Fig. 6c but these are separated between the east and west sides. The orange squares with error bars show the averaged AV profiles, which are azimuthally binned with radial 10′′ bin. The blue circles show the best-fitted AV models for the east and west sides. The AV models are derived from spherically symmetrical Plummer-like density profiles. The best-fitted characteristic radii, R0, and the fitting ranges, Rfitting of the AV models are denoted in blue solid and dash-dotted lines, respectively. The coverage of our spectral observations are denoted in black dashed lines.

4.4 Time-dependent chemical model

We adopted a pseudo time-dependent gas-phase chemical model (Pagani et al. 2009a) to reproduce the deuteration ratio, X(N2D+)/X(N2H+), and abundances of N2 H+, N2 D+, DCO+, and ortho-H2D+ in each layer. In the pseudo time-dependent approach, the physical properties (, Tkin) are kept constant while the chemicalabundances (Xspecies) evolve. We considered a chemical network, where (1) different spin states (ortho, para, and meta) of H2 and H, and their isotopologs are all included, and (2) heavy molecules are totally depleted except for CO and N2 which are partially depleted.

The initial ortho-H2/para-H2 ratio is assumed to be their statistical weights of 3:1 and H is formed later via the cosmic ray ionization of H2. After the ortho/para ratio (OPR) of H2 has dropped to values low enough to prevent destruction of H2 D+ by ortho-H2 (OPR ≲ 1% or less; Pagani et al. 1992a, 2009a), the deuterium fractionation of H enhances in the low Tkin and highly depleted environment, forming H2 D+, and subsequently D2 H+ and D. Then H and its isotopologs transfer protons/deuterons to CO and N2 and produceHCO+, DCO+, N2 H+, and N2 D+. HCO+ is not a suitable tracer for starless cores because its low rotational transition emission is optically thick and this ion is present in both the envelope and the core (which is also true for H13 CO+ and HC18O+), while ortho-H2D+, DCO+, N2 H+, and N2 D+ emissions are optically thinner and most importantly are confined to the core (Pagani et al. 2012). Since we evaluated the quasi-instantaneous conversion of CO and N2 into the aforementioned ions, we can assume that their abundances are in equilibrium with these ions locally. Therefore, the gas-phase abundances of these two neutral molecules are free parameters in our chemical model. The other free parameters are the average grain size (agr) and the cosmic ray ionization rate (ζ).

Figure 9 shows the chemical model solutions of X(ortho-H2D+), X(N2D+)/X(N2H+), and X(DCO+) for each layer toward the east and west sides. Typical ζ of 10−17 s−1 was adopted throughout the chemical model. For agr, although Steinacker et al. (2015) find that the maximum spherical grain size can reach ~0.5 μm in this source, the majority of the grains are still small and the larger grains have fluffy surfaces (e.g., Tazaki et al. 2016), therefore limiting the decrease of the total grain surface cross section. The impact on the chemistry and the freeze-out time should be negligible compared to the standard 0.1 μm grain size usually adopted in chemical models. In each layer, X(CO) and X(N2) were adjusted to make the chemical solutions of the four cations fit the observed abundances (or remain below the observed upper limits) simultaneously. We found that the time ranges for all layers meeting the observations are 0.2–2.6 Myr for the east side and 0.2–1.9 Myr for the west side.

Figure 10 shows the derived CO and N2 abundance profiles, while their values are listed in Table C.1. Although the isotopologs of CO are not considered in our chemical network, we can obtain the X(C18O) profile through a 12 C16O/12C18O abundance ratio in the layers in which C18O cannot be measured directly. In the outer region, where radiative transfer modeling is applicable, we find that the C18 O abundance is 1 × 10−7 up to a radius of 3′ (=25 200 AU). Since X(C18O) is limited to 1 × 10−7 toward the sixth eastern layer, where the 12CO abundance reaches its maximum of5 × 10−5, we find an isotopolog ratio of 12 C16O/12C18O = 500. This is compatible with the local 16O/18O ISM ratio of 557 ± 30 (Wilson & Rood 1994). As there is no reason to vary this ratio across the cloud, we obtained the inner C18 O abundance profile between the fifth eastern and the third western layers by applying the same ratio of 500 to X(12CO) profile. On the other hand, the outer 12CO abundance profile is kept as high as 5 × 10−5 because X(12CO) is poorly determined from the chemical model when both X(N2D+)/X(N2H+) and X(ortho-H2D+) are only upper limits. As a result, we were able to reproduce the C18 O spectra via the non-LTE radiative transfer calculation. Owing to the nonsymmetrical velocity shift of C18 O, we adjusted the velocity offset of the modeled spectra individually on each spectrum based on the observations. The impact on the line shape is too small to be detectable. The northern C18 O abundance profile (positive offsets) is different from the east-west profile and had to be decreased by a factor of 2–16 to fit the observations. This suggests a larger CO depletion toward the northern tail. It turns out that the observations can be fitted by our model results of C18 O (1–0) and (2–1) along the main horizontal cut (Fig. 3) and main vertical cut (southward only, Fig. 4) computed with the abundances derived from the chemical model, thereby validating our chemical model results.

thumbnail Fig. 8

Physical and abundance profiles. The Plummer-like density () profiles are derived from the combined extinction map (Fig. 6c), where the best-fit parameters are annotated. The profiles of kinetic temperature (Tkin), turbulent velocity (Vturb), abundances (Xspecies), and the abundance ratio of N2D+/N2H+ are the best-fit results from the non-LTE radiative transfer calculations.

thumbnail Fig. 9

Chemical modeling of the abundance ratio of N2D+/N2H+ and the abundances of ortho-H2D+ and DCO+ for each layer. Left and right columns: chemical solutions. Middle column: observationally derived profiles from Fig. 8. Each color represents a layer. In the chemical solution plots, the horizontal solid lines are the observationally derived X(N2D+)/X(N2H+), X(ortho-H2D+), and X(DCO+). The growing curves are the chemical solutions calculated with an average grain size (agr) of 0.1 μm and a cosmicray ionization rate (ζ) of 10−17 s−1. The two black vertical lines in each panel indicate a time range for which the model values cross the observations within their error bars and where the model curves become thicker. Thick dashed lines denote the model range that is both lower than the observed upper limit for that layer and within the global time range of the model. The red dashed vertical lines indicate the lower limit on the lifetime scale of L1512.

thumbnail Fig. 10

12CO, N2, and C18 O abundance profiles. The 12CO and N2 abundance profiles are the best-fit results from the chemical modeling. The C18 O abundance profile is obtained by assuming a 12CO/C18O abundance ratio of 500 or a constant value of 1.0 × 10−7 extending to 25.3 kAU.

5 Discussion

As shown inSect. 4, we built an onion-shell model of L 1512 to represent its physical structure (, Tkin, and Vturb are shown in Fig. 8; Vrot in Fig. 5c) and chemical abundances (N2 H+, N2 D+, DCO+, and ortho-H2D+ are shown in Fig. 8; CO, C18O, and N2 in Fig. 10). Our model reproduces all the observed spectra along the main horizontal and vertical cuts (Figs. 3 and 4) via the non-LTE radiative transfer calculations. In this section, we first compare our findings with other studies and then address the lifetime of L 1512.

5.1 Density and kinetic temperature

By integrating the Plummer profile to the edge of the cloud where it merges in extinction with the background (~320′′ from center), we find a peak of 2.3 × 1022 cm−2 toward the center of L 1512, which is higher than previously reported (0.8−1.3 × 1022 cm−2; Lee et al. 2003; Lippok et al. 2013, 2016). From the non-LTE modeling of both C18 O and N2 H+ we also find 2.1 × 1022 cm−2 in a 180′′ radius, which is a lower limit since the cloud is truncated to the C18O emission extension and is consistent with the Plummer profile. This discrepancy is due to our different approach to derive the total column density. We use dust extinction, which is independent of the dust temperature (Tdust), and N2 H+ modeling for which the uncertainty on the collisional coefficients is extremely small compared to the uncertainty on dust properties in emission at long wavelengths. On top of dust properties uncertainties, it is well known that the problem of quantifying dust in emission is degenerate. This is because of the unconstrained temperature variation along the line of sight that can hardly be disentangled from dust property variations and the nonlinearity of the blackbody intensity with temperature for the combination of wavelengths and temperatures of concern in starless cores (Pagani et al. 2015). Indeed, Lee et al. (2003) and Lippok et al. (2013, 2016) performed their estimations with dust thermal emission by assuming a constant Tdust of 10 K or using a decreasing Tdust toward the core center. However, their are less than ~50% of our result. Pagani et al. (2015) suggest that if only dust continuum data are used, Tdust can be overestimated and consequently the column density of starless cores can be underestimated. This is because warmer dust emission can easily dominate over cold dust emission in the spectral energy distribution fitting.

For the kinetic temperature, we find a central Tkin of 7.5 ± 1 K on the west side, which is better constrained that the central temperature for the east side (7.5 K). We also note that the second central layer on the east side is better constrained too ( K). Thus, we conclude that Tkin is down to 7.5 ± 1 K in the center and rises to ~10–17 K in the outer region. Because dust-gas coupling is efficient enough for cm−3 (Goldsmith 2001), we can assume that Tdust = Tkin in the core center based on ourmodeling. The central Tdust is therefore lower than Lippok et al. (2013, 2016)’s 9.8–11.3 K, which explains our peak column density discrepancy as already noticed in other sources by Pagani et al. (2015).

5.2 Cation abundance profiles

We find that N2 H+ is significantly depleted at the center of L 1512 ( cm−3) and its abundance starts to decrease at of 3 (east) to8 (west) × 104 cm−3. This is one order of magnitude lower than toward L 183. By averaging the depletion factors calculated for the east and west sides, we obtain an averaged N2 H+ depletion factor of 27, which is larger than the factor of 6 in L 183 previously found by Pagani et al. (2007) while L 183 even has a denser central cm−3. Strangely, Lee et al. (2003) suggests that N2 H+ may be depleted in L 1512 but not in another advanced and denser starless core, L 1544 ( cm−3) whereas Redaelli et al. (2019) present new observations and analysis of L 1544 for which they report a large volumetric depletion factor (>100). For N2 D+, we find a depletion factor of 4, which is also higher than the factor of 2–2.5 found in L 183 (Pagani et al. 2007) but lower than the one reported for L 1544 by Redaelli et al. (2019)5, ~15. It is surprising that the depletion factors do not correlate with the peak density among these three sources. One of the possibilities is that the dust grain size in L 1512 might be smaller than in the other cores keeping the depletion timescale shorter (Lee et al. 2003). Another possibility is the time it took to form the cores. L 183 seems to have contracted more rapidly than the two other cores (≤0.7 Myr for L 183, Pagani et al. 2013; ~1 Myr for L 1544, Redaelli et al. 2019; ≥1.4 Myr for L1512, see below Sect. 5.4), though this would have given more time for grain growth in L 1512. Many coupled dynamical and chemical models also cannot simulate such high depletion factor at cm−3 (e.g., Aikawa et al. 2001; Li et al. 2002; Pagani et al. 2013), which might be caused by their too simplistic approach. Consequently, more detailed modeling would be necessary to understand the physical and chemical evolution of L 1512.

From N2 H+ and N2 D+, we can obtain the deuteration ratio, X(N2D+)/X(N2H+), which is increasing from an upper limit of ≤0.05 at large radii to the maximum of 0.34 at the center of L 1512. The higher deuteration ratio at the center seems contradictory with the decreasing abundance of ortho-H2D+ toward the center (by a factor of ≲10). It can be explained if the ortho-H2D+ depletion is due to the further deuteration of H. Namely, H2 D+ is converted to D2 H+ and D, which is supported by our chemical modeling best-fit result, X(D)/X(ortho-H2D+) ≈ 64 and X(D2H+)/X(ortho-H2D+) ≈ 15 at the innermost layer (Fig. 11). The contrast is higher than toward L 183 (Pagani et al. 2009a, 2013) where D2 H+ is comparable to H2 D+ in abundance and D less than a factor of 10 higher than H2 D+. Moreover, in L 183, the ortho-H2D+ does not show depletion toward the center. This can be linked to the longer timescale of evolution of L 1512 as discussed in Pagani et al. (2013) Appendix B, where it is shown that the para-D2H+/ortho-H2D+ is much higher in the slow collapsing case than in the fast collapsing case (the reader should be aware that para-D2H+ is only a small fraction of total D2 H+ when comparing the two papers).

Remaining gas-phase CO can react with deuterated H and form DCO+ (Pagani et al. 2011), we find that DCO+ happens to deplete at of 3 (east) to 7 (west) × 104 cm−3, and reaches a depletion factor of 9 at the center, which is comparable to the factor of ≥ 17 found in L 183 (Pagani et al. 2012). Like in L 183, DCO+ depletes because its parent molecule, CO, freezes out onto dust grains much more than what the growth of H deuteration can compensate for Pagani et al. (2012). Likewise, N2 H+ and N2 D+ deplete because of the freeze-out of N2. CO and N2 depletion are discussed hereunder.

thumbnail Fig. 11

H isotopolog abundances derived from the chemical model constrained by the observations (in particular of ortho-H2D+, which explains why we separately show the two spin states for this species for the sake of clarity).

5.3 CO and N2 depletion

It is impossible to measure X(CO) and X(N2) directly in a dense core. This is because the emission of the low-J transitions of CO (including its isotopologs) usually becomes optically thick even for rarer isotopologs such as C18 O and C17 O, for which the emission is dominated by the envelope and N2 has no transition detectable6 in starless cores (Pagani et al. 2012).

In Sect. 4.4, we presented a chemical model to relate the four observed cations (N2 H+, N2 D+, DCO+, and ortho-H2D+) to CO and N2 in the depleted region to retrieve their abundance profiles. Figure 10 shows the abundance profiles of 12CO, C18 O, and N2. We find that N2 has a depletion factor of ~200 compared with its maximum values of 3 × 10−5. For CO isotopologs, 12CO has a depletion factor of ~220 compared with its maximum values of 5 × 10−5. Compared with the standard 12CO abundance of 1–2 × 10−4 (Pineda et al. 2010), 12CO has a depletion factor of ~430–870. C18 O would also have a depletion factor of ~220 compared with its maximum value of 1 × 10−7 in the outer region since we fix the 12 C16O/12C18O ratio of 500, and a factor7 of ~370 compared to the standard C18O abundance of 1.7 × 10−7 (Frerking et al. 1982), which are much larger than the factors of 30 and 25 found by Lippok et al. (2013)8 and Lee et al. (2003), respectively. Although both studies try to recover volumetric abundances of C18 O and not simply column-density averages, the abundance drop they derive can only be a lower limit because the line from the center is too weak to constrain the models. A factor of 30 already indicates that the contribution of the central part to the total intensity is less than ~10%, which would require a very high S/N to evaluate and still does not guarantee that the measured variation is not simply due to irregularities in the cloud shape rather than abundance variations (see Pagani et al. 2010b for a similar case).

The abundance of CO and N2 in starless cores is a matter of long debate (see, e.g., Pagani et al. 2012; Nguyen et al. 2018) since the daughter species of N2, N2 H+, is present at much higher densities than CO itself in starless cores, while CO and N2 have identical masses and similar sticking coefficients (Bisschop et al. 2006) and their freeze-out rates are thus comparable. For L 1512, 12 CO and N2 profiles are similar to each other within their error bars. This is different from the results reported by Pagani et al. (2012) who found a different N2 profile in L 183, where X(N2) is less than X(12CO) by about two orders of magnitude at the edge of the starless core. These authors suggest that N2 might still be forming from atomic N in L 183. In that case, the depletion and production rates of N2 would partly compensate each other. In contrast with L 183, we find that L 1512 has chemically evolved on a much longer timescale so that N2 chemical production from atomic N has reached steady state. This could explain the similar depletion profiles. The fact that the CO depletion in L 183 is ~3 times deeper than in L 1512 is mostly due to its peak density being 20 times higher, making collisions with dust grains much more frequent, but the difference is partly compensated by a much longer evolution time in L 1512. The L1512 case is interesting because, by showing that N2 and CO profiles are identical while N2 H+ is tracing the core and not CO, it demonstrates that the N2 /CO problem is not correctly formulated. N2 H+ cannot be directly compared to CO because it is a daughter species of N2 and its abundance therefore also depends on the abundance of H. We would need to compare HCO+ to N2 H+, which is not possible in a simple way. Also, N2 H+ is a trace species, typically at the 10−10 abundance level, at least two orders of magnitude less abundant than CO even depleted by a factor 2000 as in L 183. A correct approach is therefore to compare DCO+ and N2 D+; in most clouds, it is clearly visible that DCO+ is more abundant than N2 D+ because its emission is usually as extended as that of N2 H+ itself. There is therefore no N2/CO crisis. However, tracing the abundance of both species in starless cores remains important since it helps the understanding of the chemical evolution of the core.

5.4 Lifetime scale

Figure 9 shows that the chemical solution of each layer meets the observed abundances within ~0.2–2.6 Myr. Since the physical structure is kept constant in our pseudo time-dependent chemical modeling, we see that the chemical process is relatively accelerated in the inner layers compared with the outer layers. Therefore, we underestimate the timescale for the inner layers, which were not very dense at the beginning of contraction. In addition, only the deuteration ratio of N2 D+/N2H+ and the abundance of ortho-H2D+ are used to constrain the chemical age, while X(DCO+) is used to determine the CO abundance. It follows that our result gives a lower limit on the lifetime scale of L 1512 (~1.43 Myr) that is set by the third western layer, which is the last one in our chemical model to reach the observed N2 D+/N2H+ ratio. As a lower limit is required, this value is set at the beginning of the time range of the third western layer as shown on the top right panel (the chemical solutions for the western deuteration ratios) in Fig. 9. We also note that since the third western layer has no ortho-H2D+ detection but only an upper limit on its abundance, the time limit for this layer is the time it takes to reach X(N2D+)/X(N2H+) under the assumption of X(ortho-H2D+) at the maximum of its possible range, that is the fastest time of that layer. Therefore, L 1512 is presumably older than 1.43 Myr. We might ask whether the initial OPR of 3 is warranted. If we set an initial OPR as low as 0.01, we cannot find a good fit for both sides. On the other hand, if we use an OPR of 0.5 as the initial value, as measured in the diffuse medium (Crabtree et al. 2011) toward a few lines of sight, we can only find very marginal fits for both sides and the time remains somewhat above 0.5 Myr. This suggests that the OPR cannot be lower than ~1.

Pagani et al. (2009a) applied the same pseudo time-dependent chemical modeling on L 183 and found the lower limit on L 183 lifetime scale to be 0.15–0.2 Myr. Later, Pagani et al. (2013) conclude that the contraction of L 183 follows the free-fall timescale (~0.5–0.7 Myr). Namely, L 1512 is older than the L 183 starless core at least by a factor of ~2–3. This is in agreement with our conclusion from Sect. 5.3 that L 1512 may be chemically more evolved than L 183, while physically L 183 has reached a higher density in a shorter time. Consequently, ambipolar diffusion may have slowed down the contraction of L 1512 or even halted it to the present state, while it has had no impact on L 183. This would imply that the magnetic field is stronger in L 1512 than in L 183. This will be the subject of a follow-up study of L 1512.

6 Conclusions

We established the extinction map of L 1512 from deep JHKs and Spitzer/IRAC data. We performed a non-LTE radiative transfer calculation with an onion-shell model to reproduce the observed spectra of L1512. We obtained separately, for the east and west sides, 1D spherically symmetric profiles of the physical parameters and chemical abundances of N2H+, N2D+, DCO+, and ortho-H2D+. Afterward, we used a time-dependent chemical model to estimate the lower limit on the lifetime scale of L 1512 by fitting the chemical solutions to the observed abundances. Our chemical model captures the main reactions in the H deuteration fractionation and assumes that heavy species are totally depleted except for CO and N2, which are partially depleted. Thus, by relating the observed abundances to CO and N2 in the depleted region, we could also obtain the CO and N2 profiles and their depletion factors. We summarize our results as follows:

  • 1.

    We derived the central molecular hydrogen density of 1.1 × 105 cm−3 from the dust extinction measured from NIR and MIR maps and the central kinetic temperature of 7.5 ± 1 K from N2H+ observations. At such density, gas is thermalized with dust.

  • 2.

    The depletion factors of N2H+, N2D+, DCO+, and ortho-H2D+ are 27, 4, 9, and ≲10. Compared with N2H+, the smaller N2D+ depletion factor is due to the deuterium fractionation enhanced from a upper limit of ≤0.05 at large radii to 0.34 in the center. The depletion of ortho-H2D+ suggests that ortho-H2D+ is further deuterated to D2H+ and D.

  • 3.

    The depletion of N2H+ and N2D+ traces the freezing out of N2 toward the center of L 1512. We find that N2 has a depletion factor of ~200 compared to its maximum abundance of 4.7 × 10−5. Likewise, the depletion of DCO+ indicatesthe freeze-out of CO. We find that 12CO and C18O have a depletion factor of ~220 betweeninternal and external layers. If we compare their minimum abundance to their standard (literature) abundance, we find the depletion factor is about 400 for both isotopologs.

  • 4.

    The similarity of CO and N2 profiles suggests that L 1512 may have chemically evolved long enough and in particular that N2 chemistry has reached its steady state. This could explain that L 1512 has a higher depletion factor of N2 compared to the denser starless core, L 183.

  • 5.

    L 1512 is presumably older than 1.4 Myr. We conclude that ambipolar diffusion is the dominant core formation mechanism for this source.

In summary, our results present a precise description of density, temperature, and molecular abundances in a starless core.

Acknowledgements

S.J.L. and S.P.L. acknowledge the support from the Ministry of Science and Technology (MOST) of Taiwan with grant MOST 106-2119-M-007-021-MY3 and MOST 106-2911-I-007-504. Nawfel Bouflous and Patrick Hudelot (TERAPIX data center, IAP, Paris, France) are warmly thanked for their help in preparing the CFHT/WIRCAM observation scenario and for performing the data reduction. This work was supported by the Programme National “Physique et Chimie du Milieu Interstellaire” (PCMI) of CNRS/INSU with INC/INP co-funded by CEA and CNES and by Action Fédératrice Astrochimie de l’Observatoire de Paris. This work is based in part on observationscarried out under project numbers 152-13, 039-14, 112-15, and D01-17 with the IRAM 30 m telescope. IRAM is supported by INSU/CNRS (France), MPG (Germany) and IGN (Spain). The James Clerk Maxwell Telescope is operated by the East Asian Observatory on behalf of The National Astronomical Observatory of Japan, Academia Sinica Institute of Astronomy and Astrophysics, the Korea Astronomy and Space Science Institute, the National Astronomical Observatories of China and the Chinese Academy of Sciences (Grant No. XDB09000000), with additional funding support from the Science and Technology Facilities Council of the United Kingdom and participating universities in the United Kingdom and Canada, the JCMT data were collected during programs ID M13BC01 (SCUBA-II), M15BI046 (HARPS, and M17BP043 but data were lost due to an incorrect tuning). This research made use of the Aladin interface, the SIMBAD database, operated at CDS, Strasbourg, France, and the VizieR catalog access tool, CDS, Strasbourg, France. This research has made use of the NASA/IPAC Infrared Science Archive, which is operated by the Jet Propulsion Laboratory, California Institute of Technology, under contract with the National Aeronautics and Space Administration. This work is based in part on observations made with the Spitzer Space Telescope, which is operated by the Jet Propulsion Laboratory, California Institute of Technology under a contract with NASA. This research used the facilities of the Canadian Astronomy Data Centre operated by the National Research Council of Canada with the support of the Canadian Space Agency. Green Bank Observatory is a facility of the National Science Foundation and is operated by Associated Universities, Inc.

Appendix A: Spectral line observations

We present spectral line observations other than the main horizontal and vertical cuts and the full C18 O horizontal cut.

Figures A.1A.5 show all the single pointing observations in black lines and models of N2 H+ (1–0), N2 H+ (3–2), N2 D+ (2–1), DCO+ (2–1), and ortho-H2D+ (110 –111) in red lines, respectively. The models are reproduced with our dissymmetric onion-shell model (Sect. 4 and Table C.1). The full C18O (1–0) spectra along the main horizontal cut is shown in Fig. A.6, which is an extension of the C18 O spectra in Fig. 8. Since theC18O emission extending toward the outskirts becomes asymmetric to the core center, we adopt two slab models for reproducing the spectra at ΔRA ≤ −72′′ (west) and ΔRA ≥ 120′′ (east), respectively, instead of the above onion-shell model. To fit the observations in the outskirts, the C18 O abundance of the slab models are fixed to 10−7 (Sect. 4.4) and is found to be 3.5–5.6 × 1021 cm−2 for the western slab and 0.8–2.6 × 1021 cm−2 for the eastern slab.

thumbnail Fig. A.1

N2H+ (1–0) spectra. The x- and y-axes of the grid are the RA and Dec offsets with respect to the center of L 1512. Each cell shows the observational spectra as black, the modeled spectra as red, and the baselines as blue. The dimension of and VLSR at each cell are denoted in one of the empty cells. Only the central triplet is shown for clarity.

thumbnail Fig. A.2

N2H+ (3–2) spectra. Annotations are the same as in Fig. A.1.

thumbnail Fig. A.3

N2D+ (2–1) spectra. Annotations are the same as in Fig. A.1.

thumbnail Fig. A.4

DCO+ (2–1) spectra. Annotations are the same as in Fig. A.1.

thumbnail Fig. A.5

H2D+ (110 –111) spectra. Annotations are the same as in Fig. A.1.

thumbnail Fig. A.6

C18O (1–0) spectra along the main horizontal cut at ΔDec = −12′′. Annotations are the same as in Fig. A.1 except that the emissions are in the TMB scale.

Appendix B: N2D+–H2 collisional rate coefficients

Rate coefficients for the N2 H+–H2(J = 0) collisional system have been provided by Lique et al. (2015). Hyperfine-structure-resolved rate coefficients, based on a highly accurate ab initio potential energy surface (PES; Spielfiedel et al. 2015), were calculated for temperatures ranging from 5 to 70 K. The new rate coefficients were found to be significantly different from the N2 H+–He rate coefficients previously published (Daniel et al. 2005).

Recent studies (Dumouchel et al. 2012, 2017; Flower & Lique 2015) have shown that isotopic effects in inelastic collisions can be important, especially for H/D substitution. Hence, we decided to compute actual N2 D+–H2(J = 0) rate coefficients.

Within the Born-Oppenheimer approximation, the full electronic ground state potential is identical for N2 H+–H2 and N2 D+–H2(J = 0) systems anddepends only on the mutual distances of the five atoms involved. Then, we used for the scattering calculations the N2 H+–H2 PES of Spielfiedel et al. (2015) and the “adiabatic-hindered-rotor” treatment, which allows para-H2(J = 0) to be treated as if it were spherical. The major difference between the N2 H+–H2 and N2 D+–H2(J = 0) PESs is the position of the center of mass taken for the origin of the Jacobi coordinates. For N2 D+–H2(J = 0) calculations, we considered the effect of the displacement of the center of mass.

Since the nitrogen atoms possess a non-zero nuclear spin (I = 1), the N2 D+ rotational energy levels, such as those of N2 H+, are split in hyperfine levels, which are characterized by the quantum numbers J, F1 and F. In this case, F1 results from the coupling of the rotational level J with I1 (F1 =J +I1, I1 being the nuclear spin of the first nitrogen atom) and F results from the coupling of F1 with I2 (F =F1 +I2, I2 being the nuclear spin of the second nitrogen atom). The D atom also possesses a non-zero nuclear spin. However, in the astronomical observations, the hyperfine structure due to D is not resolved and is then neglected in the calculations.

The hyperfine splitting of the N2 D+ energy levels is very small. Considering that the hyperfine levels are degenerate, we simplified the hyperfine scattering problem using recoupling techniques (Faure & Lique 2012). Then, we performed Close-Coupling calculations (Arthurs & Dalgarno 1960) of the pure rotational excitation cross sections (neglecting the hyperfine structure) using the MOLSCAT program (Hutson & Green 1994). The N2 D+ energy levels were computed using the rotational constants of Dore et al. (2009). Calculations were carried out for total energies up to 500 cm−1. Parameters of the integrator were tested and adjusted to ensure a typical precision to within 0.05 Å2 for the inelastic cross sections. At each energy, channels with J up to 28 were included in the rotational basis to converge the calculations for all the transitions including N2 D+ levels up to J = 7. Using the recoupling technique, the hyperfine state-resolved cross sections were obtained for all hyperfine levels up to J = 7.

From the calculated cross sections, we can obtain the corresponding thermal rate coefficients at temperature T by an average over the collision energy (Ec), i.e., (B.1)

where σαβ is the cross section from initial level α to final level β, μ is the reduced mass of the system, and kB is Boltzmann’s constant.

Using the computational scheme described above, we obtained N2 D+–H2(J = 0) rate coefficients for temperatures up to 70 K. The complete set of (de-)excitation rate coefficients with J, J′≤ 7 will be made available through the LAMDA (Schöier et al. 2005) and BASECOL (Dubernet et al. 2013) databases.

Figure B.1 presents the temperature variation of the N2 H+–H2(J = 0) and N2 D+–H2(J = 0) rate coefficients for selected and transitions.

thumbnail Fig. B.1

Temperature variation of the hyperfine-resolved N2H+–H2(J = 0) (solid lines) and N2D+–H2(J = 0) (dotted lines) rate coefficients for J = 2, FJ′ = 1, F′ and J = 3, FJ′ = 2, F′ transitions.

As is shown, a very good agreement is found between N2 H+–H2(J = 0) and N2 D+–H2(J = 0) rate coefficients. Both sets of data agree within a few percents over all the temperature range. The largest deviations are seen at a low temperature (T < 10 K) that characterizes cold molecular clouds. The N2 H+ over N2 D+ rate coefficients ratio depends on the temperature showing that extrapolation techniques are not suited for the estimation of N2 D+ collisional data. The differences are due to both the center-of-mass shift on the interaction potential and the use of isotopolog specific energies of the levels.

Appendix C: Best-fit physical and abundance profiles

We present the best-fit quantities with their error ranges for each layer in the onion-shell model from Figs. 8 and 10 in Table C.1.

Table C.1

Onion-shell model parameters.

References

  1. Aikawa, Y., Ohashi, N., Inutsuka, S.-i., Herbst, E., & Takakuwa, S. 2001, ApJ, 552, 639 [NASA ADS] [CrossRef] [Google Scholar]
  2. Arthurs, A. M., & Dalgarno, A. 1960, Proc. R. Soc. London Ser. A, 256, 540 [NASA ADS] [CrossRef] [Google Scholar]
  3. Ascenso, J., Lada, C. J., Alves, J. F., Román-Zúñiga, C. G., & Lombardi, M. 2013, A&A, 549, A135 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  4. Bacmann, A., Lefloch, B., Ceccarelli, C., et al. 2002, A&A, 389, L6 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  5. Bergin, E. A., Alves, J., Huard, T., & Lada, C. J. 2002, ApJ, 570, L101 [NASA ADS] [CrossRef] [Google Scholar]
  6. Bernes, C. 1979, A&A, 73, 67 [NASA ADS] [Google Scholar]
  7. Bisschop, S. E., Fraser, H. J., Öberg, K. I., van Dishoeck, E. F., & Schlemmer, S. 2006, A&A, 449, 1297 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  8. Bohlin, R. C., Savage, B. D., & Drake, J. F. 1978, ApJ, 224, 132 [NASA ADS] [CrossRef] [Google Scholar]
  9. Cambrésy, L., Rho, J., Marshall, D. J., & Reach, W. T. 2011, A&A, 527, A141 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  10. Caselli, P., Myers, P. C., & Thaddeus, P. 1995, ApJ, 455, L77 [NASA ADS] [CrossRef] [Google Scholar]
  11. Caselli, P., Walmsley, C. M., Tafalla, M., Dore, L., & Myers, P. C. 1999, ApJ, 523, L165 [NASA ADS] [CrossRef] [Google Scholar]
  12. Caselli, P., Benson, P. J., Myers, P. C., & Tafalla, M. 2002, ApJ, 572, 238 [NASA ADS] [CrossRef] [Google Scholar]
  13. Caselli, P., Vastel, C., Ceccarelli, C., et al. 2008, A&A, 492, 703 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  14. Cox, P., Walmsley, C. M., & Guesten, R. 1989, A&A, 209, 382 [NASA ADS] [Google Scholar]
  15. Crabtree, K. N., Indriolo, N., Kreckel, H., Tom, B. A., & McCall, B. J. 2011, ApJ, 729, 15 [NASA ADS] [CrossRef] [Google Scholar]
  16. Daniel, F., Dubernet, M.-L., Meuwly, M., Cernicharo, J., & Pagani, L. 2005, MNRAS, 363, 1083 [NASA ADS] [CrossRef] [Google Scholar]
  17. Dore, L., Bizzocchi, L., Degli Esposti, C., & Tinti, F. 2009, A&A, 496, 275 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  18. Dubernet, M.-L., Alexander, M. H., Ba, Y. A., et al. 2013, A&A, 553, A50 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  19. Dumouchel, F., Kłos, J., Toboła, R., et al. 2012, J. Chem. Phys., 137, 114306 [NASA ADS] [CrossRef] [Google Scholar]
  20. Dumouchel, F., Lique, F., Spielfiedel, A., & Feautrier, N. 2017, MNRAS, 471, 1849 [NASA ADS] [CrossRef] [Google Scholar]
  21. Falgarone, E., Panis, J.-F., Heithausen, A., et al. 1998, A&A, 331, 669 [NASA ADS] [Google Scholar]
  22. Falgarone, E., Pety, J., & Phillips, T. G. 2001, ApJ, 555, 178 [NASA ADS] [CrossRef] [Google Scholar]
  23. Faure, A., & Lique, F. 2012, MNRAS, 425, 740 [Google Scholar]
  24. Federrath, C., & Klessen, R. S. 2012, ApJ, 761, 156 [NASA ADS] [CrossRef] [Google Scholar]
  25. Flower, D. R., & Lique, F. 2015, MNRAS, 446, 1750 [NASA ADS] [CrossRef] [Google Scholar]
  26. Fontani, F., Caselli, P., Crapsi, A., et al. 2006, A&A, 460, 709 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  27. Foster, J. B., & Goodman, A. A. 2006, ApJ, 636, L105 [NASA ADS] [CrossRef] [Google Scholar]
  28. Frerking, M. A., Langer, W. D., & Wilson, R. W. 1982, ApJ, 262, 590 [NASA ADS] [CrossRef] [Google Scholar]
  29. Fuller, G. A., & Myers, P. C. 1993, ApJ, 418, 273 [NASA ADS] [CrossRef] [Google Scholar]
  30. Giannetti, A., Bovino, S., Caselli, P., et al. 2019, A&A, 621, L7 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  31. Goldsmith, P. F. 2001, ApJ, 557, 736 [NASA ADS] [CrossRef] [MathSciNet] [Google Scholar]
  32. Hartmann, L., Ballesteros-Paredes, J., & Bergin, E. A. 2001, ApJ, 562, 852 [NASA ADS] [CrossRef] [Google Scholar]
  33. Hutson, J. M., & Green, S. 1994, MOLSCAT computer code, version 14 (1994), distributed by Collaborative Computational Project No. 6 of the Engineering and Physical Sciences Research Council (UK) [Google Scholar]
  34. Juvela, M., Pelkonen, V.-M., Padoan, P., & Mattila, K. 2006, A&A, 457, 877 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  35. Khersonkii, V. K., Varshalovich, D. A., & Opendak, M. G. 1987, Sov. Astron., 31, 274 [NASA ADS] [Google Scholar]
  36. Kim, G., Lee, C. W., Gopinathan, M., Jeong, W.-S., & Kim, M.-R. 2016, ApJ, 824, 85 [NASA ADS] [CrossRef] [Google Scholar]
  37. Kong, S., Caselli, P., Tan, J. C., Wakelam, V., & Sipilä, O. 2015, ApJ, 804, 98 [NASA ADS] [CrossRef] [Google Scholar]
  38. Kong, S., Tan,J. C., Caselli, P., et al. 2016, ApJ, 821, 94 [NASA ADS] [CrossRef] [Google Scholar]
  39. Körtgen, B., Bovino, S., Schleicher, D. R. G., Giannetti, A., & Banerjee, R. 2017, MNRAS, 469, 2602 [NASA ADS] [CrossRef] [Google Scholar]
  40. Körtgen, B., Bovino, S., Schleicher, D. R. G., et al. 2018, MNRAS, 478, 95 [NASA ADS] [CrossRef] [Google Scholar]
  41. Kutner, M. L., & Ulich, B. L. 1981, ApJ, 250, 341 [NASA ADS] [CrossRef] [Google Scholar]
  42. Lackington, M., Fuller, G. A., Pineda, J. E., et al. 2016, MNRAS, 455, 806 [NASA ADS] [CrossRef] [Google Scholar]
  43. Launhardt, R., Stutz, A. M., Schmiedeke, A., et al. 2013, A&A, 551, A98 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  44. Lee, C. W., & Myers, P. C. 2011, ApJ, 734, 60 [NASA ADS] [CrossRef] [Google Scholar]
  45. Lee, C. W., Myers, P. C., & Tafalla, M. 2001, ApJS, 136, 703 [NASA ADS] [CrossRef] [Google Scholar]
  46. Lee, J.-E., Evans, II, N. J., Shirley, Y. L., & Tatematsu, K. 2003, ApJ, 583, 789 [NASA ADS] [CrossRef] [Google Scholar]
  47. Lefèvre, C., Pagani, L., Juvela, M., et al. 2014, A&A, 572, A20 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  48. Lefèvre, C., Pagani, L., Min, M., Poteet, C., & Whittet, D. 2016, A&A, 585, L4 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  49. Li, Z.-Y., Shematovich, V. I., Wiebe, D. S., & Shustov, B. M. 2002, ApJ, 569, 792 [NASA ADS] [CrossRef] [Google Scholar]
  50. Lippok, N., Launhardt, R., Semenov, D., et al. 2013, A&A, 560, A41 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  51. Lippok, N., Launhardt, R., Henning, T., et al. 2016, A&A, 592, A61 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  52. Lique, F., Daniel, F., Pagani, L., & Feautrier, N. 2015, MNRAS, 446, 1245 [NASA ADS] [CrossRef] [Google Scholar]
  53. Liu, H.-L., Stutz, A. M., & Yuan, J.-H. 2019, MNRAS, 487, 1259 [NASA ADS] [CrossRef] [Google Scholar]
  54. Lombardi, M., & Alves, J. 2001, A&A, 377, 1023 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  55. Mac Low, M.-M., Klessen, R. S., Burkert, A., & Smith, M. D. 1998, Phys. Rev. Lett., 80, 2754 [NASA ADS] [CrossRef] [Google Scholar]
  56. McKee, C. F., & Ostriker, E. C. 2007, ARA&A, 45, 565 [NASA ADS] [CrossRef] [Google Scholar]
  57. Mouschovias, T. C., Tassis, K., & Kunz, M. W. 2006, ApJ, 646, 1043 [NASA ADS] [CrossRef] [Google Scholar]
  58. Myers, P. C., Linke, R. A., & Benson, P. J. 1983, ApJ, 264, 517 [NASA ADS] [CrossRef] [Google Scholar]
  59. Nguyen, T., Baouche, S., Congiu, E., et al. 2018, A&A, 619, A111 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  60. Nilsson, A., Hjalmarson, Å., Bergman, P., & Millar, T. J. 2000, A&A, 358, 257 [NASA ADS] [Google Scholar]
  61. Padoan, P., & Nordlund, Å. 1999, ApJ, 526, 279 [NASA ADS] [CrossRef] [Google Scholar]
  62. Pagani, L., Salez, M., & Wannier, P. G. 1992a, A&A, 258, 479 [NASA ADS] [Google Scholar]
  63. Pagani, L., Wannier, P. G., Frerking, M. A., et al. 1992b, A&A, 258, 472 [NASA ADS] [Google Scholar]
  64. Pagani, L., Pardo, J.-R., Apponi, A. J., Bacmann, A., & Cabrit, S. 2005, A&A, 429, 181 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  65. Pagani, L., Bacmann, A., Cabrit, S., & Vastel, C. 2007, A&A, 467, 179 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  66. Pagani, L., Vastel, C., Hugo, E., et al. 2009a, A&A, 494, 623 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  67. Pagani, L., Daniel, F., & Dubernet, M.-L. 2009b, A&A, 494, 719 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  68. Pagani, L., Steinacker, J., Bacmann, A., Stutz, A., & Henning, T. 2010a, Science, 329, 1622 [NASA ADS] [CrossRef] [PubMed] [Google Scholar]
  69. Pagani, L., Ristorcelli, I., Boudet, N., et al. 2010b, A&A, 512, A3 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  70. Pagani, L., Roueff, E., & Lesaffre, P. 2011, ApJ, 739, L35 [NASA ADS] [CrossRef] [Google Scholar]
  71. Pagani, L., Bourgoin, A., & Lique, F. 2012, A&A, 548, L4 [CrossRef] [EDP Sciences] [Google Scholar]
  72. Pagani, L., Lesaffre, P., Jorfi, M., et al. 2013, A&A, 551, A38 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  73. Pagani, L., Lefèvre, C., Juvela, M., Pelkonen, V. M., & Schuller, F. 2015, A&A, 574, L5 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  74. Pineda, J. L., Goldsmith, P. F., Chapman, N., et al. 2010, ApJ, 721, 686 [NASA ADS] [CrossRef] [Google Scholar]
  75. Redaelli, E., Bizzocchi, L., Caselli, P., et al. 2019, A&A, 629, A15 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  76. Schöier, F. L., van der Tak, F. F. S., van Dishoeck, E. F., & Black, J. H. 2005, A&A, 432, 369 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  77. Shu, F. H., Adams, F. C., & Lizano, S. 1987, ARA&A, 25, 23 [Google Scholar]
  78. Sohn, J., Lee, C. W., Park, Y.-S., et al. 2007, ApJ, 664, 928 [NASA ADS] [CrossRef] [Google Scholar]
  79. Spielfiedel, A., Senent, M. L., Kalugina, Y., et al. 2015, J. Chem. Phys., 143, 024301 [NASA ADS] [CrossRef] [Google Scholar]
  80. Steinacker, J., Pagani, L., Bacmann, A., & Guieu, S. 2010, A&A, 511, A9 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  81. Steinacker, J., Andersen, M., Thi, W. F., et al. 2015, A&A, 582, A70 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  82. Stutz, A. M., Rieke, G. H., Bieging, J. H., et al. 2009, ApJ, 707, 137 [NASA ADS] [CrossRef] [Google Scholar]
  83. Stutz, A. M., Gonzalez-Lobos, V. I., & Gould, A. 2018, ArXiv e-prints, [arXiv:1807.11496] [Google Scholar]
  84. Tafalla, M., Myers, P. C., Caselli, P., Walmsley, C. M., & Comito, C. 2002, ApJ, 569, 815 [NASA ADS] [CrossRef] [Google Scholar]
  85. Tassis, K., & Mouschovias, T. C. 2004, ApJ, 616, 283 [NASA ADS] [CrossRef] [Google Scholar]
  86. Tazaki, R., Tanaka, H., Okuzumi, S., Kataoka, A., & Nomura, H. 2016, ApJ, 823, 70 [NASA ADS] [CrossRef] [Google Scholar]
  87. van der Tak, F. F. S., Caselli, P., & Ceccarelli, C. 2005, A&A, 439, 195 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  88. Weingartner, J. C., & Draine, B. T. 2001, ApJ, 548, 296 [NASA ADS] [CrossRef] [Google Scholar]
  89. Wilson, T. L.,& Rood, R. 1994, ARA&A, 32, 191 [NASA ADS] [CrossRef] [Google Scholar]

5

It must be noted that despite their claim, including in the title, Redaelli et al. were not the first to detect and explain N2 D+ depletion.

6

The lowest N2 rotational transition J = 2–0 at 357.808 GHz is a very weak electrical quadrupolar transition, and in any case would be blocked by the telluric opacity if it were detectable.

7

The discrepancy between the two depletion factors, 430–870 of 12 CO and 370 of C18 O, is because the standard 12CO abundance of 1–2 × 10−4 and C18 O abundance of 1.7 × 10−7 are not compatible with a standard 12CO/C18O ratio of 500.

8

However, we could not reproduce the Lippok et al. (2013) model results (the emergent intensities being typically a factor of 2 too low based on their values) and the discrepancy has not been identified, but could be linked to a possible confusion between H and H2 densities in the Lippok et al. (2013) paper.

All Tables

Table 1

Observational parameters.

Table C.1

Onion-shell model parameters.

All Figures

thumbnail Fig. 1

Multipointing grids and onion-shell model overlaid with the SCUBA-II 850 μm map. The black and yellow dots in a (ΔRA, ΔDec) = (12′′, 12′′)-spacing grid show the pointings of IRAM 30-m and GBT observations, respectively. The red dots in a (ΔRA, ΔDec) = (15′′, 30′′)-spacing grid show those of JCMT HARP observations. The circles at top and bottom indicate the beam sizes of each spectral observations. The gray concentric circles with center of L1512 represent the onion-shell model with shell widths of 12′′. The 850 μm map is shown in blue-scale, and its white and black contours have the same levels as in Fig. 2g.

In the text
thumbnail Fig. 2

L1512 maps of continuum and line integrated intensity. The CFHT NIR maps at (a) J band, (b) H band, and (c) Ks band. Spitzer MIR maps at (d) IRAC1 band, (e) IRAC2 band, and (f) IRAC4 band with its contours at 5.15, 5.25, 5.35, 5.45, 5.55 MJy sr−1. (g) JCMT SCUBA-II map at 850 μm with its contours at 0, 0.05, 0.15, 0.25 mJy arcsec−2. Integrated intensity maps of (h) N2H+ J = 1–0 and (i) C18O J = 1–0 calculated within VLSR = [−3 km s−1, 16 km s−1] and [6.0 km s−1, 7.5 km s−1], respectively. The panels h and i are overlaid with contours from the 850 μm image. The central cross in each panel indicates the center of L1512 determined by the IRAC4 map. The scale bars of 0.05 pc and millimeter-wavelength beam sizes are denoted in the bottom right and bottom left corners, respectively.

In the text
thumbnail Fig. 3

Spectral observations along the main horizontal cuts at ΔDec = −12′′ (IRAM/GBT) or −16′′ (JCMT) compared to our best-fit radiative transfer model. The blue spectra show the observational data, and the red spectra show the models. Each column corresponds to different horizontal offsets from the center of L 1512 according to Fig. 1. Each row shows a spectral line, except that the N2 H+ (1–0), N2D+ (1–0) and (2–1) lines are split into three rows corresponding to their different F1 -transition groups. For C18O (2–1), our data are supplemented with OTF spectra from previous works at large offsets. The green dashed lines indicate the three σ noise level.

In the text
thumbnail Fig. 4

Spectral observations along the main vertical cuts at ΔRA = 0′′ (IRAM/GBT) or −2.2′′ (JCMT) compared to our best-fit radiative transfer model. For N2H+ (1–0) and C18O (2–1), our data are supplemented with OTF spectra from previous works at large offsets. Other annotations are the same as in Fig. 3.

In the text
thumbnail Fig. 5

Velocity structure. Mean velocity maps of (a) the F1F = 2,3–1,2 component of N2H+ J = 1–0 and (b) C18O J = 1–0 calculated within VLSR = [6.75 km s−1, 7.50 km s−1] and [6.0 km s−1, 7.5 km s−1], respectively. Contours from the 850 μm image are overlaid. The rest of the symbols in the panels a and b have the same meaning as in Fig. 2. (c) Model vertical rotation velocities, Vrot, against radius, Δr, adopted in the radiative transfer model. (d) Central velocities and their error bars of the spectra along the main vertical cuts. Two y-axes are vertical offsets to the center and declination, respectively.

In the text
thumbnail Fig. 6

Visualextinction maps and profiles. (a) Envelope AV map derived from J, H, and Ks bands with the RV = 3.1 model. (b) Core AV map derived from H, Ks, and IRAC2 bands with the RV = 5.5 B model. (c) Combined AV map merged from the previous maps with boundaries at AV = 2.5 mag. Each map is convolved with a Gaussian beam size of 30′′. (d) Comparison between AV profiles inside the strips of panel a and b. The AV values at each pixel are denoted in light blue and pink dots, while their average errors are indicated in the bottom left corner. The AV profiles are averaged in azimuth with radial 10′′ bin and shown in blue triangles and red squares with their error bars. The green dashed lines in the four planes indicate the AV threshold of 2.5 mag.

In the text
thumbnail Fig. 7

Plummer-like density profile fittings. The black dots are the AV values at each pixel in the horizontal strip shown in Fig. 6c but these are separated between the east and west sides. The orange squares with error bars show the averaged AV profiles, which are azimuthally binned with radial 10′′ bin. The blue circles show the best-fitted AV models for the east and west sides. The AV models are derived from spherically symmetrical Plummer-like density profiles. The best-fitted characteristic radii, R0, and the fitting ranges, Rfitting of the AV models are denoted in blue solid and dash-dotted lines, respectively. The coverage of our spectral observations are denoted in black dashed lines.

In the text
thumbnail Fig. 8

Physical and abundance profiles. The Plummer-like density () profiles are derived from the combined extinction map (Fig. 6c), where the best-fit parameters are annotated. The profiles of kinetic temperature (Tkin), turbulent velocity (Vturb), abundances (Xspecies), and the abundance ratio of N2D+/N2H+ are the best-fit results from the non-LTE radiative transfer calculations.

In the text
thumbnail Fig. 9

Chemical modeling of the abundance ratio of N2D+/N2H+ and the abundances of ortho-H2D+ and DCO+ for each layer. Left and right columns: chemical solutions. Middle column: observationally derived profiles from Fig. 8. Each color represents a layer. In the chemical solution plots, the horizontal solid lines are the observationally derived X(N2D+)/X(N2H+), X(ortho-H2D+), and X(DCO+). The growing curves are the chemical solutions calculated with an average grain size (agr) of 0.1 μm and a cosmicray ionization rate (ζ) of 10−17 s−1. The two black vertical lines in each panel indicate a time range for which the model values cross the observations within their error bars and where the model curves become thicker. Thick dashed lines denote the model range that is both lower than the observed upper limit for that layer and within the global time range of the model. The red dashed vertical lines indicate the lower limit on the lifetime scale of L1512.

In the text
thumbnail Fig. 10

12CO, N2, and C18 O abundance profiles. The 12CO and N2 abundance profiles are the best-fit results from the chemical modeling. The C18 O abundance profile is obtained by assuming a 12CO/C18O abundance ratio of 500 or a constant value of 1.0 × 10−7 extending to 25.3 kAU.

In the text
thumbnail Fig. 11

H isotopolog abundances derived from the chemical model constrained by the observations (in particular of ortho-H2D+, which explains why we separately show the two spin states for this species for the sake of clarity).

In the text
thumbnail Fig. A.1

N2H+ (1–0) spectra. The x- and y-axes of the grid are the RA and Dec offsets with respect to the center of L 1512. Each cell shows the observational spectra as black, the modeled spectra as red, and the baselines as blue. The dimension of and VLSR at each cell are denoted in one of the empty cells. Only the central triplet is shown for clarity.

In the text
thumbnail Fig. A.2

N2H+ (3–2) spectra. Annotations are the same as in Fig. A.1.

In the text
thumbnail Fig. A.3

N2D+ (2–1) spectra. Annotations are the same as in Fig. A.1.

In the text
thumbnail Fig. A.4

DCO+ (2–1) spectra. Annotations are the same as in Fig. A.1.

In the text
thumbnail Fig. A.5

H2D+ (110 –111) spectra. Annotations are the same as in Fig. A.1.

In the text
thumbnail Fig. A.6

C18O (1–0) spectra along the main horizontal cut at ΔDec = −12′′. Annotations are the same as in Fig. A.1 except that the emissions are in the TMB scale.

In the text
thumbnail Fig. B.1

Temperature variation of the hyperfine-resolved N2H+–H2(J = 0) (solid lines) and N2D+–H2(J = 0) (dotted lines) rate coefficients for J = 2, FJ′ = 1, F′ and J = 3, FJ′ = 2, F′ transitions.

In the text

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

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

Initial download of the metrics may take a while.